Newer
Older
#include <arbor/morph/primitives.hpp>
#include <arbor/morph/segment_tree.hpp>
#include <arbor/version.hpp>
#include <arborio/neurolucida.hpp>
#ifdef ARB_NEUROML_ENABLED
#include <arborio/neuroml.hpp>
#include "proxy.hpp"
static inline bool cable_lt(const arb::mcable& a, const arb::mcable& b) {
return std::tuple(a.branch, a.prox_pos, a.dist_pos)<std::tuple(b.branch, b.prox_pos, b.dist_pos);
}
void check_trailing(std::istream& in, std::string fname) {
if (!(in >> std::ws).eof()) {
throw pyarb_error(util::pprintf("Trailing data found at end of file '{}'", fname));
}
}
void register_morphology(py::module& m) {
using namespace py::literals;
// primitives: points, segments, locations, cables... etc.
//
m.attr("mnpos") = arb::mnpos;
// arb::mlocation
py::class_<arb::mlocation> location(m, "location",
"A location on a cable cell.");
location
[](arb::msize_t branch, double pos) {
const arb::mlocation mloc{branch, pos};
pyarb::assert_throw(arb::test_invariants(mloc), "invalid location");
return mloc;
}),
"Construct a location specification holding:\n"
" branch: The id of 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.")
"The relative position on the branch (∈ [0.,1.], where 0. means proximal and 1. distal).")
[](arb::mlocation l) { return util::pprintf("(location {} {})", l.branch, l.pos); })
[](arb::mlocation l) { return util::pprintf("(location {} {})", l.branch, l.pos); });
py::class_<arb::mpoint> mpoint(m, "mpoint");
.def(py::init<double, double, double, double>(),
"x"_a, "y"_a, "z"_a, "radius"_a,
"Create an mpoint object from parameters x, y, z, and radius, specified in µm.")
.def(py::init([](py::tuple t) {
if (py::len(t)!=4) throw std::runtime_error("tuple length != 4");
return arb::mpoint{t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>(), t[3].cast<double>()}; }),
"Create an mpoint object from a tuple (x, y, z, radius), specified 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 centred at coordinates [μm].")
[](const arb::mpoint& p) {
return util::pprintf("<arbor.mpoint: x {}, y {}, z {}, radius {}>", p.x, p.y, p.z, p.radius);
})
[](const arb::mpoint& p) {return util::pprintf("{}", p);});
py::implicitly_convertible<py::tuple, arb::mpoint>();
// arb::msegment
py::class_<arb::msegment> msegment(m, "msegment");
msegment
.def_readonly("prox", &arb::msegment::prox, "the location and radius of the proximal end.")
.def_readonly("dist", &arb::msegment::dist, "the location and radius of the distal end.")
.def_readonly("tag", &arb::msegment::tag, "tag meta-data.");
py::class_<arb::mcable> cable(m, "cable");
[](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("__str__", [](const arb::mcable& c) { return util::pprintf("{}", c); })
.def("__repr__", [](const arb::mcable& c) { return util::pprintf("{}", c); });
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
// arb::isometry
py::class_<arb::isometry> isometry(m, "isometry");
isometry
.def(py::init<>(), "Construct a trivial isometry.")
.def("__call__", [](arb::isometry& iso, arb::mpoint& p) {
return iso.apply(p);
},
"Apply isometry to mpoint argument.")
.def("__call__", [](arb::isometry& iso, py::tuple t) {
int len = py::len(t);
if (len<3) throw std::runtime_error("tuple length < 3");
arb::mpoint p{t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>(), 0.};
p = iso.apply(p);
py::tuple result(len);
result[0] = p.x;
result[1] = p.y;
result[2] = p.z;
for (int i = 3; i<len; ++i) {
result[i] = t[i];
}
return result;
},
"Apply isometry to first three components of tuple argument.")
.def(py::self*py::self)
.def_static("translate",
[](double x, double y, double z) { return arb::isometry::translate(x, y, z); },
"x"_a, "y"_a, "z"_a,
"Construct a translation isometry from displacements x, y, and z.")
.def_static("translate",
[](py::tuple t) {
if (py::len(t)!=3) throw std::runtime_error("tuple length != 3");
return arb::isometry::translate(t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>());
},
"Construct a translation isometry from the first three components of a tuple.")
.def_static("translate",
[](arb::mpoint p) { return arb::isometry::translate(p.x, p.y, p.z); },
"Construct a translation isometry from the x, y, and z components of an mpoint.")
.def_static("rotate",
[](double theta, double x, double y, double z) { return arb::isometry::rotate(theta, x, y, z); },
"theta"_a, "x"_a, "y"_a, "z"_a,
"Construct a rotation isometry of angle theta about the axis in direction (x, y, z).")
.def_static("rotate",
[](double theta, py::tuple t) {
if (py::len(t)!=3) throw std::runtime_error("tuple length != 3");
return arb::isometry::rotate(theta, t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>());
},
"theta"_a, "axis"_a,
"Construct a rotation isometry of angle theta about the given axis in the direction described by a tuple.");
// arb::place_pwlin
py::class_<arb::place_pwlin> place(m, "place_pwlin");
place
.def(py::init<const arb::morphology&, const arb::isometry&>(),
"morphology"_a, "isometry"_a=arb::isometry{},
"Construct a piecewise-linear placement object from the given morphology and optional isometry.")
.def("at", &arb::place_pwlin::at, "location"_a,
"Return an interpolated mpoint corresponding to the location argument.")
.def("all_at", &arb::place_pwlin::all_at, "location"_a,
"Return list of all possible interpolated mpoints corresponding to the location argument.")
.def("segments",
[](const arb::place_pwlin& self, std::vector<arb::mcable> cables) {
std::sort(cables.begin(), cables.end(), cable_lt);
return self.segments(cables);
},
"Return minimal list of full or partial msegments whose union is coterminous "
"with the extent of the given list of cables.")
.def("all_segments",
[](const arb::place_pwlin& self, std::vector<arb::mcable> cables) {
std::sort(cables.begin(), cables.end(), cable_lt);
return self.all_segments(cables);
},
"Return maximal list of non-overlapping full or partial msegments whose union is coterminous "
"with the extent of the given list of cables.");
// Higher-level data structures (segment_tree, morphology)
// arb::segment_tree
py::class_<arb::segment_tree> segment_tree(m, "segment_tree");
.def("reserve", &arb::segment_tree::reserve)
.def("append", [](arb::segment_tree& t, arb::msize_t parent, arb::mpoint prox, arb::mpoint dist, int tag) {
return t.append(parent, prox, dist, tag);
},
"parent"_a, "prox"_a, "dist"_a, "tag"_a,
"Append a segment to the tree.")
.def("append", [](arb::segment_tree& t, arb::msize_t parent, arb::mpoint dist, int tag) {
return t.append(parent, dist, tag);
},
"parent"_a, "dist"_a, "tag"_a,
"Append a segment to the tree.")
.def("append",
[](arb::segment_tree& t, arb::msize_t p, double x, double y, double z, double radius, int tag) {
return t.append(p, arb::mpoint{x,y,z,radius}, tag);
},
"parent"_a, "x"_a, "y"_a, "z"_a, "radius"_a, "tag"_a,
"Append a segment to the tree, using the distal location of the parent segment as the proximal end.")
.def_property_readonly("empty", [](const arb::segment_tree& st){return st.empty();},
"Indicates whether the tree is empty (i.e. whether it has size 0)")
.def_property_readonly("size", [](const arb::segment_tree& st){return st.size();},
.def_property_readonly("parents", [](const arb::segment_tree& st){return st.parents();},
.def_property_readonly("segments", [](const arb::segment_tree& st){return st.segments();},
.def("__str__", [](const arb::segment_tree& s) {
return util::pprintf("<arbor.segment_tree:\n{}>", s);});
// Function that creates a morphology from an swc file.
// Wraps calls to C++ functions arborio::parse_swc() and arborio::load_swc_arbor().
[](std::string fname) {
std::ifstream fid{fname};
if (!fid.good()) {
throw pyarb_error(util::pprintf("can't open file '{}'", fname));
}
try {
auto data = arborio::parse_swc(fid);
check_trailing(fid, fname);
return arborio::load_swc_arbor(data);
}
catch (arborio::swc_error& e) {
// Try to produce helpful error messages for SWC parsing errors.
throw pyarb_error(util::pprintf("error parsing {}: {}", fname, e.what()));
}
},
"filename"_a,
"Generate a morphology from an SWC file following the rules prescribed by Arbor.\n"
"Specifically:\n"
"* Single-segment somas are disallowed. These are usually interpreted as spherical somas\n"
" and are a special case. This behavior is not allowed using this SWC loader.\n"
"* There are no special rules related to somata. They can be one or multiple branches\n"
" and other segments can connect anywhere along them.\n"
"* A segment is always created between a sample and its parent, meaning there\n"
" are no gaps in the resulting morphology.");
m.def("load_swc_allen",
[](std::string fname, bool no_gaps=false) {
std::ifstream fid{fname};
if (!fid.good()) {
throw pyarb_error(util::pprintf("can't open file '{}'", fname));
}
try {
auto data = arborio::parse_swc(fid);
check_trailing(fid, fname);
return arborio::load_swc_allen(data, no_gaps);
}
catch (arborio::swc_error& e) {
// Try to produce helpful error messages for SWC parsing errors.
throw pyarb_error(
util::pprintf("Allen SWC: error parsing {}: {}", fname, e.what()));
}
},
"filename"_a, "no_gaps"_a=false,
"Generate a morphology from an SWC file following the rules prescribed by AllenDB\n"
" and Sonata. Specifically:\n"
"* The first sample (the root) is treated as the centre of the soma.\n"
"* The first morphology is translated such that the soma is centred at (0,0,0).\n"
"* The first sample has tag 1 (soma).\n"
"* All other samples have tags 2, 3 or 4 (axon, apic and dend respectively)\n"
"SONATA prescribes that there should be no gaps, however the models in AllenDB\n"
"have gaps between the start of sections and the soma. The flag no_gaps can be\n"
"used to enforce this requirement.\n"
"\n"
"Arbor does not support modelling the soma as a sphere, so a cylinder with length\n"
"equal to the soma diameter is used. The cylinder is centred on the origin, and\n"
"aligned along the z axis.\n"
"Axons and apical dendrites are attached to the proximal end of the cylinder, and\n"
"dendrites to the distal end, with a gap between the start of each branch and the\n"
"end of the soma cylinder to which it is attached.");
m.def("load_swc_neuron",
[](std::string fname) {
std::ifstream fid{fname};
if (!fid.good()) {
throw pyarb_error(util::pprintf("can't open file '{}'", fname));
}
try {
auto data = arborio::parse_swc(fid);
check_trailing(fid, fname);
return arborio::load_swc_neuron(data);
// Try to produce helpful error messages for SWC parsing errors.
throw pyarb_error(
util::pprintf("NEURON SWC: error parsing {}: {}", fname, e.what()));
"Generate a morphology from an SWC file following the rules prescribed by NEURON.\n"
" Specifically:\n"
"* The first sample must be a soma sample.\n"
"* The soma is represented by a series of n≥1 unbranched, serially listed samples.\n"
"* The soma is constructed as a single cylinder with diameter equal to the piecewise\n"
" average diameter of all the segments forming the soma.\n"
"* A single-sample soma at is constructed as a cylinder with length=diameter.\n"
"* If a non-soma sample is to have a soma sample as its parent, it must have the\n"
" most distal sample of the soma as the parent.\n"
"* Every non-soma sample that has a soma sample as its parent, attaches to the\n"
" created soma cylinder at its midpoint.\n"
"* If a non-soma sample has a soma sample as its parent, no segment is created\n"
" between the sample and its parent, instead that sample is the proximal point of\n"
" a new segment, and there is a gap in the morphology (represented electrically as a\n"
" zero-resistance wire)\n"
"* To create a segment with a certain tag, that is to be attached to the soma,\n"
" we need at least 2 samples with that tag."
);
py::class_<arb::morphology> morph(m, "morphology");
[](arb::segment_tree t){
return arb::morphology(std::move(t));
}))
// morphology's interface is read-only by design, so most of it can
// be implemented as read-only properties.
.def_property_readonly("empty",
[](const arb::morphology& m){return m.empty();},
.def_property_readonly("num_branches",
[](const arb::morphology& m){return m.num_branches();},
"The number of branches in the morphology.")
.def("branch_parent", &arb::morphology::branch_parent,
"i"_a, "The parent branch of branch i.")
.def("branch_children", &arb::morphology::branch_children,
"i"_a, "The child branches of branch i.")
.def("branch_segments",
[](const arb::morphology& m, arb::msize_t i) {
return m.branch_segments(i);
"i"_a, "A list of the segments in branch i, ordered from proximal to distal ends of the branch.")
.def("__str__",
[](const arb::morphology& m) {
return util::pprintf("<arbor.morphology:\n{}>", m);
});
// Neurolucida ASCII, or .asc, file format support.
py::class_<arborio::asc_morphology> asc_morphology(m, "asc_morphology",
"The morphology and label dictionary meta-data loaded from a Neurolucida ASCII (.asc) file.");
asc_morphology
.def_readonly("morphology",
&arborio::asc_morphology::morphology,
"The cable cell morphology")
.def_property_readonly("labels",
[](const arborio::asc_morphology& m) {return label_dict_proxy(m.labels);},
"The four canonical regions are labeled 'soma', 'axon', 'dend' and 'apic'.");
m.def("load_asc",
[](std::string fname) {
try {
return arborio::load_asc(fname);
}
catch (std::exception& e) {
// Try to produce helpful error messages for SWC parsing errors.
throw pyarb_error(util::pprintf("error loading neurolucida asc file: {}", e.what()));
}
},
"filename"_a, "Load a morphology and meta data from a Neurolucida ASCII .asc file.");
#ifdef ARB_NEUROML_ENABLED
// arborio::morphology_data
py::class_<arborio::nml_morphology_data> nml_morph_data(m, "neuroml_morph_data");
nml_morph_data
.def_readonly("cell_id",
&arborio::nml_morphology_data::cell_id,
"Cell id, or empty if morphology was taken from a top-level <morphology> element.")
.def_readonly("id",
&arborio::nml_morphology_data::id,
"Morphology id.")
.def_readonly("morphology",
&arborio::nml_morphology_data::morphology,
"Morphology constructed from a signle NeuroML <morphology> element.")
.def("segments",
[](const arborio::nml_morphology_data& md) {return label_dict_proxy(md.segments);},
"Label dictionary containing one region expression for each segment id.")
.def("named_segments",
[](const arborio::nml_morphology_data& md) {return label_dict_proxy(md.named_segments);},
"Label dictionary containing one region expression for each name applied to one or more segments.")
.def("groups",
[](const arborio::nml_morphology_data& md) {return label_dict_proxy(md.groups);},
"Label dictionary containing one region expression for each segmentGroup id.")
.def_readonly("group_segments",
&arborio::nml_morphology_data::group_segments,
"Map from segmentGroup ids to their corresponding segment ids.");
// arborio::neuroml
py::class_<arborio::neuroml> neuroml(m, "neuroml");
neuroml
// constructors
.def(py::init(
[](std::string fname){
std::ifstream fid{fname};
if (!fid.good()) {
throw pyarb_error(util::pprintf("can't open file '{}'", fname));
}
try {
std::string string_data((std::istreambuf_iterator<char>(fid)),
std::istreambuf_iterator<char>());
return arborio::neuroml(string_data);
}
catch (arborio::neuroml_exception& e) {
// Try to produce helpful error messages for NeuroML parsing errors.
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
throw pyarb_error(
util::pprintf("NeuroML error processing file {}: ", fname, e.what()));
}
}))
.def("cell_ids",
[](const arborio::neuroml& nml) {
try {
return nml.cell_ids();
}
catch (arborio::neuroml_exception& e) {
throw util::pprintf("NeuroML error: {}", e.what());
}
},"Query top-level cells.")
.def("morphology_ids",
[](const arborio::neuroml& nml) {
try {
return nml.morphology_ids();
}
catch (arborio::neuroml_exception& e) {
throw util::pprintf("NeuroML error: {}", e.what());
}
},"Query top-level standalone morphologies.")
.def("morphology",
[](const arborio::neuroml& nml, const std::string& morph_id) {
try {
return nml.morphology(morph_id);
}
catch (arborio::neuroml_exception& e) {
throw util::pprintf("NeuroML error: {}", e.what());
}
},"morph_id"_a,"Retrieve top-level nml_morph_data associated with morph_id.")
.def("cell_morphology",
[](const arborio::neuroml& nml, const std::string& cell_id) {
try {
return nml.cell_morphology(cell_id);
}
catch (arborio::neuroml_exception& e) {
throw util::pprintf("NeuroML error: {}", e.what());
}
},"morph_id"_a,"Retrieve nml_morph_data associated with cell_id.");
#endif