diff --git a/.gitignore b/.gitignore index ba311b4ae76ec977dbd8b100924a64bcbda71658..ee3b410fd70106ea214de3851302036bda6cda97 100644 --- a/.gitignore +++ b/.gitignore @@ -87,3 +87,4 @@ dist # generated image files by Python examples python/example/*.svg +python/example/*.csv diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp index 71a60778434e19bcb3b35cd3f8c052098bd9580d..e449df768a35107c14caa9487eeb3128835751ad 100644 --- a/arbor/include/arbor/morph/locset.hpp +++ b/arbor/include/arbor/morph/locset.hpp @@ -136,6 +136,12 @@ locset most_distal(region reg); // Most proximal points of a region. locset most_proximal(region reg); +// Translate locations in locset distance μm in the distal direction +locset distal_translate(locset ls, double distance); + +// Translate locations in locset distance μm in the proximal direction +locset proximal_translate(locset ls, double distance); + // Boundary points of a region. locset boundary(region reg); diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp index fc2188e13908efcf75f1da55980d1178cb8cf5c1..0cb46fe2d95ffe8a088d53f57bb9ea14e6e3a5b5 100644 --- a/arbor/morph/locset.cpp +++ b/arbor/morph/locset.cpp @@ -1,6 +1,7 @@ #include <algorithm> #include <iostream> #include <numeric> +#include <stack> #include <arbor/math.hpp> #include <arbor/morph/locset.hpp> @@ -116,6 +117,138 @@ std::ostream& operator<<(std::ostream& o, const terminal_& x) { return o << "(terminal)"; } +// Translate locations in locset distance μm in the proximal direction +struct proximal_translate_: locset_tag { + proximal_translate_(const locset& ls, double distance): start(ls), distance(distance) {} + locset start; + double distance; +}; + +mlocation_list thingify_(const proximal_translate_& dt, const mprovider& p) { + const auto& m = p.morphology(); + const auto& e = p.embedding(); + + std::vector<mlocation> L; + + const auto start = thingify(dt.start, p); + const auto distance = dt.distance; + + for (auto loc: start) { + auto distance_remaining = distance; + + while (loc.branch != mnpos) { + const auto branch = loc.branch; + const auto branch_length = e.branch_length(branch); + const auto dist_pos = loc.pos; + const auto prox_pos = dist_pos - distance_remaining / branch_length; + + if (prox_pos>=0) { + // The target is inside this branch. + L.push_back({branch, prox_pos}); + break; + } + if (auto parent = m.branch_parent(branch); parent==mnpos) { + // The root has been reached; return the start of the branch attached to root. + L.push_back({branch, 0.}); + break; + } + else { + loc = {parent, 1.}; + } + distance_remaining -= dist_pos*branch_length; + } + } + + return L; +} + +locset proximal_translate(locset ls, double distance) { + return locset(proximal_translate_{ls, distance}); +} + +std::ostream& operator<<(std::ostream& o, const proximal_translate_& l) { + return o << "(proximal-translate " << l.start << " " << l.distance << ")"; +} + +// Translate locations in locset distance μm in the distal direction +struct distal_translate_: locset_tag { + distal_translate_(const locset& ls, double distance): start(ls), distance(distance) {} + locset start; + double distance; +}; + +locset distal_translate(locset ls, double distance) { + return locset(distal_translate_{ls, distance}); +} + +mlocation_list thingify_(const distal_translate_& dt, const mprovider& p) { + const auto& m = p.morphology(); + const auto& e = p.embedding(); + + std::vector<mlocation> L; + + auto start = thingify(dt.start, p); + auto distance = dt.distance; + + struct branch_interval { + msize_t bid; + double distance; + }; + + for (auto c: start) { + std::stack<branch_interval> branches_reached; + bool first_branch = true; + + // If starting at the end of a branch, start traversal with its children + if (c.pos < 1) { + branches_reached.push({c.branch, distance}); + } else { + first_branch = false; + for (auto child: m.branch_children(c.branch)) { + branches_reached.push({child, distance}); + } + } + + while (!branches_reached.empty()) { + auto bi = branches_reached.top(); + branches_reached.pop(); + + auto branch = bi.bid; + auto rem_dist = bi.distance; + + auto branch_length = e.branch_length(branch); + auto prox_pos = first_branch*c.pos; + auto dist_pos = rem_dist / branch_length + prox_pos; + + if (dist_pos <= 1) { + // The location lies inside this branch. + L.push_back({branch, dist_pos}); + } + else if (m.branch_children(branch).empty()) { + // The location is more proximal than this branch, but this + // branch is terminal add the terminal location. + L.push_back({branch, 1}); + } + else { + // The location is more distal than this branch: add child + // branches to continue the search. + rem_dist = rem_dist - (1 - prox_pos)*branch_length; + for (auto child: m.branch_children(branch)) { + branches_reached.push({child, rem_dist}); + } + } + first_branch = false; + } + } + + std::sort(L.begin(), L.end()); + return L; +} + +std::ostream& operator<<(std::ostream& o, const distal_translate_& l) { + return o << "(distal-translate " << l.start << " " << l.distance << ")"; +} + // Root location (most proximal point). struct root_: locset_tag {}; diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp index b36c994d12bb826243c5133503a8f4c9ed0ff2c9..038b1364a3cdc2a395ae42597c78eebceaa28d6e 100644 --- a/arbor/morph/region.cpp +++ b/arbor/morph/region.cpp @@ -305,8 +305,8 @@ mextent thingify_(const proximal_interval_& reg, const mprovider& p) { std::vector<mcable> L; - auto start = thingify(reg.end, p); - auto distance = reg.distance; + const auto start = thingify(reg.end, p); + const auto distance = reg.distance; for (auto c: start) { auto branch = c.branch; diff --git a/arborio/label_parse.cpp b/arborio/label_parse.cpp index df45c15f0395b8519647f82987d83d1b17a6ff7a..4a2a3ce229ff8cc77b863525fe5855da89813d48 100644 --- a/arborio/label_parse.cpp +++ b/arborio/label_parse.cpp @@ -86,6 +86,10 @@ std::unordered_multimap<std::string, evaluator> eval_map { "'distal' with 1 argument: (reg:region)")}, {"proximal", make_call<arb::region>(arb::ls::most_proximal, "'proximal' with 1 argument: (reg:region)")}, + {"distal-translate", make_call<arb::locset, double>(arb::ls::distal_translate, + "'distal-translate' with 2 arguments: (ls:locset distance:real)")}, + {"proximal-translate", make_call<arb::locset, double>(arb::ls::proximal_translate, + "'proximal-translate' with 2 arguments: (ls:locset distance:real)")}, {"uniform", make_call<arb::region, int, int, int>(arb::ls::uniform, "'uniform' with 4 arguments: (reg:region, first:int, last:int, seed:int)")}, {"on-branches", make_call<double>(arb::ls::on_branches, diff --git a/doc/concepts/labels.rst b/doc/concepts/labels.rst index 4c1ee631de222cf4baafc7135e775a09a18537b6..845674e48b18c186ad4b11f6c1310ffce1631e1f 100644 --- a/doc/concepts/labels.rst +++ b/doc/concepts/labels.rst @@ -232,7 +232,7 @@ Locset expressions of branch length, so for example, on a branch of length 100 μm ``pos=0.2`` corresponds to 20 μm from the proximal end, or 80 μm from the distal end. - .. figure:: ../gen-images/location_label.svg + .. figure:: ../gen-images/location_05_label.svg :width: 300 :align: center @@ -293,6 +293,39 @@ Locset expressions On the left is the region with radius between 0.3 μm and 0.5 μm. The right shows the proximal set of this region. +.. label:: (proximal-translate ls:locset distance:real) + + The set of locations that correspond to moving each location in the ``ls`` in the proximal direction + ``distance`` μm. The locations in the output have a one to one correspondance with those in ``ls``. + + .. figure:: ../gen-images/proximal_translate_label.svg + :width: 600 + :align: center + + The proximal translation of the terminal locations (left) a distance of 10 μm using the + expression ``(proximal-translate (terminal) 10)``. + +.. label:: (distal-translate ls:locset distance:real) + + The set of locations that correspond to translating each location in ``ls`` in the distal direction + ``distance`` μm or to a terminal location, whichever is closest. + + An input location will generate multiple output locations when it is translated + past a fork point, with a new location for each child branch (see the example + below). For this reason there is not a one-to-one correspondance between locations + in the input and output sets, so the results are sorted and duplicates are removed. + + + .. figure:: ../gen-images/distal_translate_label.svg + :width: 600 + :align: center + + Two distal translations of the midpoint of branch 0 (left). + The first translation of 5 μm, ``(distal-translate (location 0 0.5) 5)``, generates a + single location on the same branch (center). + The second translation of 15 μm ``(distal-translate (location 0 0.5) 15)`` extends beyond + the end of branch 0, generating an additional location at each fork point (right). + .. label:: (locset name:string) Refer to a locset by its label. For example, ``(locset "synapse-sites")`` could be used in an expression to refer diff --git a/doc/scripts/gen-labels.py b/doc/scripts/gen-labels.py index bcb7c7100e82cf297240487f6b1b5af1c5c10dab..cf8466f64e012d729ba474345a92ba43cfde8939 100644 --- a/doc/scripts/gen-labels.py +++ b/doc/scripts/gen-labels.py @@ -159,7 +159,7 @@ tree.append(mnpos, mpoint(-3.0, 0.0, 0.0, 1.5), mpoint(-5.5,-0.2, 0.0, 0.5), tree.append(8, mpoint(-14.5,-0.1, 0.0, 0.5), tag=2) ysoma_morph3 = arbor.morphology(tree) -fn = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__), "../concepts/example.swc")) +fn = os.path.realpath(os.path.join(os.getcwd(), os.path.dirname(__file__), "../fileformat/example.swc")) swc_morph = arbor.load_swc_arbor(fn) regions = { @@ -199,6 +199,7 @@ locsets = { 'term': '(terminal)', 'rand_dend': '(uniform (region "dend") 0 50 0)', 'loc15': '(location 1 0.5)', + 'loc05': '(location 0 0.5)', 'uniform0': '(uniform (tag 3) 0 9 0)', 'uniform1': '(uniform (tag 3) 0 9 1)', 'branchmid': '(on-branches 0.5)', @@ -208,6 +209,9 @@ locsets = { 'proxint_in': '(sum (location 1 0.8) (location 2 0.3))', 'loctest' : '(distal (complete (join (branch 1) (branch 0))))', 'restrict': '(restrict (terminal) (tag 3))', + 'proximal_translate': '(proximal-translate (terminal) 10)', + 'distal_translate_single': '(distal-translate (location 0 0.5) 5)', + 'distal_translate_multi': '(distal-translate (location 0 0.5) 15)', } labels = {**regions, **locsets} diff --git a/doc/scripts/inputs.py b/doc/scripts/inputs.py index 3befe7ad031e1d13d128b41665f271ad3c235d7f..da750d9ee69598c706fc89fd20a6513ea81c3911 100644 --- a/doc/scripts/inputs.py +++ b/doc/scripts/inputs.py @@ -111,6 +111,7 @@ ls_root = {'type': 'locset', 'value': [(0, 0.0)]} ls_term = {'type': 'locset', 'value': [(1, 1.0), (3, 1.0), (4, 1.0), (5, 1.0)]} ls_rand_dend = {'type': 'locset', 'value': [(0, 0.5547193370156588), (0, 0.5841758202819731), (0, 0.607192003545501), (0, 0.6181091003428546), (0, 0.6190845627201184), (0, 0.7027325639263277), (0, 0.7616129092226993), (0, 0.9645150497869694), (1, 0.15382287505908834), (1, 0.2594719824047551), (1, 0.28087652335178354), (1, 0.3729681478609085), (1, 0.3959560134241004), (1, 0.4629424550242548), (1, 0.47346867377446744), (1, 0.5493486883630476), (1, 0.6227685370674116), (1, 0.6362196581003494), (1, 0.6646511214508091), (1, 0.7157318936458146), (1, 0.7464198558822775), (1, 0.77074507802833), (1, 0.7860238136304932), (1, 0.8988928261704698), (1, 0.9581259332943499), (2, 0.12773985425987294), (2, 0.3365926476076694), (2, 0.44454300804769703), (2, 0.5409466695719178), (2, 0.5767511435223905), (2, 0.6340206909931745), (2, 0.6354772583375223), (2, 0.6807941995943213), (2, 0.774655947503608), (3, 0.05020708596877571), (3, 0.25581431877212274), (3, 0.2958305460715556), (3, 0.296698184761692), (3, 0.509669134988683), (3, 0.7662305637426007), (3, 0.8565839889923518), (3, 0.8889077221517746), (4, 0.24311286693286885), (4, 0.4354361205546333), (4, 0.4467752481260171), (4, 0.5308169153994543), (4, 0.5701465671464049), (4, 0.670081739879954), (4, 0.6995486862583797), (4, 0.8186709628604206), (4, 0.9141224600171143)]} ls_loc15 = {'type': 'locset', 'value': [(1, 0.5)]} +ls_loc05 = {'type': 'locset', 'value': [(0, 0.5)]} ls_uniform0 = {'type': 'locset', 'value': [(0, 0.5841758202819731), (1, 0.6362196581003494), (1, 0.7157318936458146), (1, 0.7464198558822775), (2, 0.6340206909931745), (2, 0.6807941995943213), (3, 0.296698184761692), (3, 0.509669134988683), (3, 0.7662305637426007), (4, 0.5701465671464049)]} ls_uniform1 = {'type': 'locset', 'value': [(0, 0.9778060763285382), (1, 0.19973428495790843), (1, 0.8310607916260988), (2, 0.9210229159315735), (2, 0.9244292525837472), (2, 0.9899772550845479), (3, 0.9924233395972087), (4, 0.3641426305909531), (4, 0.4787812247064867), (4, 0.5138656268861914)]} ls_branchmid = {'type': 'locset', 'value': [(0, 0.5), (1, 0.5), (2, 0.5), (3, 0.5), (4, 0.5), (5, 0.5)]} @@ -120,6 +121,9 @@ ls_distint_in = {'type': 'locset', 'value': [(1, 0.5), (2, 0.7), (5, 0.1)]} ls_proxint_in = {'type': 'locset', 'value': [(1, 0.8), (2, 0.3)]} ls_loctest = {'type': 'locset', 'value': [(1, 1.0), (2, 0.0), (5, 0.0)]} ls_restrict = {'type': 'locset', 'value': [(1, 1.0), (3, 1.0), (4, 1.0)]} +ls_proximal_translate = {'type': 'locset', 'value': [(1, 0.35497750169352515), (2, 0.5160959062272675), (2, 0.6817468794150789), (5, 0.0)]} +ls_distal_translate_single = {'type': 'locset', 'value': [(0, 0.915588599565521)]} +ls_distal_translate_multi = {'type': 'locset', 'value': [(1, 0.5795163072671657), (3, 0.24228815992614555), (4, 0.20321157163712014)]} ############# regions (label_morph) diff --git a/doc/scripts/make_images.py b/doc/scripts/make_images.py index 7b1c6e6db20849217854c395c359097517153dbe..9bfa9511a35f0678fadd523650049a2816ab8d6f 100644 --- a/doc/scripts/make_images.py +++ b/doc/scripts/make_images.py @@ -279,12 +279,15 @@ def generate(path=''): label_image(inputs.label_morph, [inputs.reg_dend, inputs.reg_radlt5], path+'/region_label_examples.svg') label_image(inputs.label_morph, [inputs.ls_root], path+'/root_label.svg') label_image(inputs.label_morph, [inputs.ls_term], path+'/term_label.svg') - label_image(inputs.label_morph, [inputs.ls_loc15], path+'/location_label.svg') + label_image(inputs.label_morph, [inputs.ls_loc15], path+'/location_15_label.svg') + label_image(inputs.label_morph, [inputs.ls_loc05], path+'/location_05_label.svg') label_image(inputs.label_morph, [inputs.reg_rad36, inputs.ls_distal], path+'/distal_label.svg') label_image(inputs.label_morph, [inputs.reg_rad36, inputs.ls_proximal], path+'/proximal_label.svg') label_image(inputs.label_morph, [inputs.ls_uniform0, inputs.ls_uniform1], path+'/uniform_label.svg') label_image(inputs.label_morph, [inputs.ls_branchmid], path+'/on_branches_label.svg') label_image(inputs.label_morph, [inputs.ls_term, inputs.reg_tag3, inputs.ls_restrict], path+'/restrict_label.svg') + label_image(inputs.label_morph, [inputs.ls_term, inputs.ls_proximal_translate], path+'/proximal_translate_label.svg') + label_image(inputs.label_morph, [inputs.ls_loc05, inputs.ls_distal_translate_single, inputs.ls_distal_translate_multi], path+'/distal_translate_label.svg') ####################### regions diff --git a/python/cells.cpp b/python/cells.cpp index b6330912c9c1b2a556d876581079b41221b8e8c9..3bb1ea2b037f861cfeb4239d4364768ac1444f64 100644 --- a/python/cells.cpp +++ b/python/cells.cpp @@ -640,11 +640,11 @@ void register_cells(pybind11::module& m) { "The number of unbranched cable sections in the morphology.") // Get locations associated with a locset label. .def("locations", - [](arb::cable_cell& c, const char* label) {return c.concrete_locset(arb::ls::named(label));}, + [](arb::cable_cell& c, const char* label) {return c.concrete_locset(arborio::parse_locset_expression(label).unwrap());}, "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(arb::reg::named(label)).cables();}, + [](arb::cable_cell& c, const char* label) {return c.concrete_region(arborio::parse_region_expression(label).unwrap()).cables();}, "label"_a, "The cable segments of the cell morphology for a region label.") // Stringification .def("__repr__", [](const arb::cable_cell&){return "<arbor.cable_cell>";}) diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py index 2090f16d706ded147ae9cccdcb366a1b739731f5..b6e91ce791d7df815a2caa863320f2f18aaa8f38 100755 --- a/python/example/single_cell_swc.py +++ b/python/example/single_cell_swc.py @@ -66,7 +66,7 @@ decor.discretization(policy) # Combine morphology with region and locset definitions to make a cable cell. cell = arbor.cable_cell(morpho, labels, decor) -print(cell.locations('axon_end')) +print(cell.locations('"axon_end"')) # Make single cell model. m = arbor.single_cell_model(cell) diff --git a/scripts/run_python_examples.sh b/scripts/run_python_examples.sh index 2527a55a01659fb46420f8811328624a0778b9c1..bf4bfb31fd69192de1db9dc2febfe6cf880b144c 100755 --- a/scripts/run_python_examples.sh +++ b/scripts/run_python_examples.sh @@ -5,7 +5,7 @@ set -Eeuo pipefail if [[ "$#" -gt 1 ]]; then echo "usage: run_python_examples.sh <prefix>" - exit 1 + exit 1 fi PREFIX=${1:-} diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp index ccb3d166dbb2b069a10690d68d35e86fb311127e..3bf0076f0cba8a7f6ed0252a31eb4f2bcaebe1ec 100644 --- a/test/unit/test_morph_expr.cpp +++ b/test/unit/test_morph_expr.cpp @@ -72,6 +72,8 @@ TEST(locset, expr_repn) { auto term = ls::terminal(); auto loc = ls::location(2, 0.5); auto bdy = ls::boundary(reg::tagged(1)); + auto prox = ls::proximal_translate(ls::terminal(), 100); + auto dist = ls::distal_translate(ls::root(), 42); EXPECT_EQ(to_string(root), "(root)"); EXPECT_EQ(to_string(term), "(terminal)"); @@ -79,6 +81,8 @@ TEST(locset, expr_repn) { EXPECT_EQ(to_string(sum(root, term, loc)), "(sum (sum (root) (terminal)) (location 2 0.5))"); EXPECT_EQ(to_string(loc), "(location 2 0.5)"); EXPECT_EQ(to_string(bdy), "(boundary (tag 1))"); + EXPECT_EQ(to_string(prox), "(proximal-translate (terminal) 100)"); + EXPECT_EQ(to_string(dist), "(distal-translate (root) 42)"); } TEST(locset, invalid_mlocation) { @@ -173,9 +177,23 @@ TEST(region, thingify_named) { } } -// Embedded evaluation (thingify) tests: +// Compare locsets to within machine epsilon accuracy. +bool compare_locsets(const mlocation_list& lhs, const mlocation_list& rhs) { + if (lhs.size()!=rhs.size()) return false; + for (std::size_t i=0; i<lhs.size(); ++i) { + if (lhs[i].branch != rhs[i].branch) goto not_equal; + if (std::fabs(lhs[i].pos-rhs[i].pos) > std::numeric_limits<double>::epsilon()) goto not_equal; + } + return true; + + not_equal: + std::cerr << "mismatch:\n expected: " << rhs << "\n actual: " << lhs << "\n"; + return false; +}; +// Embedded evaluation (thingify) tests: TEST(locset, thingify) { + using pvec = std::vector<msize_t>; using svec = std::vector<mpoint>; using ll = mlocation_list; @@ -210,6 +228,25 @@ TEST(locset, thingify) { }; auto sm = segments_from_points(points, parents); + { + auto mp = mprovider(morphology(sm)); + + auto allroot = ll{{0,0}, {1,0}}; + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(root, 50), mp), (ll{{0, 0.5}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(root, 100), mp), (ll{{0, 1.0}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(root, 150), mp), (ll{{0, 1.0}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(allroot, 50), mp), (ll{{0, 0.5}, {1,0.5}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(allroot, 100), mp), (ll{{0, 1}, {1,1}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(allroot, 150), mp), (ll{{0, 1}, {2,0.5}, {3,0.25}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(allroot, 150), mp), (ll{{0, 1}, {2,0.5}, {3,0.25}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(ls::location(0, 0.1), 20), mp), (ll{{0, 0.3}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::distal_translate(ls::location(1, 0.9), 20), mp), (ll{{2, 0.1}, {3, 0.05}}))); + + EXPECT_TRUE(compare_locsets(thingify(ls::proximal_translate(ls::terminal(), 20), mp), (ll{{0, 0.8}, {2, 0.8}, {3, 0.9}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::proximal_translate(ls::location(0,0.5), 10), mp), (ll{{0,0.4}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::proximal_translate(ls::location(0,0.5), 60), mp), (ll{{0,0}}))); + EXPECT_TRUE(compare_locsets(thingify(ls::proximal_translate(ls::location(2,0.5), 60), mp), (ll{{1,0.9}}))); + } { auto mp = mprovider(morphology(sm)); diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp index fcae3e4ee2c29aca6872ebde02f6122da1d6e47f..7309e247435f8c9474ba466964fdefcaf8aeb401 100644 --- a/test/unit/test_s_expr.cpp +++ b/test/unit/test_s_expr.cpp @@ -288,6 +288,8 @@ TEST(regloc, round_tripping) { "(cboundary (join (tag 2) (region \"dend\")))", "(segment-boundaries)", "(support (distal (tag 2)))", + "(proximal-translate (terminal) 20)", + "(distal-translate (root) 20)", }; for (auto l: locset_literals) { EXPECT_EQ(l, round_trip_label<arb::locset>(l));