From 57b634d75d24b4ff5f9c9f8288d7931e3727d447 Mon Sep 17 00:00:00 2001 From: MikeG Date: Mon, 25 Apr 2022 15:32:25 +0200 Subject: [PATCH] Mixed subtree processing (#981) * A heterogeneous morphology consists of zero or more homogeneous and at least one heterogeneous neurite trees extending from the soma; * 'heterogeneous neurite trees ' is called a 'mixed subtree' for brevity * this is a breaking change with how NeuroM<=3.x works * this will fix #975 --- doc/source/heterogeneous.rst | 226 ++ doc/source/images/heterogeneous_neurite.png | Bin 0 -> 8611 bytes doc/source/images/heterogeneous_neuron.png | Bin 0 -> 32564 bytes doc/source/index.rst | 1 + neurom/core/morphology.py | 148 +- neurom/features/__init__.py | 74 +- neurom/features/bifurcation.py | 43 +- neurom/features/morphology.py | 253 ++- neurom/features/neurite.py | 350 +-- neurom/features/population.py | 34 +- neurom/features/section.py | 23 +- neurom/utils.py | 8 + tests/core/test_section.py | 2 + tests/data/swc/heterogeneous_morphology.swc | 25 + tests/features/test_get_features.py | 44 +- tests/features/test_section.py | 16 +- tests/test_mixed.py | 2157 +++++++++++++++++++ 17 files changed, 3093 insertions(+), 311 deletions(-) create mode 100644 doc/source/heterogeneous.rst create mode 100644 doc/source/images/heterogeneous_neurite.png create mode 100644 doc/source/images/heterogeneous_neuron.png create mode 100644 tests/data/swc/heterogeneous_morphology.swc create mode 100644 tests/test_mixed.py diff --git a/doc/source/heterogeneous.rst b/doc/source/heterogeneous.rst new file mode 100644 index 00000000..85e014a2 --- /dev/null +++ b/doc/source/heterogeneous.rst @@ -0,0 +1,226 @@ +.. Copyright (c) 2022, Ecole Polytechnique Federale de Lausanne, Blue Brain Project + All rights reserved. + + This file is part of NeuroM + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + 3. Neither the name of the copyright holder nor the names of + its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY + DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +.. _heterogeneous: + +Heterogeneous Morphologies +************************** + +.. image:: images/heterogeneous_neuron.png + +Definition +---------- + +A heterogeneous morphology consists of zero or more homogeneous and at least one heterogeneous neurite trees extending from the soma. +A heterogeneous neurite tree consists of multiple sub-neurites with different types (ie: basal and axon). + +A typical example of a heterogeneous neurite is the axon-carrying dendrite, in which the axon sprouts from the basal dendrite. + + +Identification +-------------- + +Heterogeneous neurites can be identified using the ``Neurite::is_heterogeneous`` method: + +.. code:: python + + from neurom import load_morphology + from neurom.core.morphology import iter_neurites + + m = load_morphology('tests/data/swc/heterogeneous_morphology.swc') + + print([neurite.is_heterogeneous() for neurite in m]) + +which would return ``[False, True, False]``, meaning the 2nd neurite extending from the soma contains multiple neurite types. + + +Sub-neurite views of heterogeneous neurites +-------------------------------------------- + +Default mode +~~~~~~~~~~~~ + +NeuroM does not take into account heterogeneous sub-neurites by default. +A heterogeneous neurite is treated as a homogeneous one, the type of which is determined by the first section of the tree. +For example: + +.. code-block:: python + + from neurom import load_morphology + from neurom.core.morphology import iter_neurites + + m = load_morphology('tests/data/swc/heterogeneous_morphology.swc') + + basal, axon_carrying_dendrite, apical = list(iter_neurites(m)) + + print(basal.type, axon_carrying_dendrite.type, apical.type) + +Prints:: + + NeuriteType.basal_dendrite NeuriteType.basal_dendrite NeuriteType.apical_dendrite + +I.E. the axon-carrying dendrite would be treated as a basal dendrite. +For feature extraction and checks, the axon-carrying dendrite is treated as a basal dendrite. +Features, for which an axon neurite type is passed, do not have access to the axonal part of the neurite. +For instance, the number of basal and axon neurites will be two and zero respectively. +A features such as ``total_volume`` would include the entire axon-carrying dendrite, without separating between basal and axon types. + +Sub-neurite mode +~~~~~~~~~~~~~~~~ + +NeuroM provides an immutable approach (without modifying the morphology) to access the homogeneous sub-neurites of a neurite. +Using ``iter_neurites`` with the flag ``use_subtrees`` returns a neurite view for each homogeneous sub-neurite. + +.. code-block:: python + + basal1, basal2, axon, apical = list(iter_neurites(m, use_subtrees=True)) + + print(basal1.type, basal2.type, axon.type, apical.type) + +In the example above, two views of the axon-carrying dendrite have been created: the basal and axon dendrite views. + +.. image:: images/heterogeneous_neurite.png + +Given that the morphology is not modified, the sub-neurites specify as their ``root_node`` the section of the homogeneous sub-neurite. +They are just references to where the sub-neurites start. + +.. note:: + Creating neurite instances for the homogeneous sub-neurites breaks the assumption of root nodes not having a parent. + + +.. warning:: + Be careful while using sub-neurites. + Because they just point to the start sections of the sub-neurite, they may include other sub-neurites as well. + In the figure example above, the basal sub-neurite includes the entire tree, including the axon sub-neurite. + An additional filtering of the sections is needed to leave out the axonal part. + However, for the axon sub-neurite this filtering is not needed because it is downstream homogeneous. + + +Extract features from heterogeneous morphologies +------------------------------------------------ + +Neurite +~~~~~~~ + +Neurite features have been extended to include a ``section_type`` argument, which can be used to apply a feature on a heterogeneous neurite. + +.. code-block:: python + + from neurom import NeuriteType + from neurom import load_morphology + from neurom.features.neurite import number_of_sections + + m = load_morphology('tests/data/swc/heterogeneous_morphology.swc') + + axon_carrying_dendrite = m.neurites[1] + + total_sections = number_of_sections(axon_carrying_dendrite) + basal_sections = number_of_sections(axon_carrying_dendrite, section_type=NeuriteType.basal_dendrite) + axon_sections = number_of_sections(axon_carrying_dendrite, section_type=NeuriteType.axon) + + print(total_sections, basal_sections, axon_sections) + +Not specifying a ``section_type`` is equivalent to passing ``NeuriteType.all`` and it will use all sections as done historically. + +Morphology +~~~~~~~~~~ + +Morphology features have been extended to include the ``use_subtrees`` flag, which allows to use the sub-neurites. + +.. code-block:: python + + from neurom import NeuriteType + from neurom import load_morphology + from neurom.features.morphology import number_of_neurites + + m = load_morphology('tests/data/swc/heterogeneous_morphology.swc') + + total_neurites_wout_subneurites = number_of_neurites(m) + total_neurites_with_subneurites = number_of_neurites(m, use_subtrees=True) + + print("A:", total_neurites_wout_subneurites, total_neurites_with_subneurites) + + number_of_axon_neurites_wout = number_of_neurites(m, neurite_type=NeuriteType.axon) + number_of_axon_neurites_with = number_of_neurites(m, neurite_type=NeuriteType.axon, use_subtrees=True) + + print("B:", number_of_axon_neurites_wout, number_of_axon_neurites_with) + + number_of_basal_neurites_wout = number_of_neurites(m, neurite_type=NeuriteType.basal_dendrite) + number_of_basal_neurites_with = number_of_neurites(m, neurite_type=NeuriteType.basal_dendrite, use_subtrees=True) + + print("C:", number_of_basal_neurites_wout, number_of_basal_neurites_with) + +Prints:: + + A: 3 4 + B: 0 1 + C: 2 2 + +In the example above, the total number of neurites increases from 3 to 4 when the subtrees are enabled (see ``A`` in the print out.) +This is because the axonal and basal parts of the axon-carrying dendrite are counted separately in the second case. + +Specifying a ``neurite_type``, allows to count sub-neurites. +Therefore, the number of axons without subtrees is 0, whereas it is 1 when subtrees are enabled (see ``B`` in the print out.) +However, for basal dendrites the number does not change (2) because the axon-carrying dendrite is perceived as basal dendrite in the default case (see ``C``.) + +features.get +~~~~~~~~~~~~ + +``features.get`` can be used with respect to what has been mentioned above for neurite and morphology features. + +.. code-block:: python + + from neurom import features + from neurom import load_morphology + + m = load_morphology('tests/data/swc/heterogeneous_morphology.swc') + + features.get("number_of_neurites", m, use_subtrees=True) + features.get("number_of_sections", m, section_type=NeuriteType.axon) + +Conventions & Incompatibilities +------------------------------- + +Heterogeneous Forks +~~~~~~~~~~~~~~~~~~~ + +A heterogeneous bifurcation/fork, i.e. a section with children of different types, is ignored when features on bifurcations are calculated. +It is not meaningful to calculate features, such as bifurcation angles, on transitional forks where the downstream subtrees have different types. + +Incompatible features with subtrees +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The following features are not compatible with subtrees: + +* trunk_origin_azimuths +* trunk_origin_elevations +* trunk_angles + +Because they require the neurites to be rooted at the soma. +This is not true for sub-neurites. +Therefore, passing a ``use_subtrees`` flag will result in an error. diff --git a/doc/source/images/heterogeneous_neurite.png b/doc/source/images/heterogeneous_neurite.png new file mode 100644 index 0000000000000000000000000000000000000000..ed4d08114e8d4bfd448bc936e05c45b08c02784c GIT binary patch literal 8611 zcmXYX2UHWy_jYJOF!Wvm8tJ`BFDjuo>74|mNC#;FLNgYSCSCfAG(+zoC4^Uj7o;e? zC{^i5A}#bU{{G*dGrKc)@7(+B-80Y3oI44Rjdf_KIjI2v0F9n5)C>S1F(=knZvlz- z5OpgI@j)K+NYCOH@j~8mjVHFL0(5PI0027YzlG$3Mx`fF$Puh%6KwA99vtTU%ncA0 z7AE28=NT{h!@;d}5P$F%K*fLO#bd@t43DU! zUoa{^q|(iyWgy{-hhjw?@`P?d985DspT|!~=u1*-=Lu=OAtxg#=eYAg_ZB}D|BIH8 z8yoL;QLmKvH%CTPnm13wi^5NoWSdW&HeH*>7Gd$ zG(KKV6vYqP7}SoiuT0T_mI-r!M)&kP%{D#agN>(=o?7_>f%+s%hMWoiGw_}v@X~__ z*_sHM#g}oAO0V>-gQ-kbw##F{tHE~t3JBzfvlQRtnVe#XYC5k-u70(!6@FRbsWHYP zjR+twdHutfAXzYNSOhy#mZtmrZ$#$hQ|S1TA3@7(Iha99DPGrc|G(QO8Pu`00XX$%YSr0IJht-8!k)-Hxr;CyDn-8fSj-S!u_kU1ndVYe_ zH;`Z!Zg{ZnxjDSXjG4jtMAXm4sBmdMFN;l&;O#sz8jNp%iF;kl z2_Dmj0&2)AeunqhWZFT?lsNYC^#TlzR+L3Kx!nWo1(eu$4?zQ51>RDff~Ulg>H>R{ zUz`s56)b*%$Dq1sHgfpCT7$+E2;F^Vt?`b~G8)uAXhtt!xgS)a{YS-xVQkkUprFsd zjp|4oRsH@5?bQ`OmD*PiVYq_y^|U8NO17j#P9+}jx&#EuQ86=Z=>V;{_TYKaLJ%~2 zQgWGG^U}04ig%z##ldnJ-q6SI=36pZN)vn>(Ju;_sM|6s=A!aCRbIYG=s4im?_L^T z=@yD&&Db&r!lbv%$8EHr$}FE=jI2CYUJ5$?u$13d!Q!MbwQ_987_|gXs`>+&XaJ49 zd2tK{ADoSLZ0^pFz1MAu8RPrL0UXQ|Xo%M>l~q(^oOM9)?*@;ntQYl;sh+Z+Y;%p# zUub|+_jnl6cx2}k%^AOU9Lh!@dwC+dhJQoqj1UD#KW|?6MbW;E6 ziuOy=9d<9Ap+K6GTb8vqBZDW>lv57Hv>Q-wzvdLi+06kPmshNtIb~%)3!}OqXh)0n z3o$p`!|&CK6i05T@u@!;2qzG7ETfcy@5Hh0OqEIDI62pGxT(-6=NJme7sSex= zw)3C^=boL*!v{$X>+xz@a&5*t8(E4o%DEVd?gYVuD{86wfm0?1G%pG)yRPCbKO-kw zy?wj2^R+I2)rLV&H}@dveJ_RkHrzJCn5bsfQ+m(TzqY?%b@0cP-odFMLiZhws=lc~ zMa*QfI)w=GU;4=v%Cy_(Wl|=oaoWkh@0t+@>6Ce3NyeyG@YVhr;X!hpxeixZ!UF_wYWwx~-ew`x=9Ef|0?DEKxus0lLx~;ANEBl>k$p6xQy#p57 z+{H5|pAu93Z!3{cG(k_i!nij491gZSW@f7yigm+2xQc{WUd%~0?Uk-c#_*%Yz}1USu)$l-j{IPUV7Cbavj zzuu@3zMSHeAwMcyMTDK0&1VGC>so#J~8Cy;5XM3 z2&EApn?BI9jHy*L-nNB8AVwmx)-i1%rK6MZmm}|{S6kyvp=GNaN7C=U{ct-2g(e4S zY&myT>!yCllhk8@(*|fF7YnPd|5H?V!kt*B@evt^fjRQHUBM=jVf=?`%FPECRS0BM z&}!wz#A3F$H+~btunD%PCJ-a?s&6}u8Md{fZwB;sJBhMu9=Q7^UQQj=D%m{4@>TxWaaIH=Rl+8r6Ek}J}wX5vgFp1j1R+`MV$M=( zj6bCdyTgWbJKb>~Bl_nOlJTu}%l}|@$sMR2k|bMtG^CP?Bz*}Wu;9RKEM2U)%(xMG z+wMTmvou4?(VSPNxDdVi6n~ZKoz%}vGpd}Gc{i3Tj4RUcgONiRDJJXe`W=E6E)*b+ zyAwd$V47A(O!2+Oq^=6E2^5m5CE1M*KiS9edXJYJFtLz?d*F*Vz+ull4 z<>yEqvOd2qV)*-E4;~$5$ z1|1j(-w!z=8z39(KZQ|TQ`Mky9AIQ^3A52u4zC5<)#4QJ_nj}uZb(A^9lVmTT`j6# z-TkUgv&TK|W>CSlme6szK!G&GQDNU;>0(obvJC4<#u88c=nS=0#@+;*REQt7gS$6R z5aT909QKs7Sc8}q{*89rXS){RB6b0?{Qe-yZ}QtGg9e*Dcvt1ShYPw3#qR=yrMM83 zl4zKYr`LWD7L$8g^jy$A))vc|JL~tpCDqEh>P!>QNIEXl6uUz zLfUN+j@US^T>q$mPlq5{)w(98?Xi#a!5aP^P8V-e7M;iYRQ9#k_&Z%74pjc6LJ?re zWAOz4Lz&ipy%X20RO37rpXOj=*s|(=@cNs#A)?it8L&6r>ykiSA0_J)DWK)Uf@C8^%?DAadhP@EP}Xr5 zYnJnRi0wwUCEMfM-NPC@gagGGxpH?hoU2@8;f2SdVV^;Rd8{d-^@0VdRkw~Ox2su* zVOmyp2Q4}Lk#njdhe@t-BIhj2C%hW-V8N8vgS9rJ8q^4xp#2zAugD(8`jdm#3gX?G z+o{6|goD+`EAZ~j-e?t#_W0m|71><$`~1wF&=A!a&&bW$c^(Tt828Ke8!!N=e0ZrH zNQ*eQGDQsiKFmD@k%x{Wdr9}{!WN25P9pp@g&S!;Fqq3>}63P9gwtbi?WLjo$uJ zRLnijyR`HX0$~KDB$X*mx?Gm`y){f|?at(fKz@PF`>D^}as7zjTR?9~7LuuaNtznK z!5WxlMKcu;nzXJ~tf~oxipapCmR+o-=uz$K7d+4~W%3k9lr7)QY!=zNi#|LYC+iY*FhccE-nB(qBznF!wxQ!Dnc{@{C4h?ZPtyeO8k48fqvPK^K}rNS7!Q@MJSR6 z`+6_3EPf6BRR8Qr+%kh`$LmkqOFTa|$!?~T56j4aT{@6L=d&M_TX*D`Mc5+N^P~1j zM?QycmJ(0q!&ZQBMKPTg8@GoW^DappyyC7CgN*LQ&^Gtn;Pk3X1Cy25i zpZ(n5s6ZhKH@@VFt@+%jOUJUgVe)wp=a}H%tgKoDqa)Bb(R^f>6AAkLAL94J6=nE7 z;`;Fl9X5)vHb4~j83cWb6V6ocHd#;b@#yAMXIkFuMEC{XeWep;y}lcW%_Ht&k*OcvwN6Fw?Cw-yFaC%J6sbn!^z_Q6u2Ds5sxtg$ z7^lNBA-{POV8>%lyH#j}x%Sy!9!X%2{;U%fLhaB_{#6`S9jk}8d41ldq0loo^O0~+ zJScK~t6M(w+0KKq2ocL-_4iX_NwAx?8B6%)r!_AOH0}+&*vVK8*FizqlT6@=mP|*w2ZhiG96xL4y|{ zSwk>!#{OW61bhKZ;si3f&y;Yznd-9XWc{wntp2B5@z&0CDVsMm$lH>_dz$FXk8)1s znA;{8u65qGXnav=Y#5Gw5p1$Xb*Xr%c?#~fyCTJD>ho{j#*LM5AcP~umfWeA(4iKe zPj5L{`j00O!pBN74JENEZA%H>F^OjR+h>G3eCCB?hk9?e?~@G z#~{u+-o$8@T!o&vhupaipRaCK&=1s(r?T2DJGr z`=vSSzORGuW@|#-!I7*^od~|V)4RA+>Oe7U^{9zkQ9k`?G<^rqf7X%9#I=76jb@n4 zW*FDYeDz2zL9jUsABdb6QNgNqH`m%Y7RTN^{3tn=pZ)QqDR941_8^LeH$#QtU@zis z-F>BK%6qZS`0rqZN?q{3|L_J-H7EvL(&_u+!@pAA#5@xyw@Z%KJ?22(i@+CuS#`KH z$7}o-&Ruv1lKgpvH_gz*5Nwde^7_$!=3lN;LOG}DZq5t7uDaj#CDz`I@-9RiyYgJO zqP-%eeE3cEBN{z$eEiRXK1S)qp%IVbUqexZpNtG?3o`6X3Bgt2CecpZv8mDDsvbkk zB7vt+nbAqsL4y^|4NRIceOijx?Y+0kx3-l~tsqvYjH5<~s%0>>`fpp!lCqN8#Bm4C zx61U}94KBbNKgA=m_jlM;qiw;b4thri}m|`p1iz1Gob&khV^ZQpZ*rlf7%4vVznHF z^*LYk&FPkn+b?SM=Fbm270Q$=Szb$gY>**sGqN(TTk*q@xsu_ECQ6{3sy43Oh{g~c zFP1mrb$;+p(jtTMVqLSOvoob$SMwt|jP!9{eBYH47p)~Xb%+IV=IHd6a#98pW3$bP z8ik~dwN6|+kN$R*W7dDRl0#TBm>>t$CA&t1?ATPf0+u+cN*+%fS@mm-cnhRxDmzOl zTD_7vETsOkVy53E9xvB(I~GP<%PYGFnDR|%gkArxkm9C!nO|031lMwwHJ|nQdeDjGFE1(JuD1R5 zBh;w*pgDyBU_7IqiW1^`Nr(B-Bu9_(es2*GLneey}?wh zx&3SRr;2&v?~hUsxdUkpl|SS7GX_xfxVU;BgzVr`oUTASTr{A@qd`!qE}?8 zE7h86IYPj+&U2Cg-<>;lxqWYdryO&yj~b&#-jI5w@&;jx02gX=!f+ViBDgTRrYrVv zH(1h>FM+82fU#!_ZYD7Z*wCP_atU}vjwnz~~@dt9lV zI&;bQ$W_x6vmK<_N%mfjEAA%dJfI}wS}Px+!&PM8zODYFX@;Se*5|elc+mMleUpMw ze`Vt}M+YUftT%6tZgjz$lXWyF>KJrAr&qe6nT@ooKPJ0Som2XTvE$irXi&joXxm&C zY|7_b9}O{ge?Ed>oEJ412PazI=S8yyfEp00SIjs?^^HG{-R)PNE32e7Zxk0j$YMvn`gOXm74qEET9o}NyAVr1_wGjHz;gx3K%@Ff$kqmKhQ}?^Mzb}-ui`*x| zHX(#(ZXZsG3YnVcpGeOh7!>u(RetO4($LhGB^@*bQkNpM{qmwN7YDY9&r!hZbNFX$ z;ns2cR8~vZA`?v*vhQii@$U5NnkM((r45Ola@1rfw=Fa|!Buofqw2*Pje#m;KBH5n zZ^~vU5j1lWA$O98;bndNz~(QnlLI_8l5-^RcCG1Nf~YOv-M*{{%9n zpY1RUF7@AvD=H=@B+TT#;twRkHnDjcYYdFU)U=$3%F!VOLS8_QiKbZ&rhRRJ23*`; z0bw;879+H`X?Wxa>1Ppzf2i6#gB^Nc*)}&6q3N$|^+Y*ttG-o^X9l&&cl^83(`2?~ z0Mz&|CjaUE`DbvQI5n$m3ypaBSyRh*f&3F+vjI+|>N(19BcyEpU5RDaExv+jd|44e zeXE`2vGbKzt-AO%X?H_*qSnxdzL#m&=+<{$nF1;I2gI0wq2m1w6_cg=vKJbcl0aC^ zl0&#~qdNQv_IFB3l~3igBi^of%2eNE9^TMM0Cg(C#e3e&tq}^Yh~#JH3%n7ZelDIfcI{7nma(_konOT4G&Uug9 z&2T+ufeYr_KYY#<4ddD>JCtrWRK-9Q9-m)Ce`ze(88Q{nUJjgT*o$L6aFA_nJ>4GD z+{T76JNHAAZA+k@_Z0!@a;^+gQyWVX!=hQxy!N!b#pjRbV@5v0V@Oq)4K=^*3}+9r zeO3bdIJ^u!SX*~Qj^JZTN&;BGv`Lji`jx;Q&;Yxk7}rd2*BN8*1kQ)7ffr>YeqnNP zcP<%zQQAsUH4>rav8sTTpuA9+`!_UwlereX(m3y=F=;bJ0@{Kr$M*p)ZmZecM@$gwT(ei zy{s6?6>eQ@zqx_Ba8c_*un z`^y(NQgxuHbKajr} zy6hLS`+7$b=k)W>i$?lAVdNT?z2WZ9(330{Y+GOCgy@yT8?e@L0<^lEX}2c!fFT4v zeCK^ww8mk7PGWINI4}pzC&J&xFxLogUt6~lwRH2)O4t)jgXbfw8frQCzp{y3)|foM zqW3YFGM%D7A-mRkKY>0|{Dx;)fW9t^*U5W%Ou7g<2*H4T>?nO-n&YX+^2<{?y=$y3 zMhCxz)@(j%pnN}S{xS(@0T2MdYBr5%Lgf){wV7~Uk_(wR z8Fc9tAJ?sus0Qm{mm=nJ^ZOz9<0Jv&=C6H~MyvzQ?dJtAZ4XJUasO&*#w8$rCF2v}f+zqc7vI*7{Ed~$&~dG9 z9VHjCsz@XMXRoT1?^tI(KK=<)?#J;NpU9OR<-mS&o#Wa49t>GI7MofZJJG?Ncbde>+SHbu zgL&26WFu?m35JqckIH3Gtuzv;pRF!`sya91ScX4jL0OJkzwR7c@cZI`yqsK70ajZb zVkSgAdPmD?dtOU* zo?jPIkif57WA?QnxxW!51s`@890;d-sjgGgj*4KHe+v|_xdWgISx%CFs* zC?3~Z)1OgdHX`+4@3uLiDs~Hq^K^|Bld%`q|L0Ja7oF4R7r!YnEzT7!x#Ne+@(cbRZZ20s`V*hb28BikUadB#H~FI3O5I4yt4FK@ z7Z#EJXFi?>2f#>+E~AzOP_{H7`;0K{ve&Ow-Wt&%23cf))uIW3dB9McoeyPkxpbDw zb&s;n9-`FD2irfqZ|U@(=Y2`BQA$2#%b$p4B~$6Hne^nI%R}GG0J=29^=SbQJkrNh z4+Bf{D^uz@M?74wlnxI^`> zZZD&0`gi&p6IJFpp%m|ypkW=XCN>$HO!3YsRgo*}Y^+v;1oLCeo^`FaR3TOve_V}q zUEyEJB{i(o{4nvgMj?O`r`8R=wuJKvsTRqJdBQBFPh|*Nnev;Vk6rBGxbASBi7=W_ znOa^u(~1|$qsY~){nDRRXH!kLeB|bi2$BSKoA#LxZ#iXQEIMDafsX-F=33nE@jPEa->lsW~#+}gu+yk_C zY|7?sVaQ+TwdJki5|$&b#D_5`b@cjw_YRe=-PMCcW_|a~=&S%2o&qIKFSA!Otl$lP z6;&ElhqBPle^XL-FmaA`j8!e)>nE##yAa=W98ak>+KS@(nfLik&c@d4<~ zuhC}d3P5l71tUsg8vD-(@cC#;>S7wYlmRVr81>JFjtn?gnP*5vMI(KiSxS6O;>-mz zU&OJOPUH`jUEZE^V5rZ|)4wU#MVnar#3Ehr+{ZLSLC<}2hBUwCIC#*#lHd8toD69K z?!Aew=MQ#tiCL9*$sGKSVYf(~NGaVRDBU343QGyXg3^mfcXuoZh#)21NbZu-%~I0cy(l4F()BI= zzxO>a=WurUotfXZX`0{{SgCD1EP006xR06>#?1VlZ-E0m8% zePFuCDS;oMetaK&{EYe=#}TCK0ss)2Jp7^+Nf%n79#Xo#e($Q~VD9Q+;%o-+@bKWa zvbS+DHE}fKc5t>x+n1mM02lyDuVlfV>AOoFzWP&hh@+JTI;XRWj89J)Oh_K(+FH(q z@Ea*JxGyH6+p*?G%H@ork;rWoWjM3eZREPg z#oo2l^JiSyN6D2D#j%!m9yDx;fB&R!rkiRgv%QuK_mkYMs@#~F<}cuzCt}M7A02<^ z4%Ad$4<@?#0+~aWMpSYtvuW=N%4jRIAqw-J)!06Kg*~zggvV1PpFg@CBD55&gE7Six5J1>YlAc)Qc#ucgT2@%fn&K zd1!D=_MH!1BLJV4?l-ir`f&Nk4#?b6I`+-`66<+vgU0~cq`{EIm73T@IoKOZ)cpFz z=w0X!>FXMqJV$~U?^;e5njx-!2>#s5uByJ1n(8M!}^09`+MQN$~j>22lc+T<~KDu(ybGA{I&qz zTa}&%O(eN!1$DXCzfh0x;u4-w{qcSRbTb0Jy|- zf*g$dVVS=e(EgsD*7PfYRt3Jy?)Y%G44M$G^4lZYEo-lAAMpyLVsiSLu!&#sWOo&% z%Jd!QqVN1Kt4mP3PlYO*J^xq%e@Ynr8&oLpN3yi=I5=&{Mo*Phrqq)uZAN(pG_ekB z-O<@h`nC6NAo-n(E>JrjCBBzB3IkvxSpJm};Q=N@Qp1|BuIIfkxVFqfOQ0To^r{0W3!vkV|E>~cnh zlP+6>0ox*3p>m7gH?*)qEAfgt#;o>Y$s3+G1y_3phQ_2V*+4gwynGEf?my|rKPuc) ztHz{>lHAOBRJy);QSQWp3q&LSIJCih?@Bez(LM+W1x!QQ_)_@F*M7|_i^>MVhpyQj z9-!W@UjtsUUKO@?3wx=zL*Hp`d$U2-&r-SS!i(D8P2&=M9Z&%9POtONd+xGg@aO4l zcyr+c63bg;r_LR#-cSz44l(%!Z&D&%8W$nJ8 z?2+5+e+o+rJsLnVGyE|Lc#Cl~A>d>mF*nC`Rjf2T&?VZWInslHPs45V!ks^XqtZ_; ze?0QJ(A1XVTe4TMvVVy7`Br4K)&3Df1r2X8CTFg|-T9n8y?=I!Qp7k34?`L|etv}kG6HhGnsWPamQ%%b_`^5fKha}`K9+EG5V=BhV$K1l--8UGm)#e4Bc zV=s_%I{W3gZ&^Fdmd?FY$o4{&8&`hghkRf_v;NRk5vlTN$iM!~C)3y8Vz)F@@Xv-R zj#_-Yo0X=R{lvnz6XsZdXK*f)?T##M%Viy(R65b-ds|A_li(mTSaH9DJswa^TEZ52 zS>BzA{~l5Vy*6cb-<=z~OSE6@w2il;nm#$QZ&$N=eD+r?l&9}!P5gGwp^tppg|^Ye zgT6-#c)dg0jR26G$2Wyk?{6CRqpFgJ(D2MdUJ%M@r7pi8_;dYlV3U#2V~z(K zKsWl(_t2W+)`wR*K}gH$kj|FD-@xKpzodOg(Vp=)+O{p(fj#fL#;MDEN2P-DoXHon zol86~D|jC;i2s1WTRAT6BiVl4fYtW(4f5WY&#Q<(9I-!;{&d9}bOd6QLCcnJW| zw>@(O{RDD2WN$Gdb${DbMQOa=kqfx$|~DvXfHn~eu^g} z#HuuZo1Ct14=RmTx}U=%+gAdSrvDbpz^}AGA6SZpx1m@QIjzT@pdLh^snRz-YOQP* z_pn|3?f6a8tYK0;8jX^UOfh*&7FjdnGd6;8B3L4QrBW$*$HIP+g_@LPm_;s;k<=yN z|CXA>n&Q%gYs%u1`MhGbFlaCI&ZLZX|3u8im2?i7yHmha0R|IbT z*T+?&nSq`(f^A{)M3vboBt(93+t`M}PcaL(5WGRVw$F9A^00r^9}FV9oImTqq`fA- zHB-1Bp24=x3P&)mfCg`V<*JOweJFI}0ilHdT8hyvhh|;_1+wdUW+nu?J~&|a2fpMr zo|yX%ma62XxG2%TRgX7&3m5)}3)+aeC#qC8CM}Z7oi{kMM#E<1CHdQPcx37@RWZON zBUWkv{PE`pm44Uy(k4+Y)-2rfGDd_{WUf@o3U+7q^Kk)6SQ5Qyl$(S1uIifT6#Kg{UN zbt>Ar({(yQ8L^7%T)69qJ9}*#b!}O;uZ;j@Hd6d5l3m9L%BgQH97;Dtv+C2k^wzLe zj))K~c48?)Jm-Pv*SaxkCT6?SOM&4CF^noE?zn!&e=_w!(zt9*D-xFW#M26L{Vt+P z=4Sk+jpO!?*%HtcLlf!4OGsKE&$&gP6k)tTOKS?FdxbIi{X6PJfwSt`R;3vf4ggN= z$)c^g2W8xiHh*V>^ zI6%$=L79xuQ@zhC<`B-n5QiDQvzl}gb(mDeu}4dxz8~+e%6KNVQ=g^3f*0VQWdzd8pbzvB zE&Zkbwrt`YJ+HxVo}i7}m&kBiquve*W(ER#eUBg4R7jybg}+>ovL-)>Oqi_}^}mHv z27UwJ%3!T0AVL8(na?}9`otZzQuJ=hC;5Dp1fOeLfAei$k~W3BexbhWJrh8X9EsV{ zgiK``uLhI@cs|)XP~oH|Z2yDmGqYeCu*5t&4*rMz*c&EZIs~Mw$ZRmc#aJr;<<|r& z&z`2KNTCD$S@JDI<1kl8paYZuE4HGM|2au}aB+P&82~E77l^y6X(_EM>hh(TfE z3{him+D!5WJ_b|OCVbEPqL(s!`Dnit$Fr;syUFPMTcAKJ82WePS{Mq16hP0fHQwZl z<)gvbz0WHlqk4nL{76H0vVxjoze3XixUQ%ze`jjv!p@I%E>4(8O8*;|mI9NO3|%nJ zbDJHC7I@VyU1uWq$JNgXwzbVo$L~w5AaADt$}?);z0KZTqYs+DUN6Qo6vO}R;dUz$ zp+L}8c1M&}zGd6CHCq|WNj0|Q-|}|k_$`K|)n1>=76)O|vy0`LWRvrw;Xt8e?~i|| zZIFMJz&yR~qF@u3_C?`97p`54ycF^NV?23aBlblS^v%jn28}rB1sx3S7(-UWy9w&S z|5hizar~`jl_nO8!p%$M8Wp@bD)mB=Xb0b10bjciYk^uc%k^aK2djYz;*jf~R_Vf2sS)7q+a;HuV(Cf0U0CT<-^w6fI zM9K?4{Hb@YS6R!7Rl0)ng&+9>LsanGsd;5}1N2o+U(YqxZcB6e5+uk zUR>3&@#dE-3K>J3$5KrkGF%P>rXXCIP+%4P?ba5dc{{SlNo{pRhtb57>_F9wfiNIE zYZi*E7n*f=bU+0`pQWdm#k%PuYsIX-NWRL!yDSdL0iFZ^#jz-KgRF#g%$dTiAOZGB zdwQRk_Q^olaYgjti&psxV~lXEb^{BfY;Zt!dpcl)p@KYa`w$SQ_)8tgEx^0fbDiET z-v42q^K4aIw_eS1X2@4$5K_*UJm9x03$d?u^d=35wZDlrQCUDrRcIj$HcZbj?gtQO zbbVn5X6>PWKi5QF6hm=2&j^7BfCy3i-AuZTBk#np&6<^QYlOEit(Qb3belmHHMQJ9R&eL#x z>N+9Y-vPpv##+YUeny%;ET1_*2IL+^f?vy`9Je<9v{x=nN7bd_Ux_;Vt;C6E`*YVK zzQORMy0WYoq-GHwA2{!HxlE1lX)7z^;7Va8^8hMPJhd$oZ&~7YaB#*4=AZ9(@xO9r zYIG)_AAH>U%xsqgIG@q*Q>a+U;Pe*?KTnLlWYoO=(ryH7=1+lt3(`aC3QE6&-kV%+ z)_Z!;5A!-Y;ZLm)D9tacXkOVNo7&wmw566;-M~O&e)Vn;jr|Wd#8<$OU;viL zsHXEI348Dx(vczaK8dWgVy^gUw7SWO=Bn;{s#)p+-m7t>TbkN0V)|p`&Negz5KZ!)iB zn8cE;d%tU4Y0qZW*w+Gmn~ETH~x`qe3@)q+csxh zcwm9#g_y94xb-aPuR9$aR?F`b25EWtKG|JbxNAOB1JFBgh#K$QA!Im)K2oRmCZFR; z;Ga1Lrf^c>**fI#rGG~{>(N^p|0J4`p9v_4?IsRNQy*YOx_zW-uc(u6x6&AX7rSwX zEa?FZs95>%_`Xrup~6)otD+7JlN@Z3+RGgI2&co6`jPy1rj%v5@^FsgXho1{WdiGq z=!bV`?{@(W!kd$CC3!q%-rc@B6xEcoz8Pbsk+C+f4sfpy4he>mLsz~%eq)RlC0x-4 zZ4wmGzLY_^nQfZ#>pbH9vO>)M^Zx4k_`Y%ezHFrQ*Cjn47KGIFZJ#w8$?Y^U=TqCM z|0+E>WNX*ybutoqFL4t+QEq0oICKjDuG&TYMK_vk_s3Ohatmt!X8A5q5qyDv&%==PZ8t&ph!b)=}Vc);7yDTrixr~Hm$;bs}^%fv)_SLTn3Y% z@37d541~52m}htI(z~agk7`dIVw52_HCT534lWEArs*t=qTpcYdUM#DKw*tLW6nmgg0atFsu7aZFUDs;x=|=vTJ_cRqoVArmxAj7crPUoF9ni~O@yua1xL`$0i}xH7!{mH zexhK*y|lZP8Etp&OJ%i8G)wTP%^nxZ)p0##I?!HF$N+F;H+?t5zX4?f8$UX8hV+>K z`Yk>$92{gYfu`FUwS6Cye!679ZTUubmU6bQba|^bK`$A&sK_qj2b}Fh5DK9-jVcCM zJ6>Yf7%+cW*mR@*K~OzDbFDco8?Ihej``-}pbvx|bGVGf{}CwKEBeb$pmp%2X+=5I zhKq)R$zM%Ur54U|my&~*5rSt^(@w;GfVB)L!ym@sPoJ_S3Pa^-CSQ=8DnP>3ou_Sy z+LLHbXO9(9txW~@&@{S>6rpwaH>?ULujUG^K6hskcmfN;hA;6dJTQigF)1_qnVe{? zgl70s-{3}@A5zqUy{9Jl18dcVR~!BY;TM>2j_kXEpL#AW5u(?IBd=IO6A!92Lkj~r zVh^rrK%0FS_=7`!NAcUVs082J)WPI;h)gURvCna@A}^YdZ9wMs@-6$on`PEcRjlPB z2ylw=vA)Q{Npv5f55g|8`5Ng^0$!e=0Z;Kk;e zpeLi3;#YW3mhPDcz|H;+eeF8HcKOt;6YRBOv4`OnxxkzBm(O-)T;_vAoo-aJB#r7A@8iXCG%05z4OL`NSj0JTDz4m_);9-%Ts2M%E<#+_dB&-`s;Ows7>X!C=A8m5V*WTTPLPt2q%*pd79 z7MM!Mf5?_x+AaJdfTC;d^%qUI`P8ENf~$k=sbT#lsybmgQNPnJoKB2G_*G3 zz6i1ZGHOQV0*zj>LAl6H;zn$wv>2s%p+*bei~R)pIw+7QIcoH1*p`m)#~IVv>>A)X>KUWfmcc{Xu@DRJ4!y(V}TgfE!hhajc2#i3r3L*%4P~rC^kxfq#~ce*f3KR z%5NVMPAK}q9J76u0XNs#b!atb{c8cc4qW3O4WVrI6FX9t+iak*$|)+o%+h@7EneMBYtR0h^E$FApnRYuo}CvZ z_Ty#m{!BWVoj7))q_a~4pv$vawG{U1uw5K77grCTn)Xf%OGBc?X{}l@1z+9B3we}T z@WG>=iyk0~(?+C_IqYQaIv`cZ?s(oiFt_S5gc!;LXW)z4At8u3Ep}eSPVwHBzYA2X`!l^R^X#FGYHF540}e4-J&|& zm{c$#U zW4uO-vW8 z^R}T6D#^$PZO_Rp_PkMBabT4nepuwWE}KyPHKvGr=?K5FlnsyV@B`9i4~Dy3po@>E zeq_O(X;;?qgW>%xwLCyadIbE{gLN(Vq-}vI<&c+0_vZexM7`_t3OObu1uTE+5Ln1{ z7C%=~gjYLa{@N_^B7@Ex+Ad+W+Km}bXv%4yRH*^I&kJkdQn+Eyic{xY)}repH0?Sl z?nx|^JU0#FTQdzlug2!_WL~aKHk=p`Odh8sy||Bv6;?*Atv529lY_37KLS$(z$V6&v1G428pGfNI9oG=J zovWenv~7V(Phr_m1hb&!;<*>NM)3z=yiifw_~|={5A=zOrnNpK0PIlLTaXQPiZz8| z25#`-{^r{FM5L@pHM`icMd$vyBCEaG{E+aJt(-H*H7}$b9b0?!xG7TFsbN`C24;49 z2RNN2%5uBMJi`$1qFp{Dy4k}UM00SLR9nc;_~nR;eKqShaM%-c$>6~xHcZudB?V|Od}n3w4AF9 z1e_i_=5{aIvo8HC{D?U#5D3%mD&JoyE&#vzq%YYxB)kG91Ac?c`gH{kP1nSK2||%{ zL1sq@Vnp+QWtb5~*#N<^j0Isd7`&OK&*L z$8KsK-w=MyL;oZmd&w0)<^T1m%F-3jhpsj>8mt>iM34A@OoRGyu^V{ArSV<{SK8I{$#D7V7boB{F)|;HGqq2UjXs&b1J9tP7o_l0FMZ^-2 z;!KzfoYA#G&0>EBjphR$R>+q2nFJm{?Uo@3i1;Ry?xK!|@TuoSG| zI8>lD-7MVm;Gtb>MOU#d+7=YqdxLRgCpkDxx3RLi34r{N;&ru1v=1MS&pl3t=7nI# z>cn?);R#mMHvN9t{$qbxv8)i>hRPSRgOuWov488@>M|V(1n)_mX<;bNOf;ouUu@!qk&-*hNm0h@&GEWI>|8UOMO|7gb8F7Vx)S^X% zgj9CxFwqK@5!U_HRdsyDzeh6gQmWVvjd=KEk!Dw+JSAFp*KIL%c5G+C^m#m3*Vk91 z9TL3@NjjJOXek>0KeYgQiVw7Qz&R7zwoD$S6gAq%LprZz5CL5C_6y3oBkK#04uBoW z!McuUuRyGo&aB4I`&0IqB>XoWafc83Hv3(X;kihbo#Z4c-MkUYLe zJkNZFne8UIE*lnU7Dy|zQdT!jXb8A{4XJAGqN;-r%dDj$pDkN(ZNV-y~i%qsY^C z#FGfG(kdHWPxRSJXm{3*FTFV3`+Kq6jnF#k-Fr7@P*y#{LmHU~7HILKVKr2o3r}f` zg3dt0J;o%;Zry!sS;GE)4oR?gf_PpRgNHH5^la@`daAThY3N^hq3^`Yiv5Jga&p#t zsvQ0EpB18lU}QW*y((Ep#YSE+dVu%s*&lYI9>Do0y2c7U-(4Ae9shdDQ z?_k_-4VkRaQ9k9+Q^|SUk8=KWRlkN{@62*8k4Z`086v3Crek|Lo&!&!&53ppKHVL> zga@C9Dd%LR3Sx%UmSYh_!S@J_&=RWEyNBv67_P{Y~a$6IYv`z&y=zv|QuJcfr-t z;>ne}!k7Gw(|3COter}sq&7JV&F?W$K6$Go?w}`Sa00pB(HrNP<~f8eB>IT=U{~&Y zX?kUSbi)E7xz*quj0X2Uu&nQ}A^6AVO^WwR>zT@=tfjKrSZ25irMie?C$kcRnK>!9 zNF$+$jiG+N7tJm4;{jgPc5^_V>)G0N&Vc>*wZ`&y%IreF@Qe4n3|nUC!B z_1^17ozFAFb>yBnI%!I(y`xZm@m`O78kO}&C8jwB?vDwy2eizPpHE?#=NarFOgazF zV1sDqyMmEN3BzAWG?)aGyok{_&%UvNTWBL)0Tl!X$iJJ9kexI7hx50Vr8(QfP2tju zL31-q4fYHiM(?}b?MgYLOFvw|A#GFDYNxCmTSW$Dw1Q5kQkeGnM@Y!6U02YZAtUN)lv*v^EVlhA*S{&voFuR>WNwo%}IZ zEO-;lUjA_UpKnv#k{p63*EN($^SzN_2z$!iA$rxD94i(0$~YCO&PS8o#~z0(D=xXTn1)(6tgT0RvRQXU9TJN=Y!E62*v#P{<4&Pe>}&Jz+vJ3 zC1Eq9f%5wE2ZBf$wtGeW#C3dEOQ0>`U~g%0ySX)81?Kxeo1otK8kh(;Z<->uTM-l9 z!^CqL81bmd)eJ#p10LkK4l=)coLx(}Imc6dV?aNg`yZKz8m%=TcH3_y$^6c5(GQn@ z*xDmXHY^^D1$rns7^{MO=MLUpiNuWesGkbvRGsb8+S}PT=7Tj)4}$?o&rHkBzq=%A zivIsuQP(y2V5BcD1eao_H#JUK_&D~8AsLX^2sK0tqPC@_O!MtORe(;f+y=<_RNgPC2+`E!0ezS=f|^t z@&he<3^l-Hpem*Ux(ifbOhjg#?Y0xFOe%njdb6`NS&ZvIA*y9|@TBv{P`KI=jw5^A z2S@HhJV(TfOchwkKh|onq-h4Az5WRr4aHuv&PbRa%zSLEktK{g|G^rK<^vE0crXez znO=Lo$>&%2rxJeZ)C(aQZBe&xP4tg^Ygn+~7rmHjXVdzfbp?VnGF6cDMP~r@Q2`=v)|JZ4-t32ATN?EX3XZ8I_|!Xn78}{2ys$jwR6T(;ryeqj!IiAx90&*Es;%#Q-nO|ARyoQp^ywb#*doocNb!q#JGTW z#XWTYP8(yOQY#OC|F!qEbmZi;5jNh)u~34w;=k6Fslo7PA0SfWZ*h>!xbor)Gs;N7E{4TYhCR2zs6X%7k!6j zjnw-19dat#m65E_ zer5LdwakB_c%mly%rqOTeyDWcX0<9~mum9;urU9qKe(e)vxTrZ5>G|@sYwc08bw5* z#K2O!5&IVZWl$|o;oF7zRj`h;)VW7-gZ=WsAU_hj7Y}^>tW#4hy!sWY*y@3ftkzV! z@n^&=O2cSVT6l?<{N`K;@f@*cOQL-G73zfuSlsR_lLyb*=XsBAssvIcW7FUoVBVl_ zowljuJ=eSq{l_M2Ls~8r88fsyfM^zH$XT>rnnQN!+wP@kfb}>X?q=48ejrsn`E1{k zH1x(29*G(p`H$WHycD&S*bPW-e`(xzO^Fyn&~dnJHsaze zOA&k93Lfe+Q#AG)(to5HJqS4~b|1l5Yuf+dF}Q?iSWBXS{ZLa^@!56sXro9@t0Foz z!sR#GE(38LO1QpfN0|foIaRnc9g{tExHPq!t~Z1nf)%$zzBnusw99~K#l|7?0K8E! z2T}oc2}H+{;Q)eiZ)?^%CTblJ6)Ho2^&_{C6`Zve1g8<0Ir z!m;bKG#QgUDzJ=b+ zX8Hl0PZD^D77>An%{ko(q-SbB5Ywzl29$q$tU(N#aB%!8%{XBhgjp088s$TWV;zcK z9LSfm=_E`P6OegU^O-cquj#Rdx_oNpIfJ_HY!BA=>5eJ;1^*nomCyB`(>fTe*dk)) zcS>`)qS%@9-;&F)HRYb}laOzj$~9{3S^lQIf3NLpXOLerIwQ6Vkft-hQujHMnX1PFhmZHCtrB{u(Iqus1@AI6V>875Efgi&e zrUH%bY2s(*d0Y3J2>M3R7ln8d4y&YjwLb6!s(mK8CJ7>?_N$-IwwS77Ng|SBoWKFl zo0awA_Flr&;28Q4=sjZwKDoPV%)?ES8lXA+-t|=w@^7sc5R1 zSmh-1L{Sq>?3jD}BmD*aURrk{MiNs)#sJjG-iMO^L_;IE@re@SpI(DJccDWO$49XF z<+rpnRfPO^=A!{0D2VE#PT zD4L@cGjhiP4PX$G3&U8!KHFJaY{Mv%MIWh(Z5~aFg)Rhqadj4>?o+gLSwJguT(FLz znN*1>xF{<@^Zz1)w=<#S2ES<(pOHXPYD(*%Tmw5e z#Zj06!FGDlz&M6>p7Se=;8>okK-mq~Ho87yNlMrhChYREc0u}&t{=s5{yBD@#Wkqz zwq(>9`G{HL76%sf2Fe2OBR@=WUkkO1hvUgksxOC6iAk|r5bs=jv?C)g?tXju^M zb)o&_t9$WIu%dptdl1zmJ>x+Va#{fnJXQKwi5`4&%kx5;qmAxzRVV01#9Gdu*9#PQ!Be+l zv1lsJ$d(kqD6}@q1|qn8!u)*Gi&+~z@8lPC-nfas0=wBdQQrCoMUUILULWGMo~>`a zK20R`36q;WHHeki`y{(1(V+A_;iXITXvWLm7_s+l4Bh4en$Sv=C#6Q3yM;!c-Adx7 z`*`*h^KCX2jM66>TpGFk99J;}UCSnC$o|GOM|utqN8|U}NCkKat0)1a7En(=LYy}a znP!q^mAyNd2eq)gV6b9mfc{p?aZzPRz`paybLV#?g8@#(JQknsyAvI^V&G^LhDC`> zx_%XpMTm;*i-Rh{IVSdVewV%0&W*6*xVW=>1s@HP83M!|Loj`?s^F(@tK~qnw@7oN&VD``*x9X94b2b}^nOrbEl$)z}$oq1ADIsaFh$Rp`XGFV-01a)cTN(#b z)aWoXG-3DszIhh6Ft=cErodckLdK8QJqOAZ0offBBWOV-I!e|g(bO?-dCVE5pVYl_ zt{U0d@Roj^x>+1UFSYPYDE}R!*FqhGDU(C=S}^tmlQ4r7R#i-_V35iBPbQ(d=%&-o zm_C6A>fv?dMH>o$J4*){#H65};3B-j6W1`4HAR1jI;rh)LCAvP-4>DlBB{~>V@+Qc zL#O8YcGyR@xZISxbAOS&FsB7v5n@WKj=Vkc=4fS0>8}TX`H>RssA`$?;h&wY+6+wU zyPs7Zo$T}|qh-y#(p0w({A@If>7( zmpq8wZ}6j28n}3|=X1r8A3qPq@4DYYpA`Pi9pNGhyZoJdh;xiJ87tZgF1)n#na4(B zASNL(u>X;^EciBiBF%h&JciGnzR*uG^R5Jwl5GDGGw)6Ojo-c5>-H9}1Wk4=)zu+p z()@D9jO^nvwB97o+19+Drd-TfEK{QES#Qahqa}jyyj)qvUr7>E&{~C8jwlu|Nzqv; zOH=hRa#8(rk#(hjMNBa~NYUPloOY1{76vb`06FV#FBQe!OGK$MdJ)5ZFk?S+$)XQi zE0RRtrEUBv!n|Tj!2jq5EvM&9&Ut0rB3dkK(lczPv1eZqNGHuyZJaND&SBqH96QXj z)xY$chtX8Rh|+tP6}e)69j$aK_Yz^Xf7u+l*8NrPXR-wP78<=P_)-*hKZDa7c;^y> z1N>RPm&PReiKVG$U;bx7q$V%!JC4h-Z<5B<&+bQ>ezcu9iM<(} z`OX!EGA-#=4E}XCzXCY6D@Br+j9ar5k7W;Dkk)5J|3ydwab96FG}xQQ1iDheUaT^L zh|fo24Sr}T`L$~^2n?nNuVtw*2>o1TLAi!|ld_(f%bC zF`Gn)YC!5fc)ra(=COu!*kB_ut1S9~T?{}-7G6+8l`_XdME2qn@K0OCHPe?kSzt&T z$}sE4yd`YxA4&f_nW)$F#o1czcegie3F z?eHuJ~?5m^_M* z%jV9i{|R%9IKic?R_Kk|4FRluC8+E~fJea}Jk)E?CoTv`nXj%I6Wy29v^ch+C`BCg zS|}dYQRH>mCPhTpJd~Y&FdFZEaS0Cop>OHsk>@1_kslzAfv{7I-bP4D+#BT8!tnZb z2J+0;BHW*8+K2{D`Yvzzt_|$#FdyDL>QnNYp!Il*-Gdn%UZ}3NAh)?sXh6Z#T&1h8aRo!_L_t>#qZk)mYv#ae@!onSqX6G5xN4%bj2GHIdZSUen70MD1 zTe~xvtl#rf(4Q=`4iLq*gHn54^ZM?dGJ{&4!b}#a`(FFUo}rDjf)o+^yU{m$ZDDvK z^po+%(B_qmnAF!Ww96(hVym1-VpDsOu6a!amp&QMc`=yQG4$pRF;|Qf^%Tir)qu&6 z+BOCD?)89P#H~xLPLyM6599vDOss1n_MdwB6@p$HK|)Hm-_~_4G#8}{Qsi5lA_&Qs zk1wc?b1k`#QxCrQN~Lcb1y{wSpY+Jhi1iu%=0d*JoH(His{S(g%IKgDl>hjvqL}os zOWpEtC@Q?UHD)@}<;Fh!Ry6Qy7dklUj7Rh-#AoK9uB=lRVF0(tZ`8A@ux zcL~H~+j5tMVXxkXK2nvrbElUxXYoeEkvsoxVI`c#=(6lFgU*j}`BY$b&xi216OVb; zH~5TL+?H0al}1S#%`4!A((S$PSFlBGMhrL(QgQn@FZZ_$b2foO<197hgkVJ+xTh5h zge@Oa9b(qI@&~|dNw|Ek*x(X+(v300bh!=9PCybdTXMcOT1Zo>(HwkbUig@uD>(3m zeLlME2NzMYJbQpov;Msl44iHnb3c;pOWtG?d~O31D5ZfM4nff?;4??I}csU1#8Wkf-N-bdYoIK zHNR?(W2wg~2cr}s*s#l%z%v15$&ktH%q{yiMRd0q-{YtkfBvw`NYeH(5#q^ehW5vt z-+x1$ZO^CpP6YEd{>n8{j;l!v#%sPF_*Db#-f2UVC7EWR)rH6MRLRHiSWe+^81MHq zP5~#Jlwakf24d#-j7cT^eKo%sXa6PhUt#__swyCcB;h0m;z~F#TvJZK zRTuXq7x*nTSIUrM@1$knkHnqcM7uv7bz7L&0R27FF;&)Q%;^Bdpp$G2Q&gSx&k1>> zh^fEB+0~}gj6zki^D-9Y6GI^O_HfY9rDU&7%hP|=?5fTh5x%V_s%f~s;@fg0$x3)a1K|nl(P;D^+m=blM>Ww3{rdA_o@XS1UO8s{h0tH|z_T|{5Zuz}o%!_*l!ixcv zrE>w@xw56)%f2L}T>9X{T)#jIYuWoOi$2j4OLJxs3pmwf{@l!GXB?5ILlZnn5mx$o zoJ8!HKTe!ZVnseCDs1Vi0wnOhenAN24}YkcNO*DaY-~i1G@s$o6~%=^Jl5Q+xsvC5 zk+S1TV$>7DOD{o(Z8Uv``@4K>^3t6ML6a{ z)3XtwWGlNTu^_7tuSrZOz0!@D8O6*NI> zUqi(Ofxu^4xa#J=X{Pv@StZb5A^0AzsJ~zq^^D%#I~IwLB!cWNxUDu9F}6HKURK+Q z?41~Qmg{2?iKx4qp)S}altC{X3ei*PN4l))d zp^Wwlb>Ua~>AMr=Tn{-+*bDzX<^)>O-lVN}meSQx$~dX#o8MA4kgAlQCSCn-U{2JZ z*Mycs@#bo(CzxDvg=zDyf(v8khsHn-H_NJTG-CdX34+!1U^hpc(-f*No!7t2i81py z0AW9<6HLt_=OCG9;;Mhk)fF4^2$_ESNvZu^$6X*#c@2L*sObHN{NGik*Blh)kC6(O z5`%vyvy{awngd-MFKoKZ`h>qJti3)E!9MEr2JNQldx-|{RCMrX2|y($WGe6Q?)f-o`zjgZ9dRdO4Qt7#8 z%T{|c@%`aeFWgesHqD7PW`K2Kq;t3NvDLOTE+y`d*^DUHVU=fU+bJYv&n0P6BvwuG zSUhtkqOqWXGy+D&V3=qgtOg_`tY7@g_lQ?FJcYDLC=@b9k+FK*u| z0RJ3B{M8$8@7?j`aRgQrMQ@G>U%wSG{Rvn5H~i;(XAPD2#aUr<=lqV6j>K{I?@QEf z718yZ*46ezrvIr6$2N5EdHEIJt&EN`DOZG@%2&R9D|~Bv1AK|#!QV><*whDI!b`Uw zkj9b{<_!FQC!T`1tVnO05~~HS87KU#!+!ez%XL3JiIaW4?5+IwgjaYz2C(OZL_=nP zRKA(&SOcOxiUrhWpgI;cc^0E>HcCyEJp5tJ+nawebWs7dE$=sy=hfSrCi~4ngP|_C zj=XycT7S-4qXw-A_~Eh^_VO=;JpmWuyriIyb(UBr#t`$#5RJJ)`d+xaPU^QPmtr1D zPvQP!)`aKO=u8n~FYXqvD32*m?2sN$lb$ENaXD9~hCSUId@Uo;t*bImiLr~3PVOcP zK0}4fVgdw39s?FAVtYHRqf^z(^U5n(-wtH=OAFzLSQsN?x#pw+sJa#}xkY`kXb=gp z)Mo}|&_8tcLM0E^OHFdtNguB2aULJe%h~7!r~0_Wt_Io5{>0on2pWFN5eq%mK0x_$ z-C|erZ&8_VH}jLx%DcFCw_xjJ*vom;(o(g_wMKA9HZsh8G$A=8trP1+z*z^hvR#3& zXQ*m`?LiBMI;YG2mV|fTAMTl` z9=`tZ_&w?}QG9ZNm(e4z$w`{dd3Di1e3|GhgZWEFubP+wUZ zBrXL*RfLXMgs_E-FM6vp$JMklSyk0{tTr;Wnrw8 z4HBfo!+lrV6E~DlLiF`S^^~m9^_+fAPx(*AsE`Gp2Q+)n(2TIA6%3ywh*M?5U5!1? zhY)SKoS!x4lHjsJVuRq;1~f}_?yrpuX&oK!^UkFZr3yY-=iUC6C;nrotzZYq$9Ya+ zwC!#erFm<3ZBpQiuW=HAHHAA@npUf7I&E>ft{cg!o}tMWIvj z>cjkpI#*&$?f7oL&v;6opz(K3mdd=Ye!V4%tBt5z+*V~`!?m8+3q~meMJ-_khG-%M z>#iE1te^9Sr0>jT-BMmq)Jtp!{p~o1GWH6VwMEu*hS+mJo`|ixKgfs` zg;DD|AgrMKB;T2Z^07uKu;{W6eFbw^F4g-)Pj7FF6XW2OZl*s&uk>jy_vd~tR22@A z)1!?nne<@?l4gyc!h{4E8sm!oUrlcv)>iYp4w4N%-mar>p8@B3a?{yBRxb7uGK&d%I>@19o}a(Na{OzHj6 zlZG;{&M{U2tq9mvnW`qPacm}M2HVo^$|JsHq^#zKZvG4moRa^+YSfJ)`->4+oM05= zXX(7MMgb)Hi5d9iL|@B1h46p>{kOAGpuw>H6aKU-iTe}Dr6p3#Fe+>!8qD^@4yW@m zzmBGKLdI?VGYqqvHL)PGBg$xLcSubEXiTAj)b-FO9>PIWq6NQlQVkM=0~OVju`j=z zFOI#l{&3<*^Y6`uN*@+q;|PuNuTa`cW&CVzNOF6u(pBvyc#k8DE6LPE*Ce@y^c+{d z65IHe(Wt+S6!rTaN%(-z^YQHz>JSsM)e#Ne>9Js6Icw5S6w{5 z>NKmI`-Vkr>ZeaUsO(g9^SJZy1;5j%dG6}2p&jL`v~qOzj#7p{ zq{(TM{gHGBjb#hNODRAvEod0%zN77wwElDJ{pHIKJpTazsuaXLusO#gyf4C8?CyL$ zH3^(~Jw*`5{TzrxavyVGXU)W%>N6sDmCv-)N+pAaUxV5^6h>tZ_xG2=gSM7`Y@N!U z?k_hV4Hx9L$y~^6$;BFUTAqtuSAN5>FCGMNev^chK{obQh^Askca%gw z6UN1_52v3c1s1sUJ7#3&Di(&ko@tZ}D3#f)Ui^X}W!Yxg=AYBAOGiINb$~8*GFNS+ z>kytV%%_;XF>VKslaDXJN~VD`8n6W=0vhudolKqa=JczQq%5xdl4?TY0r4SbLNJOH zF~R)^75R$;$=)8|PtJJSDS~=TjcC-dZE+a?G(1;gMpZ*A!~k*Od~;H^{a%s12Z}$U z1ORfxdDgLmT8a-ig(SP_4bx)Y0?X6mK19Syi7dukT!(y#fR3N*!hxFmH$86;W_{i# zQ)GlmK}ASjZXUzEe0E2K&@*v|?7>>KlzcTxIRccos(1_2n*55mN?Xx31pKvY#(j~a zxL!F8y;8I()2#svW#d4Jp?qrO^X@+x@L9j-7hyII#Ubx2rk5lNkqA!xm+c=K_;sTt0XZs+W!HFelo#E> zS1Ll~2Jm49RoH>jY~0r>D~B`2@mv_b%=29hJ@{^v8J%x|vii97vHuz=QdV|a{cLe1 zqC4DX^m}!Rk#-WyQ~Ma!p7&2}2yMBE7(U7!$v!2OaIZapgL8&lPkNe?PYi)0cgnp|#?Z zKumWbDW%U6YIqe84fsL0tXF~9GFdjM&HY@hG#g#|o4m?TF$fyyiu3&bk<|jA--m&# z7#K*ct+6P@Iqm$@cANuofS2r6wN53?_B50`SYKga9?T{9+tZx4%Ic31%7c4D%)}eB zq1Jq^9B-y#S9iRX!6D%7i%yHGkB@>0FE!W>7JLf&RK`UEzF31xZ$03-JiMPMUuS>+ ziIrlMlmy$~VbSk@>6Msn+vPcN2wiXvD9LPjNy^|Kc7 zV|c81>LjLn=mg&i}EoN4G7sAmW6>~VM@=H5)HFMIb5l> zr`~Xy-r)p`tbk=w(kY<2(1>z}k*UYfM5hC+bjNhzavEpn`gJNCwx>VjfKwsn3u{x1 zsB`EIfY=W>IRUN5+rwD6^Z6CzWu=c)Yh}HnL{uZZn=};bj(*6Z-$z8*ROOA;t^fjJ zihlSibe8GqC;iMLpOH%|&geXwML*lwE73>>Pz_KgS{{9(Sr)%@H-c=ezxu?e7!cLy6tX{2jfNbn?11Nj$o_j*DP3XF0%W3NlyuT#rea7WP@H-(!;vnc2(5dO1s0iq#;&(f zSuOq8pV~DYUhSVp*8n(pbd9@B0_RrwQ!gyyL0rhms-1AB==~n;j$$Je+@+naXUklt zlIiA_DYJsd2^C^VUCIi=uKX9n&efw7$`8vy21B|C=iIr{2~XPhfe- z$IKJN^MwAS+fw7VJa1i2yI=lpg7XK8sF>0V-5cLN^9W_&B32VRP4bPB{t+~D{w&C; zL5qx(PjM&a*DG%fP48Cl)3;!b(8RR*E6n+9o2vqMEJygW$kTX!ov?Q z1GAsb&YV8aN#%*AUih5C{W^8KI%MrBxM6O@p)j>p5<}MMp z6(I%<+bSE9B^b{gZX=-IQLQD#z(qDt`Th~<N-;wgN8_H8_#AUNgr z5P7Y%kcU>jggf?bDhFzbx)r^1Iq@5b^VfhjRC#WOS-fEW32oNpv%=yAxs=Ha*}ma! zhp21X`oQ7&qD?5Q4IKonsCfO&9IAz9jclA4Bs#;nBmMh|`tz3P_!(^JSQ9Hbsl@@& zZ5C+*2s4Wg^2gtHMMI62uZ$HA|y_mO5MXg7AIU)rdRK{lirFl{w1 z3R{6T_Eg~jrs#%pIvA@|^Y`Yfds-rvXj>=Skd5(~XU)0TI<+md|tWsiI7+_}bA2CRjjq#U$S{{?+e_$m_{kA9t?xBu{aqpV=9 zgv*;rD8+Zxe8rjRn9OdhR5r@t`HcHj6GNAO3)+*WRCoPLI83w-dL&v!0oRz2W(U?& zYVA%~z=)xj!%u#+2sTs z%sSeYkxs>(6`Xnt+}lfY_${Gg8jNHIH7d|+O)-%ZHA^026f>((+Pu*0Zn zAD|6$i*fTcC%t{_8u+*y8fMbWv15qp0H_m-8;Q0G=n8+f`<7qZhf3Z_XX-NULze5M1y$w{=K6B%D zNJ``Up#?q}ym*VqbCTsjHS@5a3gm9$MlsPa5Zsyw%PONS1i3t(WI!Sr~*2` zo+qjvsr}GyURbqmocsnFu5wso%3F6FDf_aU7w6!e$>V{IiR8ADy^V{e6eA~KAklK# zWyd;6$YY+eorJnHJ@dqg%B)<-()4FFAr|0d@VN>yaG)Cy@XF&UMDO@HY6EHpAg-i< zKY$()D$gcf$>7C_`C3#3P65LD`~yVO(ir@Mt=?G^A?ti7WST|B=VhY9O|C=` z;^WMvbg3vg%Q>d(R zV!d?j*(6P}*eM>=E!j7Jr#m(L&eiQ0VzK%FC#7B3a{MKirQ z>(jZr2W`A(0?YxPbF?&HJ^=$ia&U4uYI6kJ97n?sv{>_ZWFTOdI=+61oATrKu{BA& z$aB1~57#-PeJzqC#6xW`Ud6KLz#yf(=s2LU`S?vKlui&1`(X!Cx|zBzf|8i%u(aAG zdl!%x>qt>5Z3)%%jgmpnQAH$G3omJZyrH!`2q9{TS-%+VZ<=S_emW;vpq2ieis~3i>oHPt*NGbz;_~2 z>ed|EP{->{u5HdHd7l(EyoWO5W~$urb`h~dY*ge_Qv|SP_~ZhUkGx(B*z+!`*fgu#mplp{nsP?xLBYY2L_-|0cZipGWWim9iz`ueVD2D;?$L+{E^%UF@ipK2 zo+h>~gl>9>d+4BMZY%63&e^W_ZFE=jU+}Gl>A6#PPg6ehAc(k;s6RE;&)>nyG8e+k znJ=m)$6CcK8Ts)>M51bV< zfA5~xdXq6J^!`^UonUUTf%C2$#kLdzDnl z`xg6>FKN}r9lfn^ro1%bOya5Ctgni4!tlN9`7uKKyU1f{aCs2#YGGyo+hgt~&W1ZB z`>>elgvw5oZkc!((ClCyY{@NHKyW5GrvEX z>i%;!50#fEQ~8a{2^#&@HhGTCvTYP~c*kCVb)!97GKEP(k%xqRd`3|zn_@be+ElrO zJA?@9Yjy%vA$nT+eMRu(JZ>)oPYRe%j z-j;k64rf#Hvv`}F=gq2@x$?zhk=zLvdtJe>+F{(*ABHDw<`hci00sPTY`$+M1w^e} zbJ9lZ1rA>bHkJ4TSiOYW z#s`(=z2SgP6^A_Y)tA?UK1#c1@9Adx2MtfAF8O#wiGN^Y`x7c`ZPH#??*_=82kk^f z{nZH#PZ8!V&l#V_OVk+Q*!*JnAK`vi6+$rG4Do){d6I(7 z)Dz>Tu&=d7Y4-xaD=Ie5ji2(_f7hBx%U$g$^FW$%1ZjbWv*5uC&%Uj}&A$)e_8~Q3 z5kMDEgVTcZ66*$-l7NoB-nZnumgp>aYH4@TXkEU|B5Ca%2AHJxrqk5`$>GuwTBVM_ z8@pYy$N(g9+f~{lLpj&We`!Acc&>T-wFq2v=4|iYvE$g=pPR6hgF9}6g5FTr)&3)b5QmU^Ca4%<(oMhK&y^s3^=TW zW8Dh#dDZDE#uKKbq9bzYwyC6}k>V_=o{=27;mhNUgyMMTeTOnP`X}M-fM8$~2b7Hd z4cviKB)-cU+K;#yC4}5}EkM)Rf?1?&&O1AlATucCCo{bk(iiwMU#>02aBR9|E1C-^ zQqo!Shzp2Xn2Bfii1$g~0;Cx~$BT*yKHuflS3nWT)H46_4;t1@xHv^U)HPKzMjD*_ z6iRNxTkWGqQRe2q)kBTX5i|6aXq6`Qw$g^$>CIfB+4!W~?VLKyR2e7{Y=lRp`*Gg| zHc67X*{0WR&i$upaz1)&t-e0VwByd|?GP8Q&X!u$S`#s+9Z;yNi)V}?bDB(v&178X+gi#RuQe38H?}vZiJQu&?(=DZA!C+kR|l4k(Z&ZzG(iJ1xbWJY^6D{3hWIR!OW}Y6Vo(Du&@p zGF(QZMLdZdL+ymX7th3*`&9KgLn-C-kA<`Lr?4f-^2M36XyvI+N3A(QZ@lWbHk=jr zdYUG3-pT50ItR-E7sZ58-`J7oF#%87nIfnjp-PiY9V9yCBK25llsC6nmkpk~_aVcex2BS{M;{Xo#v)U@Li~G5^T61XyTMG!x=tF7xKZ?xMs*$Q z{GA4<{A<-xeCgP}P2*m_f_+t9_L4-1aF5;Z;|0eKh1Woi%`T7{kAhl)HoX)F|=pSKOubARd6X+)0VG z3qhS61|ey*os+Awqa=-LIgyMys$=f9%1Uo@?>u_?`C@kNtW>~1y)BNNZW&J@`%vd! zMTsF1i}@BXps7a}q;JHX$GL#SVbG6@wgfrtj=j&{Ua1IK zgTcadp;%acbgl6TV+nJ2!91Ihxf2CzL51*yBAhmY#nhch{-_BU8_OGh%Z0oYpkmhV zNr71-yBr$}dI(hsWr|bL;5|fp!x)|TweE$djG$&u$xjE~P^=Zz4WHbT`xfkaw$oRf zsnhg{2rFm^J53s{xieHOGWH{NQU}F06-q;@FOWdIuQ8WU`^ECMvn|*lhfed=>|@t- zwjT_ndJ>a(0jd?Nn%}OEsS%Pp@Soi&8$Pz1x-iK|L z*)5l(MxA_Yv$YlYwEMgJ^H(uGEB93m8U40&J3Jxm{F5v4{9L;HHy;A%BH_|K9bH4B zti-9_1vpClh^8IMicq{``ggC!Mmfty2i%Kl(k;Hn!yPA29_MaMQJRx1Iz=5wLD!Lu zRa76uV&=?EJ1*(WXTPO7$IefBvC=zUSvzl{%Yw3?$*(w_&KXH|?>s#6lO57m?2b>{dQ9*wgKzq+bkZbcn`ZnVG zmguR0HD5+W@eTP8zE4`H!CJG8n>tr?TJW1kiS!Q)pptZYNk(R77$aC_o=od}CpGi& zP-Xxr_~?S__!VsjFSzX2HFu5#t+;EQL2Umu)r;8V%@n@G2<_mIFeeR6wiW669F50! z0pK&Z-8KYQq!@VvXEHh3IV?+-OK^5FrZa$Vh)Ka#e7v&vs!dGOI#snq>uNv&3Wg_h z$dt*wEBRCsCu+V!VB<6`lMuM{J8vZSl^SI2RfsM&%jMIzo9$=gL}hgtBtSS};rAWP z-%P*XhYn=E&}W2q3i`ZZl-pnTn;b4UXY1XF)-}KV3%39^$3la;#XHibfAe4drwzqM z^3WX(P|~fRJ{D}!^Y$>~ZAfd#@k}b=)cZae102&i!v*?Y1SxaUw0TjsAJj?C0HDk=aYlB)@v;w6y>vDaA7Y z-}siQ0UPbj&_VO4ob#l-z{7)?li|V^( zGTyKCPK}*EFQ+|Ma1*2VABWr0gH$~1Oa8NnE%p*3n_a}uW3)dRKR{HRS(W5fItXXM zg_(E-|2{dk7u7D6!rr?!(oE`gh%3ZVHud`hxyc-1O`HfUz)@SOb5Zl*Uom>pBJ(`@ zkLFkXXVNGU_){Mx*$rY+tn)&_S}*JDP}EiIe}{Ue1rOcj4ldQ^l?TgzV&EQfBm>0^ZGq} zJ1`55>uH(j6OYyLqBbK0rG`zWnpIIiZE1L8sIh04Te3pTle=|11Y~jtEWD2c+GJ8o zGDjzj%6Hb%*A4Y{KHB#~GMx6+utZF5`XQj_iCsRWjRn2GoIyQKW}!Pc&k7&dEa)S*toA+p_U;EbN#AO2>!g+>!h9 zAVd~Yu9CJ93-s)?&OFjB<2!Fsgf2>+&FJ}pNtw5Kq+HT|LBo@f-_vjs@5>b@Rc*~f++8X#P?8G!Q zB2@e9;j?j=vbjdiD8*YQWx__Ht-l&IEgzAWy#7a|^NYLJ3; zcZ{K<_-n>E@oC$Pc@Og=JfuM$36Lq~(HwFGX$&*+H=`zL@`rrJ!dDFKY6X*4eg5pxH+pZ*I<^_s>(8EE?zGcqQLhq0_S=isK4h-E zlX*4Vy>^>LRFk-`#m@AdaWl*-zIQGFXC7FStj)(NNpdN*$zer%E@tZ;0{C(o7^es@ za`H)hi4$rU`^qtOu^FPUxiL5ON(NQm`9b9oS}(E~c>n&a$2x(q3lF-&DU%8i4&mo? z+Oqrgsk!TEU^kI68LD-kjD)D z-1wute>+V(uV7{@Jj<0-$fo=~Kz^Jv%&rbC_RBxScb^Q9Khz@rWc%r5J=&RrG6O&c zd}aK(yv2^qd~f3N>=4?X5EekQ$hY@vNnH3)1+~{8nVR0!{CnS&tx8zRFu{ozC53^m zLWwN~()#cTZOm8`l)T?Sf>NObdgJQH8EC--pK1UZaU9{B_<+5by^|L7_XLCj(YfrC zIi)&_V_i_q!wU{37C@ZCx{m(<<*-VZFh=wk=Az#^Lp+xU zlZ2gYnRu9o1g_eJbi*kM3YxYUGBEJHk;D|uGr_kpT+sK8%{o?d0qVm;K@#A^HzLQR z;M7w&ilh{D2jZo$_Ui|{(?b&NN7lsoJhnM=fyw|7Wv#l}afpMhGeqxq3Pj=Lt*Gej z{v(z|Jvm@7JJgH~LLJ3=Ns~1gV&{GXhWt18)$WKLo_(&EiCH8@Xt&6$2}O~=^5)%e z{7;q^^@9Qqf0;j6?6`_=ruvD`(8(uD0?!iPwB8H#q-yU<^89K+4ql{}vz$6id}K4Z z*?(Y8|AC+n=k-gk7$deZdomD~!tj6O{S`73TNz`LXDDCs_i0udhF(H?TNRyN&`C!n z@|ZcMxxarc>^-&r^OK_EWo@U}ykkfgxL{ALSfl~#Dajj^v7+Q|5k&e+oPyxeKmyCM78lRXs54r&@-m5KdT87w5LdV=OqJCrS4;Z(XQQ)nnVu6u;ORb=hL-HFOB4G zVu}$G!r(YE^^M+K`k#-q-v)84o}_0UTFB}3zc^(S5=H*5{;UjE^I@aTS53kRr!gP? zw9ej4ADZ*)iO2KZcXh${F4mki1MfPxEbQzJOJ?}GxI;c)D1AQ-^KA+TPMn%sn$T3f zPXx?y-X=Lp%=Yw{wbWhOOf#}Tr6^=$L-ZZ>oD1|bmuHu!EiCn%#F^q6Y7%(?Z(|5| z{NZlT-?fE>ew`_^osGx)4Q~Dtbw<_;!Olex%7yccMsE40LUk&I*KGp31YQTZ-4g;9=9@q zOwq1Kq^Ya32hoa52zF^og(YNe&*AWf9vjl0bBJ$*eICYb0-kywX2<)@t8X%T#p*nM z9UbDrWml-Q_7xtBWUs?4jyWVqIp+i~<7_jb_%RSjgz#!8l5f?GFVIgWHZF+$>d~Ch zjIalj5HHeEF0%QsHajFo9(5nrhS&@7c#-$B;6qr1yf^lSJwsyT-pOySvcT5(-p@9f z4Nn?=6MJqBa0{e0RQ)iK`sqWI=gPTV3rRas947qh-6%d?l$FzBLsO+i1U56rmyObi z7{h;O7YlFPVdhT7%})w5x6*3HvP|h!wb2d8M(-kEi@fcCQx&iZNyjmk zKWHQNJT^AUD&u`bU_o9JSZj*!@1wom&araa08%c?O^1~y{0Jejes0q^JGxK2-h|r} zRpsW;f!){|wyK85CZaVaw3AvUeA)|;W|C+hsgF7|I68_jzcG z$xC_hf+|~|K-ZsfQe^RW#CZ1XRCUfLhi9d%OP;9-#?@{&!vg>5Pks+K+2@Cv5c)`) z)`iV9N{ytoLlSW+gp)qBwz$gzBqy=Z<*fz&Jchye{C%t0?d{>ELiprLW{AO z{|oN&mpAQ@cqP4tu^jChFy4Bcd1p1U&R90^A(Q+jou@A^&P--;xb# z$?L|6cJ!#87|BBRIVsGVBSyA=kjp}JRrSteW=2vs6R1ngr}1u*3s8t`9Kzhh zuyGU%6|EJE{_ayTLB*jzKwbW^Xc47EOyBsQv`^k+C{sn(jWaZ?{{JK}RHHcH3^+-z zf%l-W!5UM6e{~R58P_`Qh>;9a++h5Skt8ZijYOEVerKGiJpX<+_xdH8&hF;T{D!11 z(=lzg%R4(cbbwuY-L4JXc_(c+$rlraJ(9-R!o{r@@A1hK70}Or43+Ps@9eXmUh3qg zfd3qRiZQmdgg21P^F?M_K$pC`=;;(_FRaMvtS(GvL{B=-r(qEJioUPLC=uZmnOy$s z{_=9kop!b&Yu#PTazF13$70K^81v3pKVe1qT097hO$9NOy{H(Ydf4oIu^(mCNrZl) z9n-rIMzd@Q93~z5l@P`Nd`s8YvtILct?lni`yWc3GWLteV;g$Ra<-2HGjoh1mmBsX zkdMR>jBAGb9c$gtaRE07I}roYkX9rB|eA~C^zGlQS;*AC)0*m|)Z?E86b0C+y# zcTOR`eJ-X8o=m)}U^`&u71y!K+SzXBUBY_b>?N60BZl8LM**iugE{zIOB_AIG#rJI zS=2m14qhu6LuANJ3cd@@Sgi0tRCHIPnH|n+%VAuiIt9)GCW1g5alFy)o-OADM@G=O zFbSwu#@>P2?r9xm6JVR{mOLktuBU12s?j$rvGy39Uf73-3!qv{!|DRJGwzSab;SY=3KtuwH53!TS(dW{qlQ8qWo` zRvWS^&SPWu#Ye-UEp84=_c=KQ+L(ieYK2bne}@ex($y0Yo4>G0S6VHYFFkK&AeNew zEm%0T54(O_)Yo%z2zjAwO`Ax#VA&PE(%U+tE830cP629#(K`q$wIS%Y>5LyewIz@g zzm+GIwqhLM-pSjtSPq)D@~8S`CCjOxccs7Q z0*N?pVQ($M)lz&`3=LrCJ+j}wRMxo#?w$=ZIA+-~&Hk33E!*^%rasNk)oEj{)e5ch z=q^!|LDO{RLbkg>kVo#U8)eatUtvkeERwIg*UmIRx=^4Es*N zs}x|(|6OCuxc3Gk+CNt__g5x^7}4OHwD<8%8#yS4QJCzdJR~wiWU^@O8L@DxHw}#M z%*Ko_Ma+exS*)MVn9oO&`H~|m=Ny;s{4B{%%d6w4fUW zZTLN6X1We52;&$w2Nvcs{B%l6z!)^faIXVUic4Vh0GOO(Yu0lQBlF~u4}YNN!$12l z)zlDKXtcN7Rk}v-QDM1nO;4iI-pXXFo`h3BooTKDvaOsZ#$DTA=)!$27laKiZ4~2- zVeB8IREj`7hh`@HF0~p8g4zNbfzDqoO3(e?IKUOcq*IJ|ro<8>XokBlWvS~oLLFDJ z5SIyc<53X~PmHaf)>6Mwn>!lGAYKB3vIjpvB{g`mFW(S@ibG>(3lUD`)H(e^o2fYj z=WqBHB3}L!Jx$jzJct}a&`yPM(~No8@MX(7@7%()m5BbRLs@E=mSZv zV`se(fGUlQ{+2#_s?&<6FXHWHBUX-EWHvp`)q2@0{1A|poRrnOo#2hf{ez^a>rV%+ zhB-BqvK&pe$0=wwR3cRC`Rd*7`FAe|1$BuDH;g-fCVo~C{pKmBe@t(X0I(zCd9T|A z3R3ZfMG1J6L()@j&|Sx9k=C>L8}}!vFX~f3ezu>zf8P z8+KIGc)K~r%dBwYX-d9B@q>~8O@;_MZWgCsIT8+fLm#WxxY2(Ge&$cTEcP%cLWG1H z-^_7idVz9j#O0Mn%7tlv>)Z5-90cj+*#h3_=kWLh7#ZTmE6o&wGz{X4a(lDZCScDj z*L|<>$pJW>m*aYyv&jF|(vE4OJ%>Cz`j0z~+f3)-`lRs0xy>ABw%r2JsZ-zQ&}W>U zj@JaC?b=@7DYC3&r%Bx_b^*EnD(4OU&CQ>VIZTQFKKHZ;e`d<+!OYjagb68{DELU3rt&e| z?tftbLmB5{w8dHs(#R(v<7zrN08^Gn06qWoQ%fhno;XkcigRhSp%QV>CAB4%H0s9C znr^%^$$Ju7(-?wzM2L!>T}k|C{8KZHFD_AV1{37@UmiYC1pBd*ar2O|kto)lUj&OO zy*L-r^IDM@-z_Y07K&}*)Verg8q5*!l9T9ccB*j{Q(0^b``%%k8fF3`sc(GB3mAKvc}(~7 zeu?g|Kh(+DvTlO_mBKtD^eSHHP@>|j8o)&iCTt&epSSW(g{PZ@eeW=ziRt7P!Zg`+ zk6znQ9qSOL5t1WKyQPKX5;DJ85;A0C?EA~NZJ<9ib^TJY`{!&GX!|wn@#^DHqHQ|; z%APbEyh2=Q@ZZO-c-s*$kh}d~J2>tXz_ao?HKTvO|CTUcZTZn&*I%h)I`6mq=eAJ; zKUL&tY|gLlAHmYm-<76PK&>22)*BbEDILj(3h^FzFH_MzC_T!f3vw%#;3r8T*52tn zaiL;4z%2FSAoN8J!P#&|AtQi>SdCDg`ql`n`+;FnPe@xx-&<4ihzt-k`2MPFEKCCB z2h{_N5e1G2m-@|Z-9w9sY$%dwJQ>$a;<>+NP~Vb*why|mgU$I79k>8JpgcHaGXa!l-=uI{f-o`zy8d#MOg;#Plmj$IA7<3=35hLyk$y7#hR-@970 zlnA2ZVh)O49J~%B|5^;y#bqRX!@A>aqFo&44Wv#OwpGK2m?rkI^HJ@0BvB@Ab~Y(v|a7fuBa z>qQT=$?xx3pdO)SHNW|Jl@o}agnOWH%w;`VYtrEy7gi`+XyF+Cf}U>Dg4}-SS35Xb zPf1%aYdTDa>7}0eWNs>8EB-Veq;mbWr>V%k_%(59+?_S=oGj3ukd)Xl^`SSs$mUIOj8#Ftn)SkIEXd+Sj5y?U20p&>$zDEy*!jREef^O5v+vnIVvJ)NRL@%t(h=aW1VI7J_hy5a+U*y= z_Ow7?qiGGFb(lJ!iX?l!x+lcq)Nt}g8r3T{bs~2X{(n;SPO`-<4e!TpO|#-KS=-NY z8E>IuauwdLpD%?b(mX2fi1_mF!qKnLZWy-F%?&nhWT?@?5_R}tvBiUkis+E7T76`K zapA8pmh(Fa0H|0v2#Yw1xvVdOv^68Rjp~$khNC8eI_qma=jMpVU0cJGknuj3z?0?jW2I{=rMbZar2K&Pe3~F9Yu-%WJhH)s19^MPxBK+xzZicoMWpUSuv9h=?#qetnUdh{%3)WJH#Iusn?<0)t2r`ab}PF4 z$09~0BLWF4)-dUmNeDwUVT}*W3$(X_LrM?;b$xZU>&mgMjqb>SNe}*sZ!hB#V*Jce zvHmtXpiirxYm<_dn6RjY0fW$$^5*mPs@NY*T@Ol~3)5^vdhZ>d+s;aO zH&!1pMm5yYj$<}HBg&UQyx*v3*GNuO;PpwKY0^}z(R}v6){{|siNw$M_$+vS{3na2 zV>UfTKah@qro3cibajhC7JnL(2(q20Et~b;sM~+$h5-|?ts@f#Y#jz)2oJr`-r#gi z40~M?>p@Clz>nDWC|qxwF{Xl!Zbi#yI_pXOf?J-$YHptZN~5Ok0B6S!XZ2;J%R|Ujg|ffNa&d+`{P06 zsKn=0n_erc7`Fs`c~TSqj_?HF_`yw;UF zrk(oF*$<|Exs*sQgXuj5kqOn7EN?pgbUp7Rg<*?ZEjT!wlfH|?rNR9YMF2|By>MKv z%2%|{v$2ly5FVR8#-G1^w%OK}L@tw%P{=Dxwu?)$7@@1YxWj?r7Vsu~-^9%?|Ew>a z`cp&RvI?1#_wqGD5>XqT;Kd3D880J8olj`*oP|*x{jFusLC+8SkcX*bh=ahLn7X`o z`>*e%P2amO^;kwnC1be99HkB%h&QD(3RrwyY2`dIzQC69)&78j%>*?3L-o_#?9Tro zL|y~1+RY~r`gqF;5?BgMP|?wA;iwZQu&_EdGhKZt6?*#|IX4>z#TFmAE8-iWf~;TnR<^1%cvM}W8b z%(!lcWMgzq&m^*3Brzh)ydhr-nPzm9Jfl<@pdeGPNf_BCdm;TuY0&@bU?7P`^|D#q1Ti z)IZZlmn=&-@og$FelFCYGP~ux!BK6R;~i;d9-${zhW*7zwxMH6Z)p@Jp;Q2#ZZ2t% zIthQ;;d|KWY#QkCesjy>nJ>@@;d65*vrB%OL70(gJ%`zYd`R%X_1w!Harl9I@jwj?(oY@~EUcONJ1wL}4^Z@u zlF&@#1|Ltx5^J_BZ8TgdWNCh|livOz5ffpZ15Jtbq77MX{y|WfYubMU!Aa)l4QuEv z1^~c>ORBNNQ~$C5YVO6W{NHs-IzKeK`|tg%Sn5Ai3-c>bQ<<3zF(5C$aKsm^42Z?3(Wrihq0Z* z!s$?H+Wh}BslS3iKwm4hiQu$JyW0DIPpsg=lmM0J;hFU8^L 1 @@ -112,12 +116,24 @@ def ipostorder(self): children.pop() yield cur_node - def iupstream(self): - """Iterate from a tree node to the root nodes.""" - t = self - while t is not None: - yield t - t = t.parent + def iupstream(self, stop_node=None): + """Iterate from a tree node to the root nodes. + + Args: + stop_node: Node to stop the upstream traversal. If None, it stops when parent is None. + """ + if stop_node is None: + def stop_condition(section): + return section.parent is None + else: + def stop_condition(section): + return section == stop_node + + current_section = self + while not stop_condition(current_section): + yield current_section + current_section = current_section.parent + yield current_section def ileaf(self): """Iterator to all leaves of a tree.""" @@ -211,7 +227,36 @@ def __repr__(self): NeuriteType.undefined: 4} -def iter_neurites(obj, mapfun=None, filt=None, neurite_order=NeuriteIter.FileOrder): +def _homogeneous_subtrees(neurite): + """Returns a list of the root nodes of the sub-neurites. + + A sub-neurite can be either the entire tree or a homogeneous downstream + sub-tree. + """ + it = neurite.root_node.ipreorder() + homogeneous_neurites = [Neurite(next(it).morphio_section)] + + for section in it: + if section.type != section.parent.type: + homogeneous_neurites.append(Neurite(section.morphio_section)) + + homogeneous_types = [neurite.type for neurite in homogeneous_neurites] + + if len(homogeneous_neurites) >= 2 and homogeneous_types != [ + NeuriteType.axon, + NeuriteType.basal_dendrite, + ]: + warnings.warn( + f"{neurite} is not an axon-carrying dendrite. " + f"Subtree types found {homogeneous_types}", + stacklevel=2 + ) + return homogeneous_neurites + + +def iter_neurites( + obj, mapfun=None, filt=None, neurite_order=NeuriteIter.FileOrder, use_subtrees=False +): """Iterator to a neurite, morphology or morphology population. Applies optional neurite filter and mapping functions. @@ -240,8 +285,13 @@ def iter_neurites(obj, mapfun=None, filt=None, neurite_order=NeuriteIter.FileOrd >>> mapping = lambda n : len(n.points) >>> n_points = [n for n in iter_neurites(pop, mapping, filter)] """ - neurites = ((obj,) if isinstance(obj, Neurite) else - obj.neurites if hasattr(obj, 'neurites') else obj) + if isinstance(obj, Neurite): + neurites = (obj,) + elif hasattr(obj, "neurites"): + neurites = obj.neurites + else: + neurites = obj + if neurite_order == NeuriteIter.NRN: if isinstance(obj, Population): warnings.warn('`iter_neurites` with `neurite_order` over Population orders neurites' @@ -249,14 +299,28 @@ def iter_neurites(obj, mapfun=None, filt=None, neurite_order=NeuriteIter.FileOrd last_position = max(NRN_ORDER.values()) + 1 neurites = sorted(neurites, key=lambda neurite: NRN_ORDER.get(neurite.type, last_position)) + if use_subtrees: + neurites = flatten( + _homogeneous_subtrees(neurite) if neurite.is_heterogeneous() else [neurite] + for neurite in neurites + ) + neurite_iter = iter(neurites) if filt is None else filter(filt, neurites) - return neurite_iter if mapfun is None else map(mapfun, neurite_iter) + + if mapfun is None: + return neurite_iter + + if use_subtrees: + return (mapfun(neurite, section_type=neurite.type) for neurite in neurite_iter) + + return map(mapfun, neurite_iter) def iter_sections(neurites, iterator_type=Section.ipreorder, neurite_filter=None, - neurite_order=NeuriteIter.FileOrder): + neurite_order=NeuriteIter.FileOrder, + section_filter=None): """Iterator to the sections in a neurite, morphology or morphology population. Arguments: @@ -272,6 +336,8 @@ def iter_sections(neurites, neurite_order (NeuriteIter): order upon which neurites should be iterated - NeuriteIter.FileOrder: order of appearance in the file - NeuriteIter.NRN: NRN simulator order: soma -> axon -> basal -> apical + section_filter: optional section level filter. Please note that neurite_filter takes + precedence over the section_filter. Examples: @@ -282,13 +348,14 @@ def iter_sections(neurites, >>> filter = lambda n : n.type == nm.AXON >>> n_points = [len(s.points) for s in iter_sections(pop, neurite_filter=filter)] """ - return flatten( - iterator_type(neurite.root_node) - for neurite in iter_neurites(neurites, filt=neurite_filter, neurite_order=neurite_order) - ) + neurites = iter_neurites(neurites, filt=neurite_filter, neurite_order=neurite_order) + sections = flatten(iterator_type(neurite.root_node) for neurite in neurites) + return sections if section_filter is None else filter(section_filter, sections) -def iter_segments(obj, neurite_filter=None, neurite_order=NeuriteIter.FileOrder): +def iter_segments( + obj, neurite_filter=None, neurite_order=NeuriteIter.FileOrder, section_filter=None +): """Return an iterator to the segments in a collection of neurites. Arguments: @@ -297,6 +364,7 @@ def iter_segments(obj, neurite_filter=None, neurite_order=NeuriteIter.FileOrder) neurite_order: order upon which neurite should be iterated. Values: - NeuriteIter.FileOrder: order of appearance in the file - NeuriteIter.NRN: NRN simulator order: soma -> axon -> basal -> apical + section_filter: optional section level filter Note: This is a convenience function provided for generic access to @@ -306,7 +374,8 @@ def iter_segments(obj, neurite_filter=None, neurite_order=NeuriteIter.FileOrder) sections = iter((obj,) if isinstance(obj, Section) else iter_sections(obj, neurite_filter=neurite_filter, - neurite_order=neurite_order)) + neurite_order=neurite_order, + section_filter=section_filter)) return flatten( zip(section.points[:-1], section.points[1:]) @@ -314,6 +383,35 @@ def iter_segments(obj, neurite_filter=None, neurite_order=NeuriteIter.FileOrder) ) +def iter_points( + obj, + neurite_filter=None, + neurite_order=NeuriteIter.FileOrder, + section_filter=None +): + """Return an iterator to the points in a population, morphology, neurites, or section. + + Args: + obj: population, morphology, neurite, section or iterable containing + neurite_filter: optional top level filter on properties of neurite neurite objects + neurite_order: order upon which neurite should be iterated. Values: + - NeuriteIter.FileOrder: order of appearance in the file + - NeuriteIter.NRN: NRN simulator order: soma -> axon -> basal -> apical + section_filter: optional section level filter + """ + sections = ( + iter((obj,)) if isinstance(obj, Section) + else iter_sections( + obj, + neurite_filter=neurite_filter, + neurite_order=neurite_order, + section_filter=section_filter + ) + ) + + return flatten(s.points[:, COLS.XYZ] for s in sections) + + def graft_morphology(section): """Returns a morphology starting at section.""" assert isinstance(section, Section) @@ -368,7 +466,9 @@ def length(self): The length is defined as the sum of lengths of the sections. """ - return sum(s.length for s in self.iter_sections()) + # pylint: disable=import-outside-toplevel + from neurom.features.neurite import total_length + return total_length(self) @property def area(self): @@ -376,7 +476,9 @@ def area(self): The area is defined as the sum of area of the sections. """ - return sum(s.area for s in self.iter_sections()) + # pylint: disable=import-outside-toplevel + from neurom.features.neurite import total_area + return total_area(self) @property def volume(self): @@ -384,7 +486,13 @@ def volume(self): The volume is defined as the sum of volumes of the sections. """ - return sum(s.volume for s in self.iter_sections()) + # pylint: disable=import-outside-toplevel + from neurom.features.neurite import total_volume + return total_volume(self) + + def is_heterogeneous(self) -> bool: + """Returns true if the neurite consists of more that one section types.""" + return self.morphio_root_node.is_heterogeneous() def iter_sections(self, order=Section.ipreorder, neurite_order=NeuriteIter.FileOrder): """Iteration over section nodes. diff --git a/neurom/features/__init__.py b/neurom/features/__init__.py index 924de8d1..65356970 100644 --- a/neurom/features/__init__.py +++ b/neurom/features/__init__.py @@ -36,9 +36,11 @@ >>> ap_seg_len = features.get('segment_lengths', m, neurite_type=neurom.APICAL_DENDRITE) >>> ax_sec_len = features.get('section_lengths', m, neurite_type=neurom.AXON) """ + +import inspect import operator from enum import Enum -from functools import reduce +from functools import reduce, partial from neurom.core import Population, Morphology, Neurite from neurom.core.morphology import iter_neurites @@ -64,15 +66,30 @@ def _flatten_feature(feature_shape, feature_value): return reduce(operator.concat, feature_value, []) -def _get_neurites_feature_value(feature_, obj, neurite_filter, kwargs): +def _get_neurites_feature_value(feature_, obj, neurite_filter, use_subtrees, **kwargs): """Collects neurite feature values appropriately to feature's shape.""" kwargs.pop('neurite_type', None) # there is no 'neurite_type' arg in _NEURITE_FEATURES - return reduce(operator.add, - (feature_(n, **kwargs) for n in iter_neurites(obj, filt=neurite_filter)), - 0 if feature_.shape == () else []) + + return reduce( + operator.add, + ( + iter_neurites( + obj, + mapfun=partial(feature_, **kwargs), + filt=neurite_filter, + use_subtrees=use_subtrees, + ) + ), + 0 if feature_.shape == () else [] + ) + + +def _is_subtree_processing_applicable(feature_function): + """Returns true if feature's signature supports the use_subtrees kwarg.""" + return "use_subtrees" in inspect.signature(feature_function).parameters -def _get_feature_value_and_func(feature_name, obj, **kwargs): +def _get_feature_value_and_func(feature_name, obj, use_subtrees=False, **kwargs): """Obtain a feature from a set of morphology objects. Arguments: @@ -87,43 +104,72 @@ def _get_feature_value_and_func(feature_name, obj, **kwargs): # pylint: disable=too-many-branches is_obj_list = isinstance(obj, (list, tuple)) if not isinstance(obj, (Neurite, Morphology, Population)) and not is_obj_list: - raise NeuroMError('Only Neurite, Morphology, Population or list, tuple of Neurite,' - ' Morphology can be used for feature calculation') + raise NeuroMError( + "Only Neurite, Morphology, Population or list, tuple of Neurite, Morphology" + f"can be used for feature calculation. Got: {obj}" + ) neurite_filter = is_type(kwargs.get('neurite_type', NeuriteType.all)) res, feature_ = None, None if isinstance(obj, Neurite) or (is_obj_list and isinstance(obj[0], Neurite)): + # input is a neurite or a list of neurites if feature_name in _NEURITE_FEATURES: - assert 'neurite_type' not in kwargs, 'Cant apply "neurite_type" arg to a neurite with' \ - ' a neurite feature' + + assert 'neurite_type' not in kwargs, ( + 'Cant apply "neurite_type" arg to a neurite with a neurite feature' + ) + feature_ = _NEURITE_FEATURES[feature_name] + if isinstance(obj, Neurite): res = feature_(obj, **kwargs) else: res = [feature_(s, **kwargs) for s in obj] + elif isinstance(obj, Morphology): + # input is a morphology if feature_name in _MORPHOLOGY_FEATURES: + feature_ = _MORPHOLOGY_FEATURES[feature_name] + + if _is_subtree_processing_applicable(feature_): + kwargs["use_subtrees"] = use_subtrees + res = feature_(obj, **kwargs) + elif feature_name in _NEURITE_FEATURES: + feature_ = _NEURITE_FEATURES[feature_name] - res = _get_neurites_feature_value(feature_, obj, neurite_filter, kwargs) + res = _get_neurites_feature_value(feature_, obj, neurite_filter, use_subtrees, **kwargs) + elif isinstance(obj, Population) or (is_obj_list and isinstance(obj[0], Morphology)): # input is a morphology population or a list of morphs if feature_name in _POPULATION_FEATURES: feature_ = _POPULATION_FEATURES[feature_name] + + if _is_subtree_processing_applicable(feature_): + kwargs["use_subtrees"] = use_subtrees + res = feature_(obj, **kwargs) elif feature_name in _MORPHOLOGY_FEATURES: feature_ = _MORPHOLOGY_FEATURES[feature_name] + + if _is_subtree_processing_applicable(feature_): + kwargs["use_subtrees"] = use_subtrees + res = _flatten_feature(feature_.shape, [feature_(n, **kwargs) for n in obj]) elif feature_name in _NEURITE_FEATURES: feature_ = _NEURITE_FEATURES[feature_name] res = _flatten_feature( feature_.shape, - [_get_neurites_feature_value(feature_, n, neurite_filter, kwargs) for n in obj]) + [ + _get_neurites_feature_value(feature_, n, neurite_filter, use_subtrees, **kwargs) + for n in obj + ] + ) if res is None or feature_ is None: raise NeuroMError(f'Cant apply "{feature_name}" feature. Please check that it exists, ' @@ -132,7 +178,7 @@ def _get_feature_value_and_func(feature_name, obj, **kwargs): return res, feature_ -def get(feature_name, obj, **kwargs): +def get(feature_name, obj, use_subtrees=False, **kwargs): """Obtain a feature from a set of morphology objects. Features can be either Neurite, Morphology or Population features. For Neurite features see @@ -147,7 +193,7 @@ def get(feature_name, obj, **kwargs): Returns: List|Number: feature value as a list or a single number. """ - return _get_feature_value_and_func(feature_name, obj, **kwargs)[0] + return _get_feature_value_and_func(feature_name, obj, use_subtrees=use_subtrees, **kwargs)[0] def _register_feature(namespace: NameSpace, name, func, shape): diff --git a/neurom/features/bifurcation.py b/neurom/features/bifurcation.py index 423a3ef2..0bc25479 100644 --- a/neurom/features/bifurcation.py +++ b/neurom/features/bifurcation.py @@ -29,11 +29,12 @@ """Bifurcation point functions.""" import numpy as np + +import neurom.features.section from neurom import morphmath -from neurom.exceptions import NeuroMError from neurom.core.dataformat import COLS from neurom.core.morphology import Section -from neurom.features.section import section_mean_radius +from neurom.exceptions import NeuroMError def _raise_if_not_bifurcation(section): @@ -156,8 +157,8 @@ def sibling_ratio(bif_point, method='first'): n = bif_point.children[0].points[1, COLS.R] m = bif_point.children[1].points[1, COLS.R] if method == 'mean': - n = section_mean_radius(bif_point.children[0]) - m = section_mean_radius(bif_point.children[1]) + n = neurom.features.section.section_mean_radius(bif_point.children[0]) + m = neurom.features.section.section_mean_radius(bif_point.children[1]) return min(n, m) / max(n, m) @@ -182,7 +183,35 @@ def diameter_power_relation(bif_point, method='first'): d_child1 = bif_point.children[0].points[1, COLS.R] d_child2 = bif_point.children[1].points[1, COLS.R] if method == 'mean': - d_child = section_mean_radius(bif_point) - d_child1 = section_mean_radius(bif_point.children[0]) - d_child2 = section_mean_radius(bif_point.children[1]) + d_child = neurom.features.section.section_mean_radius(bif_point) + d_child1 = neurom.features.section.section_mean_radius(bif_point.children[0]) + d_child2 = neurom.features.section.section_mean_radius(bif_point.children[1]) return (d_child / d_child1)**(1.5) + (d_child / d_child2)**(1.5) + + +def downstream_pathlength_asymmetry( + bif_point, normalization_length=1.0, iterator_type=Section.ipreorder +): + """Calculates the downstream pathlength asymmetry at a bifurcation point. + + Args: + bif_point: Bifurcation section. + normalization_length: Constant to divide the result with. + iterator_type: Iterator type that specifies how the two subtrees are traversed. + + Returns: + The absolute difference between the downstream path distances of the two children, divided + by the normalization length. + """ + _raise_if_not_bifurcation(bif_point) + return ( + abs( + neurom.features.section.downstream_pathlength( + bif_point.children[0], iterator_type=iterator_type + ) + - neurom.features.section.downstream_pathlength( + bif_point.children[1], iterator_type=iterator_type + ), + ) + / normalization_length + ) diff --git a/neurom/features/morphology.py b/neurom/features/morphology.py index ddee3aa1..c36525c2 100644 --- a/neurom/features/morphology.py +++ b/neurom/features/morphology.py @@ -46,11 +46,14 @@ import warnings from functools import partial +from collections.abc import Iterable import math import numpy as np from neurom import morphmath -from neurom.core.morphology import iter_neurites, iter_segments, Morphology +from neurom.core.morphology import ( + iter_neurites, iter_sections, iter_segments, iter_points, Morphology +) from neurom.core.types import tree_type_checker as is_type from neurom.core.dataformat import COLS from neurom.core.types import NeuriteType @@ -63,9 +66,27 @@ feature = partial(feature, namespace=NameSpace.NEURON) -def _map_neurites(function, morph, neurite_type): +def _map_neurites(function, morph, neurite_type, use_subtrees=False): return list( - iter_neurites(morph, mapfun=function, filt=is_type(neurite_type)) + iter_neurites( + obj=morph, + mapfun=function, + filt=is_type(neurite_type), + use_subtrees=use_subtrees, + ) + ) + + +def _map_neurite_root_nodes(function, morph, neurite_type, use_subtrees=False): + neurites = iter_neurites(obj=morph, filt=is_type(neurite_type), use_subtrees=use_subtrees) + return [function(neurite.root_node) for neurite in neurites] + + +def _get_points(morph, neurite_type, use_subtrees=False): + return list( + iter_points(morph, section_filter=is_type(neurite_type)) + if use_subtrees + else iter_points(morph, neurite_filter=is_type(neurite_type)) ) @@ -92,35 +113,39 @@ def soma_radius(morph): @feature(shape=()) -def max_radial_distance(morph, neurite_type=NeuriteType.all): +def max_radial_distance(morph, origin=None, neurite_type=NeuriteType.all, use_subtrees=False): """Get the maximum radial distances of the termination sections.""" - term_radial_distances = _map_neurites(nf.max_radial_distance, morph, neurite_type) - - return max(term_radial_distances) if term_radial_distances else 0.0 + term_radial_distances = _map_neurites( + partial(nf.max_radial_distance, origin=origin), + morph, + neurite_type=neurite_type, + use_subtrees=use_subtrees, + ) + return max(term_radial_distances) if term_radial_distances else 0. @feature(shape=(...,)) -def number_of_sections_per_neurite(morph, neurite_type=NeuriteType.all): +def number_of_sections_per_neurite(morph, neurite_type=NeuriteType.all, use_subtrees=False): """List of numbers of sections per neurite.""" - return _map_neurites(nf.number_of_sections, morph, neurite_type) + return _map_neurites(nf.number_of_sections, morph, neurite_type, use_subtrees) @feature(shape=(...,)) -def total_length_per_neurite(morph, neurite_type=NeuriteType.all): +def total_length_per_neurite(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Neurite lengths.""" - return _map_neurites(nf.total_length, morph, neurite_type) + return _map_neurites(nf.total_length, morph, neurite_type, use_subtrees) @feature(shape=(...,)) -def total_area_per_neurite(morph, neurite_type=NeuriteType.all): +def total_area_per_neurite(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Neurite areas.""" - return _map_neurites(nf.total_area, morph, neurite_type) + return _map_neurites(nf.total_area, morph, neurite_type, use_subtrees) @feature(shape=(...,)) -def total_volume_per_neurite(morph, neurite_type=NeuriteType.all): +def total_volume_per_neurite(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Neurite volumes.""" - return _map_neurites(nf.total_volume, morph, neurite_type) + return _map_neurites(nf.total_volume, morph, neurite_type, use_subtrees) @feature(shape=(...,)) @@ -132,13 +157,12 @@ def trunk_origin_azimuths(morph, neurite_type=NeuriteType.all): The range of the azimuth angle [-pi, pi] radians """ - def azimuth(neurite): + def azimuth(root_node): """Azimuth of a neurite trunk.""" return morphmath.azimuth_from_vector( - morphmath.vector(neurite.root_node.points[0], morph.soma.center) + morphmath.vector(root_node.points[0], morph.soma.center) ) - - return _map_neurites(azimuth, morph, neurite_type) + return _map_neurite_root_nodes(azimuth, morph, neurite_type, use_subtrees=False) @feature(shape=(...,)) @@ -151,22 +175,22 @@ def trunk_origin_elevations(morph, neurite_type=NeuriteType.all): The range of the elevation angle [-pi/2, pi/2] radians """ - def elevation(neurite): + def elevation(root_node): """Elevation of a section.""" return morphmath.elevation_from_vector( - morphmath.vector(neurite.root_node.points[0], morph.soma.center) + morphmath.vector(root_node.points[0], morph.soma.center) ) - - return _map_neurites(elevation, morph, neurite_type) + return _map_neurite_root_nodes(elevation, morph, neurite_type, use_subtrees=False) @feature(shape=(...,)) -def trunk_vectors(morph, neurite_type=NeuriteType.all): +def trunk_vectors(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Calculate the vectors between all the trunks of the morphology and the soma center.""" - def vector_to_root_node(neurite): - return morphmath.vector(neurite.root_node.points[0], morph.soma.center) - - return _map_neurites(vector_to_root_node, morph, neurite_type) + def vector_from_soma_to_root(root_node): + return morphmath.vector(root_node.points[0], morph.soma.center) + return _map_neurite_root_nodes( + vector_from_soma_to_root, morph, neurite_type, use_subtrees=use_subtrees + ) @feature(shape=(...,)) @@ -247,6 +271,7 @@ def trunk_angles_inter_types( source_neurite_type=NeuriteType.apical_dendrite, target_neurite_type=NeuriteType.basal_dendrite, closest_component=None, + use_subtrees=False, ): """Calculate the angles between the trunks of the morph of a source type to target type. @@ -274,8 +299,12 @@ def trunk_angles_inter_types( If ``closest_component`` is not ``None``, only one of these values is returned for each couple. """ - source_vectors = trunk_vectors(morph, neurite_type=source_neurite_type) - target_vectors = trunk_vectors(morph, neurite_type=target_neurite_type) + source_vectors = trunk_vectors( + morph, neurite_type=source_neurite_type, use_subtrees=use_subtrees + ) + target_vectors = trunk_vectors( + morph, neurite_type=target_neurite_type, use_subtrees=use_subtrees + ) # In order to avoid the failure of the process in case the neurite_type does not exist if len(source_vectors) == 0 or len(target_vectors) == 0: @@ -310,6 +339,7 @@ def trunk_angles_from_vector( morph, neurite_type=NeuriteType.all, vector=None, + use_subtrees=False, ): """Calculate the angles between the trunks of the morph of a given type and a given vector. @@ -329,7 +359,7 @@ def trunk_angles_from_vector( if vector is None: vector = (0, 1, 0) - vectors = np.array(trunk_vectors(morph, neurite_type=neurite_type)) + vectors = np.array(trunk_vectors(morph, neurite_type=neurite_type, use_subtrees=use_subtrees)) # In order to avoid the failure of the process in case the neurite_type does not exist if len(vectors) == 0: @@ -357,6 +387,7 @@ def trunk_origin_radii( neurite_type=NeuriteType.all, min_length_filter=None, max_length_filter=None, + use_subtrees=False, ): """Radii of the trunk sections of neurites in a morph. @@ -379,10 +410,6 @@ def trunk_origin_radii( * else the mean radius of the points between the given ``min_length_filter`` and ``max_length_filter`` are returned. """ - if max_length_filter is None and min_length_filter is None: - return [n.root_node.points[0][COLS.R] - for n in iter_neurites(morph, filt=is_type(neurite_type))] - if min_length_filter is not None and min_length_filter <= 0: raise NeuroMError( "In 'trunk_origin_radii': the 'min_length_filter' value must be strictly greater " @@ -405,11 +432,17 @@ def trunk_origin_radii( "'max_length_filter' value." ) - def _mean_radius(neurite): - points = neurite.root_node.points + def trunk_first_radius(root_node): + return root_node.points[0][COLS.R] + + def trunk_mean_radius(root_node): + + points = root_node.points + interval_lengths = morphmath.interval_lengths(points) path_lengths = np.insert(np.cumsum(interval_lengths), 0, 0) valid_pts = np.ones(len(path_lengths), dtype=bool) + if min_length_filter is not None: valid_pts = (valid_pts & (path_lengths >= min_length_filter)) if not valid_pts.any(): @@ -419,6 +452,7 @@ def _mean_radius(neurite): "point is returned." ) return points[-1, COLS.R] + if max_length_filter is not None: valid_max = (path_lengths <= max_length_filter) valid_pts = (valid_pts & valid_max) @@ -430,34 +464,40 @@ def _mean_radius(neurite): ) # pylint: disable=invalid-unary-operand-type return points[~valid_max, COLS.R][0] + return points[valid_pts, COLS.R].mean() - return _map_neurites(_mean_radius, morph, neurite_type) + function = ( + trunk_first_radius + if max_length_filter is None and min_length_filter is None + else trunk_mean_radius + ) + + return _map_neurite_root_nodes(function, morph, neurite_type, use_subtrees=use_subtrees) @feature(shape=(...,)) -def trunk_section_lengths(morph, neurite_type=NeuriteType.all): +def trunk_section_lengths(morph, neurite_type=NeuriteType.all, use_subtrees=False): """List of lengths of trunk sections of neurites in a morph.""" - def trunk_section_length(neurite): - return morphmath.section_length(neurite.root_node.points) - - return _map_neurites(trunk_section_length, morph, neurite_type) + return _map_neurite_root_nodes(sf.section_length, morph, neurite_type, use_subtrees) @feature(shape=()) -def number_of_neurites(morph, neurite_type=NeuriteType.all): +def number_of_neurites(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Number of neurites in a morph.""" - return len(_map_neurites(lambda n: n, morph, neurite_type)) + return len(_map_neurite_root_nodes(lambda n: n, morph, neurite_type, use_subtrees)) @feature(shape=(...,)) -def neurite_volume_density(morph, neurite_type=NeuriteType.all): +def neurite_volume_density(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Get volume density per neurite.""" - return _map_neurites(nf.volume_density, morph, neurite_type) + return _map_neurites(nf.volume_density, morph, neurite_type, use_subtrees) @feature(shape=(...,)) -def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None): +def sholl_crossings( + morph, neurite_type=NeuriteType.all, center=None, radii=None, use_subtrees=False +): """Calculate crossings of neurites. Args: @@ -478,11 +518,11 @@ def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None center=morph.soma.center, radii=np.arange(0, 1000, 100)) """ - def _count_crossings(neurite, radius): + def count_crossings(section, radius): """Used to count_crossings of segments in neurite with radius.""" r2 = radius ** 2 count = 0 - for start, end in iter_segments(neurite): + for start, end in iter_segments(section): start_dist2, end_dist2 = (morphmath.point_dist2(center, start), morphmath.point_dist2(center, end)) @@ -499,13 +539,28 @@ def _count_crossings(neurite, radius): center = morph.soma.center if radii is None: radii = [morph.soma.radius] - return [sum(_count_crossings(neurite, r) - for neurite in iter_neurites(morph, filt=is_type(neurite_type))) - for r in radii] + + if isinstance(morph, Iterable): + sections = filter(is_type(neurite_type), morph) + else: + if use_subtrees: + sections = iter_sections(morph, section_filter=is_type(neurite_type)) + else: + sections = iter_sections(morph, neurite_filter=is_type(neurite_type)) + + counts_per_radius = [0 for _ in range(len(radii))] + + for section in sections: + for i, radius in enumerate(radii): + counts_per_radius[i] += count_crossings(section, radius) + + return counts_per_radius @feature(shape=(...,)) -def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None): +def sholl_frequency( + morph, neurite_type=NeuriteType.all, step_size=10, bins=None, use_subtrees=False +): """Perform Sholl frequency calculations on a morph. Args: @@ -514,6 +569,7 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None step_size(float): step size between Sholl radii bins(iterable of floats): custom binning to use for the Sholl radii. If None, it uses intervals of step_size between min and max radii of ``morphologies``. + use_subtrees: Enable mixed subtree processing Note: Given a morphology, the soma center is used for the concentric circles, @@ -525,62 +581,61 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None If a `neurite_type` is specified and there are no trees corresponding to it, an empty list will be returned. """ - neurite_filter = is_type(neurite_type) - if bins is None: min_soma_edge = morph.soma.radius - max_radius_per_neurite = [ - np.max(np.linalg.norm(n.points[:, COLS.XYZ] - morph.soma.center, axis=1)) - for n in morph.neurites if neurite_filter(n) + if use_subtrees: + sections = iter_sections(morph, section_filter=is_type(neurite_type)) + else: + sections = iter_sections(morph, neurite_filter=is_type(neurite_type)) + + max_radius_per_section = [ + np.max(np.linalg.norm(section.points[:, COLS.XYZ] - morph.soma.center, axis=1)) + for section in sections ] - if not max_radius_per_neurite: + if not max_radius_per_section: return [] - bins = np.arange(min_soma_edge, min_soma_edge + max(max_radius_per_neurite), step_size) + bins = np.arange(min_soma_edge, min_soma_edge + max(max_radius_per_section), step_size) - return sholl_crossings(morph, neurite_type, morph.soma.center, bins) + return sholl_crossings(morph, neurite_type, morph.soma.center, bins, use_subtrees=use_subtrees) -def _extent_along_axis(morph, axis, neurite_type): +def _extent_along_axis(morph, axis, neurite_type, use_subtrees=False): """Returns the total extent of the morpholog neurites. The morphology is filtered by neurite type and the extent is calculated along the coordinate axis direction (e.g. COLS.X). """ - it_points = ( - p - for n in iter_neurites(morph, filt=is_type(neurite_type)) - for p in n.points[:, axis] - ) - try: - return abs(np.ptp(np.fromiter(it_points, dtype=np.float32))) - except ValueError: - # a ValueError is thrown when there are no points passed to ptp + points = _get_points(morph, neurite_type, use_subtrees=use_subtrees) + + if not points: return 0.0 + return abs(np.ptp(np.asarray(points)[:, axis])) + @feature(shape=()) -def total_width(morph, neurite_type=NeuriteType.all): +def total_width(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Extent of morphology along axis x.""" - return _extent_along_axis(morph, axis=COLS.X, neurite_type=neurite_type) + return _extent_along_axis(morph, COLS.X, neurite_type, use_subtrees) @feature(shape=()) -def total_height(morph, neurite_type=NeuriteType.all): +def total_height(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Extent of morphology along axis y.""" - return _extent_along_axis(morph, axis=COLS.Y, neurite_type=neurite_type) + return _extent_along_axis(morph, COLS.Y, neurite_type, use_subtrees) @feature(shape=()) -def total_depth(morph, neurite_type=NeuriteType.all): +def total_depth(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Extent of morphology along axis z.""" - return _extent_along_axis(morph, axis=COLS.Z, neurite_type=neurite_type) + return _extent_along_axis(morph, COLS.Z, neurite_type, use_subtrees) @feature(shape=()) -def volume_density(morph, neurite_type=NeuriteType.all): +def volume_density(morph, neurite_type=NeuriteType.all, use_subtrees=False): """Get the volume density. The volume density is defined as the ratio of the neurite volume and @@ -589,12 +644,7 @@ def volume_density(morph, neurite_type=NeuriteType.all): .. note:: Returns `np.nan` if the convex hull computation fails or there are not points available due to neurite type filtering. """ - # note: duplicate points are present but do not affect convex hull calculation - points = [ - point - for point_list in iter_neurites(morph, mapfun=sf.section_points, filt=is_type(neurite_type)) - for point in point_list - ] + points = _get_points(morph, neurite_type, use_subtrees) if not points: return np.nan @@ -604,12 +654,14 @@ def volume_density(morph, neurite_type=NeuriteType.all): if morph_hull is None: return np.nan - total_volume = sum(iter_neurites(morph, mapfun=nf.total_volume, filt=is_type(neurite_type))) + total_volume = sum(total_volume_per_neurite( + morph, neurite_type=neurite_type, use_subtrees=use_subtrees) + ) return total_volume / morph_hull.volume -def _unique_projected_points(morph, projection_plane, neurite_type): +def _unique_projected_points(morph, projection_plane, neurite_type, use_subtrees=False): key = "".join(sorted(projection_plane.lower())) @@ -623,9 +675,7 @@ def _unique_projected_points(morph, projection_plane, neurite_type): f"Please select 'xy', 'xz', or 'yz'." ) from e - points = list( - iter_neurites(morph, mapfun=sf.section_points, filt=is_type(neurite_type)) - ) + points = _get_points(morph, neurite_type, use_subtrees) if len(points) == 0: return np.empty(shape=(0, 3), dtype=np.float32) @@ -634,23 +684,24 @@ def _unique_projected_points(morph, projection_plane, neurite_type): @feature(shape=()) -def aspect_ratio(morph, neurite_type=NeuriteType.all, projection_plane="xy"): +def aspect_ratio(morph, neurite_type=NeuriteType.all, projection_plane="xy", use_subtrees=False): """Calculates the min/max ratio of the principal direction extents along the plane. Args: morph: Morphology object. neurite_type: The neurite type to use. By default all neurite types are used. projection_plane: Projection plane to use for the calculation. One of ('xy', 'xz', 'yz'). + use_subtrees: Enable mixed subtree processing Returns: The aspect ratio feature of the morphology points. """ - projected_points = _unique_projected_points(morph, projection_plane, neurite_type) + projected_points = _unique_projected_points(morph, projection_plane, neurite_type, use_subtrees) return np.nan if len(projected_points) == 0 else morphmath.aspect_ratio(projected_points) @feature(shape=()) -def circularity(morph, neurite_type=NeuriteType.all, projection_plane="xy"): +def circularity(morph, neurite_type=NeuriteType.all, projection_plane="xy", use_subtrees=False): """Calculates the circularity of the morphology points along the plane. The circularity is defined as the 4 * pi * area of the convex hull over its @@ -661,16 +712,17 @@ def circularity(morph, neurite_type=NeuriteType.all, projection_plane="xy"): neurite_type: The neurite type to use. By default all neurite types are used. projection_plane: Projection plane to use for the calculation. One of ('xy', 'xz', 'yz'). + use_subtrees: Enable mixed subtree processing Returns: The circularity of the morphology points. """ - projected_points = _unique_projected_points(morph, projection_plane, neurite_type) + projected_points = _unique_projected_points(morph, projection_plane, neurite_type, use_subtrees) return np.nan if len(projected_points) == 0 else morphmath.circularity(projected_points) @feature(shape=()) -def shape_factor(morph, neurite_type=NeuriteType.all, projection_plane="xy"): +def shape_factor(morph, neurite_type=NeuriteType.all, projection_plane="xy", use_subtrees=False): """Calculates the shape factor of the morphology points along the plane. The shape factor is defined as the ratio of the convex hull area over max squared @@ -681,16 +733,17 @@ def shape_factor(morph, neurite_type=NeuriteType.all, projection_plane="xy"): neurite_type: The neurite type to use. By default all neurite types are used. projection_plane: Projection plane to use for the calculation. One of ('xy', 'xz', 'yz'). + use_subtrees: Enable mixed subtree processing Returns: The shape factor of the morphology points. """ - projected_points = _unique_projected_points(morph, projection_plane, neurite_type) + projected_points = _unique_projected_points(morph, projection_plane, neurite_type, use_subtrees) return np.nan if len(projected_points) == 0 else morphmath.shape_factor(projected_points) @feature(shape=()) -def length_fraction_above_soma(morph, neurite_type=NeuriteType.all, up="Y"): +def length_fraction_above_soma(morph, neurite_type=NeuriteType.all, up="Y", use_subtrees=False): """Returns the length fraction of the segments that have their midpoints higher than the soma. Args: @@ -707,7 +760,11 @@ def length_fraction_above_soma(morph, neurite_type=NeuriteType.all, up="Y"): raise NeuroMError(f"Unknown axis {axis}. Please choose 'X', 'Y', or 'Z'.") col = getattr(COLS, axis) - segments = list(iter_segments(morph, neurite_filter=is_type(neurite_type))) + + if use_subtrees: + segments = list(iter_segments(morph, neurite_filter=is_type(neurite_type))) + else: + segments = list(iter_segments(morph, section_filter=is_type(neurite_type))) if not segments: return np.nan diff --git a/neurom/features/neurite.py b/neurom/features/neurite.py index 6aac5e9f..7c972779 100644 --- a/neurom/features/neurite.py +++ b/neurom/features/neurite.py @@ -48,124 +48,151 @@ import numpy as np from neurom import morphmath -from neurom.core.morphology import Section +from neurom import utils +from neurom.core.types import NeuriteType +from neurom.core.morphology import Section, iter_points from neurom.core.dataformat import COLS from neurom.features import NameSpace, feature, bifurcation as bf, section as sf from neurom.morphmath import convex_hull +from neurom.core.types import tree_type_checker as is_type + feature = partial(feature, namespace=NameSpace.NEURITE) L = logging.getLogger(__name__) -def _map_sections(fun, neurite, iterator_type=Section.ipreorder): +def _map_sections(fun, neurite, iterator_type=Section.ipreorder, section_type=NeuriteType.all): """Map `fun` to all the sections.""" - return list(map(fun, iterator_type(neurite.root_node))) + check_type = is_type(section_type) + + def homogeneous_filter(section): + return check_type(section) and Section.is_homogeneous_point(section) + + # forking sections cannot be heterogeneous + if ( + iterator_type in {Section.ibifurcation_point, Section.iforking_point} + and section_type != NeuriteType.all + ): + filt = homogeneous_filter + else: + filt = check_type + + return list(map(fun, filter(filt, iterator_type(neurite.root_node)))) @feature(shape=()) -def max_radial_distance(neurite): +def max_radial_distance(neurite, origin=None, section_type=NeuriteType.all): """Get the maximum radial distances of the termination sections.""" - term_radial_distances = section_term_radial_distances(neurite) + term_radial_distances = section_term_radial_distances( + neurite, origin=origin, section_type=section_type + ) return max(term_radial_distances) if term_radial_distances else 0. @feature(shape=()) -def number_of_segments(neurite): +def number_of_segments(neurite, section_type=NeuriteType.all): """Number of segments.""" - return sum(_map_sections(sf.number_of_segments, neurite)) + return sum(_map_sections(sf.number_of_segments, neurite, section_type=section_type)) @feature(shape=()) -def number_of_sections(neurite, iterator_type=Section.ipreorder): +def number_of_sections(neurite, iterator_type=Section.ipreorder, section_type=NeuriteType.all): """Number of sections. For a morphology it will be a sum of all neurites sections numbers.""" - return len(_map_sections(lambda s: s, neurite, iterator_type=iterator_type)) + return len( + _map_sections(lambda s: s, neurite, iterator_type=iterator_type, section_type=section_type) + ) @feature(shape=()) -def number_of_bifurcations(neurite): +def number_of_bifurcations(neurite, section_type=NeuriteType.all): """Number of bf points.""" - return number_of_sections(neurite, iterator_type=Section.ibifurcation_point) + return number_of_sections( + neurite, iterator_type=Section.ibifurcation_point, section_type=section_type + ) @feature(shape=()) -def number_of_forking_points(neurite): +def number_of_forking_points(neurite, section_type=NeuriteType.all): """Number of forking points.""" - return number_of_sections(neurite, iterator_type=Section.iforking_point) + return number_of_sections( + neurite, iterator_type=Section.iforking_point, section_type=section_type + ) @feature(shape=()) -def number_of_leaves(neurite): +def number_of_leaves(neurite, section_type=NeuriteType.all): """Number of leaves points.""" - return number_of_sections(neurite, iterator_type=Section.ileaf) + return number_of_sections(neurite, iterator_type=Section.ileaf, section_type=section_type) @feature(shape=()) -def total_length(neurite): +def total_length(neurite, section_type=NeuriteType.all): """Neurite length. For a morphology it will be a sum of all neurite lengths.""" - return sum(_map_sections(sf.section_length, neurite)) + return sum(_map_sections(sf.section_length, neurite, section_type=section_type)) @feature(shape=()) -def total_area(neurite): +def total_area(neurite, section_type=NeuriteType.all): """Neurite surface area. For a morphology it will be a sum of all neurite areas. The area is defined as the sum of the area of the sections. """ - return sum(_map_sections(sf.section_area, neurite)) + return sum(_map_sections(sf.section_area, neurite, section_type=section_type)) @feature(shape=()) -def total_volume(neurite): +def total_volume(neurite, section_type=NeuriteType.all): """Neurite volume. For a morphology it will be a sum of neurites volumes.""" - return sum(_map_sections(sf.section_volume, neurite)) + return sum(_map_sections(sf.section_volume, neurite, section_type=section_type)) @feature(shape=(...,)) -def section_lengths(neurite): +def section_lengths(neurite, section_type=NeuriteType.all): """Section lengths.""" - return _map_sections(sf.section_length, neurite) + return _map_sections(sf.section_length, neurite, section_type=section_type) @feature(shape=(...,)) -def section_term_lengths(neurite): +def section_term_lengths(neurite, section_type=NeuriteType.all): """Termination section lengths.""" - return _map_sections(sf.section_length, neurite, Section.ileaf) + return _map_sections(sf.section_length, neurite, Section.ileaf, section_type) @feature(shape=(...,)) -def section_bif_lengths(neurite): +def section_bif_lengths(neurite, section_type=NeuriteType.all): """Bifurcation section lengths.""" - return _map_sections(sf.section_length, neurite, Section.ibifurcation_point) + return _map_sections(sf.section_length, neurite, Section.ibifurcation_point, section_type) @feature(shape=(...,)) -def section_branch_orders(neurite): +def section_branch_orders(neurite, section_type=NeuriteType.all): """Section branch orders.""" - return _map_sections(sf.branch_order, neurite) + return _map_sections(sf.branch_order, neurite, section_type=section_type) @feature(shape=(...,)) -def section_bif_branch_orders(neurite): +def section_bif_branch_orders(neurite, section_type=NeuriteType.all): """Bifurcation section branch orders.""" - return _map_sections(sf.branch_order, neurite, Section.ibifurcation_point) + return _map_sections( + sf.branch_order, neurite, Section.ibifurcation_point, section_type=section_type + ) @feature(shape=(...,)) -def section_term_branch_orders(neurite): +def section_term_branch_orders(neurite, section_type=NeuriteType.all): """Termination section branch orders.""" - return _map_sections(sf.branch_order, neurite, Section.ileaf) + return _map_sections(sf.branch_order, neurite, Section.ileaf, section_type=section_type) @feature(shape=(...,)) -def section_path_distances(neurite): +def section_path_distances(neurite, iterator_type=Section.ipreorder, section_type=NeuriteType.all): """Path lengths.""" - - def pl2(node): - """Calculate the path length using cached section lengths.""" - return sum(n.length for n in node.iupstream()) - - return _map_sections(pl2, neurite) + return _map_sections( + partial(sf.section_path_length, stop_node=neurite.root_node), + neurite, + iterator_type=iterator_type, section_type=section_type + ) ################################################################################ @@ -173,120 +200,125 @@ def pl2(node): ################################################################################ -def _map_segments(func, neurite): +def _map_segments(func, neurite, section_type=NeuriteType.all): """Map `func` to all the segments. `func` accepts a section and returns list of values corresponding to each segment. """ - return [ - segment_value - for section in Section.ipreorder(neurite.root_node) - for segment_value in func(section) - ] + return list(utils.flatten(_map_sections(func, neurite, section_type=section_type))) @feature(shape=(...,)) -def segment_lengths(neurite): +def segment_lengths(neurite, section_type=NeuriteType.all): """Lengths of the segments.""" - return _map_segments(sf.segment_lengths, neurite) + return _map_segments(sf.segment_lengths, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_areas(neurite): +def segment_areas(neurite, section_type=NeuriteType.all): """Areas of the segments.""" - return _map_segments(sf.segment_areas, neurite) + return _map_segments(sf.segment_areas, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_volumes(neurite): +def segment_volumes(neurite, section_type=NeuriteType.all): """Volumes of the segments.""" - return _map_segments(sf.segment_volumes, neurite) + return _map_segments(sf.segment_volumes, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_radii(neurite): +def segment_radii(neurite, section_type=NeuriteType.all): """Arithmetic mean of the radii of the points in segments.""" - return _map_segments(sf.segment_mean_radii, neurite) + return _map_segments(sf.segment_mean_radii, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_taper_rates(neurite): +def segment_taper_rates(neurite, section_type=NeuriteType.all): """Diameters taper rates of the segments. The taper rate is defined as the absolute radii differences divided by length of the section """ - return _map_segments(sf.segment_taper_rates, neurite) + return _map_segments(sf.segment_taper_rates, neurite, section_type=section_type) @feature(shape=(...,)) -def section_taper_rates(neurite): +def section_taper_rates(neurite, section_type=NeuriteType.all): """Diameter taper rates of the sections from root to tip. Taper rate is defined here as the linear fit along a section. It is expected to be negative for morphologies. """ - return _map_sections(sf.taper_rate, neurite) + return _map_sections(sf.taper_rate, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_meander_angles(neurite): +def segment_meander_angles(neurite, section_type=NeuriteType.all): """Inter-segment opening angles in a section.""" - return _map_segments(sf.section_meander_angles, neurite) + return _map_segments(sf.section_meander_angles, neurite, section_type=section_type) @feature(shape=(..., 3)) -def segment_midpoints(neurite): +def segment_midpoints(neurite, section_type=NeuriteType.all): """Return a list of segment mid-points.""" - return _map_segments(sf.segment_midpoints, neurite) + return _map_segments(sf.segment_midpoints, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_path_lengths(neurite): +def segment_path_lengths(neurite, section_type=NeuriteType.all): """Returns pathlengths between all non-root points and their root point.""" pathlength = {} - def segments_pathlength(section): + def segments_path_length(section): if section.id not in pathlength: - if section.parent: - pathlength[section.id] = section.parent.length + pathlength[section.parent.id] - else: - pathlength[section.id] = 0 + + pathlength[section.id] = ( + 0.0 + if section.id == neurite.root_node.id + else section.parent.length + pathlength[section.parent.id] + ) + return pathlength[section.id] + np.cumsum(sf.segment_lengths(section)) - return _map_segments(segments_pathlength, neurite) + return _map_segments(segments_path_length, neurite, section_type=section_type) @feature(shape=(...,)) -def segment_radial_distances(neurite, origin=None): +def segment_radial_distances(neurite, origin=None, section_type=NeuriteType.all): """Returns the list of distances between all segment mid points and origin.""" - pos = neurite.root_node.points[0] if origin is None else origin - - def radial_distances(section): - """List of distances between the mid point of each segment and pos.""" - mid_pts = 0.5 * (section.points[:-1, COLS.XYZ] + section.points[1:, COLS.XYZ]) - return np.linalg.norm(mid_pts - pos[COLS.XYZ], axis=1) - - return _map_segments(radial_distances, neurite) + origin = neurite.root_node.points[0, COLS.XYZ] if origin is None else origin + return _map_segments( + func=partial(sf.segment_midpoint_radial_distances, origin=origin), + neurite=neurite, + section_type=section_type + ) @feature(shape=(...,)) -def local_bifurcation_angles(neurite): +def local_bifurcation_angles(neurite, section_type=NeuriteType.all): """Get a list of local bf angles.""" - return _map_sections(bf.local_bifurcation_angle, - neurite, - iterator_type=Section.ibifurcation_point) + return _map_sections( + bf.local_bifurcation_angle, + neurite, + iterator_type=Section.ibifurcation_point, + section_type=section_type, + ) @feature(shape=(...,)) -def remote_bifurcation_angles(neurite): +def remote_bifurcation_angles(neurite, section_type=NeuriteType.all): """Get a list of remote bf angles.""" - return _map_sections(bf.remote_bifurcation_angle, - neurite, - iterator_type=Section.ibifurcation_point) + return _map_sections( + bf.remote_bifurcation_angle, + neurite, + iterator_type=Section.ibifurcation_point, + section_type=section_type, + ) @feature(shape=(...,)) -def partition_asymmetry(neurite, variant='branch-order', method='petilla'): +def partition_asymmetry( + neurite, variant='branch-order', method='petilla', section_type=NeuriteType.all +): """Partition asymmetry at bf points. Variant: length is a different definition, as the absolute difference in @@ -295,46 +327,58 @@ def partition_asymmetry(neurite, variant='branch-order', method='petilla'): :func:`neurom.features.bifurcationfunc.partition_asymmetry` """ if variant not in {'branch-order', 'length'}: - raise ValueError('Please provide a valid variant for partition asymmetry,' - f'found {variant}') + raise ValueError( + "Please provide a valid variant for partition asymmetry. " + f"Expected 'branch-order' or 'length', got {variant}." + ) if method not in {'petilla', 'uylings'}: - raise ValueError('Please provide a valid method for partition asymmetry,' - 'either "petilla" or "uylings"') + raise ValueError( + "Please provide a valid method for partition asymmetry. " + f"Expected 'petilla' or 'uylings', got {method}." + ) + + # create a downstream iterator that is filtered by the section type + it_type = utils.filtered_iterator(is_type(section_type), Section.ipreorder) if variant == 'branch-order': return _map_sections( - partial(bf.partition_asymmetry, uylings=method == 'uylings'), + partial(bf.partition_asymmetry, uylings=method == 'uylings', iterator_type=it_type), neurite, - Section.ibifurcation_point) + iterator_type=Section.ibifurcation_point, + section_type=section_type + ) - asymmetries = [] - neurite_length = total_length(neurite) - for section in Section.ibifurcation_point(neurite.root_node): - pathlength_diff = abs(sf.downstream_pathlength(section.children[0]) - - sf.downstream_pathlength(section.children[1])) - asymmetries.append(pathlength_diff / neurite_length) - return asymmetries + return _map_sections( + partial( + bf.downstream_pathlength_asymmetry, + normalization_length=total_length(neurite, section_type=section_type), + iterator_type=it_type, + ), + neurite, + iterator_type=Section.ibifurcation_point, + section_type=section_type + ) @feature(shape=(...,)) -def partition_asymmetry_length(neurite, method='petilla'): +def partition_asymmetry_length(neurite, method='petilla', section_type=NeuriteType.all): """'partition_asymmetry' feature with `variant='length'`. Because it is often used, it has a dedicated feature. """ - return partition_asymmetry(neurite, 'length', method) + return partition_asymmetry(neurite, 'length', method, section_type=section_type) @feature(shape=(...,)) -def bifurcation_partitions(neurite): +def bifurcation_partitions(neurite, section_type=NeuriteType.all): """Partition at bf points.""" - return _map_sections(bf.bifurcation_partition, - neurite, - Section.ibifurcation_point) + return _map_sections( + bf.bifurcation_partition, neurite, Section.ibifurcation_point, section_type=section_type + ) @feature(shape=(...,)) -def sibling_ratios(neurite, method='first'): +def sibling_ratios(neurite, method='first', section_type=NeuriteType.all): """Sibling ratios at bf points. The sibling ratio is the ratio between the diameters of the @@ -342,25 +386,28 @@ def sibling_ratios(neurite, method='first'): 0 and 1. Method argument allows one to consider mean diameters along the child section instead of diameter of the first point. """ - return _map_sections(partial(bf.sibling_ratio, method=method), - neurite, - Section.ibifurcation_point) + return _map_sections( + partial(bf.sibling_ratio, method=method), + neurite, + Section.ibifurcation_point, + section_type=section_type, + ) @feature(shape=(..., 2)) -def partition_pairs(neurite): +def partition_pairs(neurite, section_type=NeuriteType.all): """Partition pairs at bf points. Partition pair is defined as the number of bifurcations at the two daughters of the bifurcating section """ - return _map_sections(bf.partition_pair, - neurite, - Section.ibifurcation_point) + return _map_sections( + bf.partition_pair, neurite, Section.ibifurcation_point, section_type=section_type + ) @feature(shape=(...,)) -def diameter_power_relations(neurite, method='first'): +def diameter_power_relations(neurite, method='first', section_type=NeuriteType.all): """Calculate the diameter power relation at a bf point. Diameter power relation is defined in https://www.ncbi.nlm.nih.gov/pubmed/18568015 @@ -368,13 +415,18 @@ def diameter_power_relations(neurite, method='first'): This quantity gives an indication of how far the branching is from the Rall ratio (when =1). """ - return _map_sections(partial(bf.diameter_power_relation, method=method), - neurite, - Section.ibifurcation_point) + return _map_sections( + partial(bf.diameter_power_relation, method=method), + neurite, + Section.ibifurcation_point, + section_type=section_type, + ) @feature(shape=(...,)) -def section_radial_distances(neurite, origin=None, iterator_type=Section.ipreorder): +def section_radial_distances( + neurite, origin=None, iterator_type=Section.ipreorder, section_type=NeuriteType.all +): """Section radial distances. The iterator_type can be used to select only terminal sections (ileaf) @@ -383,29 +435,32 @@ def section_radial_distances(neurite, origin=None, iterator_type=Section.ipreord pos = neurite.root_node.points[0] if origin is None else origin return _map_sections(partial(sf.section_radial_distance, origin=pos), neurite, - iterator_type) + iterator_type, + section_type=section_type) @feature(shape=(...,)) -def section_term_radial_distances(neurite, origin=None): +def section_term_radial_distances(neurite, origin=None, section_type=NeuriteType.all): """Get the radial distances of the termination sections.""" - return section_radial_distances(neurite, origin, Section.ileaf) + return section_radial_distances(neurite, origin, Section.ileaf, section_type=section_type) @feature(shape=(...,)) -def section_bif_radial_distances(neurite, origin=None): +def section_bif_radial_distances(neurite, origin=None, section_type=NeuriteType.all): """Get the radial distances of the bf sections.""" - return section_radial_distances(neurite, origin, Section.ibifurcation_point) + return section_radial_distances( + neurite, origin, Section.ibifurcation_point, section_type=section_type + ) @feature(shape=(...,)) -def terminal_path_lengths(neurite): +def terminal_path_lengths(neurite, section_type=NeuriteType.all): """Get the path lengths to each terminal point.""" - return _map_sections(sf.section_path_length, neurite, Section.ileaf) + return section_path_distances(neurite, iterator_type=Section.ileaf, section_type=section_type) @feature(shape=()) -def volume_density(neurite): +def volume_density(neurite, section_type=NeuriteType.all): """Get the volume density. The volume density is defined as the ratio of the neurite volume and @@ -416,46 +471,59 @@ def volume_density(neurite): .. note:: Returns `np.nan` if the convex hull computation fails. """ - neurite_hull = convex_hull(neurite.points[:, COLS.XYZ]) - return neurite.volume / neurite_hull.volume if neurite_hull is not None else np.nan + neurite_volume = total_volume(neurite, section_type=section_type) + + def get_points(section): + return section.points[:, COLS.XYZ].tolist() + + # note: duplicate points included but not affect the convex hull calculation + points = list( + utils.flatten(_map_sections(get_points, neurite, section_type=section_type)) + ) + + hull = convex_hull(points) + + return neurite_volume / hull.volume if hull is not None else np.nan @feature(shape=(...,)) -def section_volumes(neurite): +def section_volumes(neurite, section_type=NeuriteType.all): """Section volumes.""" - return _map_sections(sf.section_volume, neurite) + return _map_sections(sf.section_volume, neurite, section_type=section_type) @feature(shape=(...,)) -def section_areas(neurite): +def section_areas(neurite, section_type=NeuriteType.all): """Section areas.""" - return _map_sections(sf.section_area, neurite) + return _map_sections(sf.section_area, neurite, section_type=section_type) @feature(shape=(...,)) -def section_tortuosity(neurite): +def section_tortuosity(neurite, section_type=NeuriteType.all): """Section tortuosities.""" - return _map_sections(sf.section_tortuosity, neurite) + return _map_sections(sf.section_tortuosity, neurite, section_type=section_type) @feature(shape=(...,)) -def section_end_distances(neurite): +def section_end_distances(neurite, section_type=NeuriteType.all): """Section end to end distances.""" - return _map_sections(sf.section_end_distance, neurite) + return _map_sections(sf.section_end_distance, neurite, section_type=section_type) @feature(shape=(...,)) -def principal_direction_extents(neurite, direction=0): +def principal_direction_extents(neurite, direction=0, section_type=NeuriteType.all): """Principal direction extent of neurites in morphologies. Note: Principal direction extents are always sorted in descending order. Therefore, by default the maximal principal direction extent is returned. """ - return [morphmath.principal_direction_extent(neurite.points[:, COLS.XYZ])[direction]] + points = list(iter_points(neurite, section_filter=is_type(section_type))) + + return [morphmath.principal_direction_extent(np.unique(points, axis=0))[direction]] @feature(shape=(...,)) -def section_strahler_orders(neurite): +def section_strahler_orders(neurite, section_type=NeuriteType.all): """Inter-segment opening angles in a section.""" - return _map_sections(sf.strahler_order, neurite) + return _map_sections(sf.strahler_order, neurite, section_type=section_type) diff --git a/neurom/features/population.py b/neurom/features/population.py index 98cbd4de..4720e3eb 100644 --- a/neurom/features/population.py +++ b/neurom/features/population.py @@ -45,15 +45,18 @@ from neurom.core.dataformat import COLS from neurom.core.types import NeuriteType +from neurom.core.morphology import iter_sections from neurom.core.types import tree_type_checker as is_type from neurom.features import feature, NameSpace -from neurom.features.morphology import sholl_crossings +from neurom.features import morphology as mf feature = partial(feature, namespace=NameSpace.POPULATION) @feature(shape=(...,)) -def sholl_frequency(morphs, neurite_type=NeuriteType.all, step_size=10, bins=None): +def sholl_frequency( + morphs, neurite_type=NeuriteType.all, step_size=10, bins=None, use_subtrees=False +): """Perform Sholl frequency calculations on a population of morphs. Args: @@ -62,6 +65,7 @@ def sholl_frequency(morphs, neurite_type=NeuriteType.all, step_size=10, bins=Non step_size(float): step size between Sholl radii bins(iterable of floats): custom binning to use for the Sholl radii. If None, it uses intervals of step_size between min and max radii of ``morphs``. + use_subtrees (bool): Enable mixed subtree processing. Note: Given a population, the concentric circles range from the smallest soma radius to the @@ -72,13 +76,27 @@ def sholl_frequency(morphs, neurite_type=NeuriteType.all, step_size=10, bins=Non neurite_filter = is_type(neurite_type) if bins is None: + + section_iterator = ( + partial(iter_sections, section_filter=neurite_filter) + if use_subtrees + else partial(iter_sections, neurite_filter=neurite_filter) + ) + + max_radius_per_section = [ + np.max(np.linalg.norm(section.points[:, COLS.XYZ] - morph.soma.center, axis=1)) + for morph in morphs + for section in section_iterator(morph) + ] + + if not max_radius_per_section: + return [] + min_soma_edge = min(n.soma.radius for n in morphs) - max_radii = max(np.max(np.linalg.norm(n.points[:, COLS.XYZ], axis=1)) - for m in morphs - for n in m.neurites if neurite_filter(n)) - bins = np.arange(min_soma_edge, min_soma_edge + max_radii, step_size) + + bins = np.arange(min_soma_edge, min_soma_edge + max(max_radius_per_section), step_size) return np.array([ - sholl_crossings(m, neurite_type, m.soma.center, bins) + mf.sholl_crossings(m, neurite_type, m.soma.center, bins, use_subtrees=use_subtrees) for m in morphs - ]).sum(axis=0) + ]).sum(axis=0).tolist() diff --git a/neurom/features/section.py b/neurom/features/section.py index 259699f3..e9c5ed6d 100644 --- a/neurom/features/section.py +++ b/neurom/features/section.py @@ -33,6 +33,7 @@ from neurom import morphmath as mm from neurom.core.dataformat import COLS from neurom.core.morphology import iter_segments +from neurom.core.morphology import Section from neurom.morphmath import interval_lengths @@ -41,9 +42,14 @@ def section_points(section): return section.points[:, COLS.XYZ] -def section_path_length(section): - """Path length from section to root.""" - return sum(s.length for s in section.iupstream()) +def section_path_length(section, stop_node=None): + """Path length from section to root. + + Args: + section: Section object. + stop_node: Node to stop the upstream traversal. If None, it stops when no parent is found. + """ + return sum(map(section_length, section.iupstream(stop_node=stop_node))) def section_length(section): @@ -137,6 +143,13 @@ def segment_midpoints(section): return np.divide(np.add(pts[:-1], pts[1:]), 2.0).tolist() +def segment_midpoint_radial_distances(section, origin=None): + """Returns the list of segment midpoint radial distances to the origin.""" + origin = np.zeros(3, dtype=float) if origin is None else origin + midpoints = np.array(segment_midpoints(section)) + return np.linalg.norm(midpoints - origin, axis=1).tolist() + + def segment_taper_rates(section): """Returns the list of segment taper rates within the section.""" pts = section.points[:, COLS.XYZR] @@ -213,6 +226,6 @@ def section_mean_radius(section): return np.sum(mean_radii * lengths) / np.sum(lengths) -def downstream_pathlength(section): +def downstream_pathlength(section, iterator_type=Section.ipreorder): """Compute the total downstream length starting from a section.""" - return sum(sec.length for sec in section.ipreorder()) + return sum(sec.length for sec in iterator_type(section)) diff --git a/neurom/utils.py b/neurom/utils.py index 90ab2a4b..a2a04eab 100644 --- a/neurom/utils.py +++ b/neurom/utils.py @@ -136,3 +136,11 @@ def str_to_plane(plane): def flatten(list_of_lists): """Flatten one level of nesting.""" return chain.from_iterable(list_of_lists) + + +def filtered_iterator(predicate, iterator_type): + """Returns an iterator function that is filtered by the predicate.""" + @wraps(iterator_type) + def composed(*args, **kwargs): + return filter(predicate, iterator_type(*args, **kwargs)) + return composed diff --git a/tests/core/test_section.py b/tests/core/test_section.py index 25fc4817..93708504 100644 --- a/tests/core/test_section.py +++ b/tests/core/test_section.py @@ -45,6 +45,8 @@ def test_section_base_func(): assert_almost_equal(section.area, 31.41592653589793) assert_almost_equal(section.volume, 15.707963267948964) + # __nonzero__ + assert section def test_section_tree(): m = nm.load_morphology(str(SWC_PATH / 'simple.swc')) diff --git a/tests/data/swc/heterogeneous_morphology.swc b/tests/data/swc/heterogeneous_morphology.swc new file mode 100644 index 00000000..d3b26ba5 --- /dev/null +++ b/tests/data/swc/heterogeneous_morphology.swc @@ -0,0 +1,25 @@ +# Created by MorphIO v3.3.3 +# index type X Y Z radius parent +1 1 0.000000000 0.000000000 0.000000000 0.500000000 -1 +2 3 -1.000000000 0.000000000 0.000000000 0.100000001 1 +3 3 -2.000000000 0.000000000 0.000000000 0.100000001 2 +4 3 -3.000000000 0.000000000 0.000000000 0.100000001 3 +5 3 -3.000000000 0.000000000 1.000000000 0.100000001 4 +6 3 -3.000000000 0.000000000 -1.000000000 0.100000001 4 +7 3 -2.000000000 1.000000000 0.000000000 0.100000001 3 +8 3 0.000000000 1.000000000 0.000000000 0.100000001 1 +9 3 1.000000000 2.000000000 0.000000000 0.100000001 8 +10 3 1.000000000 4.000000000 0.000000000 0.100000001 9 +11 3 1.000000000 4.000000000 1.000000000 0.100000001 10 +12 3 1.000000000 4.000000000 -1.000000000 0.100000001 10 +13 2 2.000000000 3.000000000 0.000000000 0.100000001 9 +14 2 2.000000000 4.000000000 0.000000000 0.100000001 13 +15 2 3.000000000 3.000000000 0.000000000 0.100000001 13 +16 2 3.000000000 3.000000000 1.000000000 0.100000001 15 +17 2 3.000000000 3.000000000 -1.000000000 0.100000001 15 +18 4 0.000000000 -1.000000000 0.000000000 0.100000001 1 +19 4 0.000000000 -2.000000000 0.000000000 0.100000001 18 +20 4 0.000000000 -3.000000000 0.000000000 0.100000001 19 +21 4 0.000000000 -3.000000000 1.000000000 0.100000001 20 +22 4 0.000000000 -3.000000000 -1.000000000 0.100000001 20 +23 4 1.000000000 -2.000000000 0.000000000 0.100000001 19 diff --git a/tests/features/test_get_features.py b/tests/features/test_get_features.py index b0ac1be0..65ab30d5 100644 --- a/tests/features/test_get_features.py +++ b/tests/features/test_get_features.py @@ -644,23 +644,33 @@ def test_section_radial_distances_origin(): def test_number_of_sections_per_neurite(): - nsecs = features.get('number_of_sections_per_neurite', NEURON) - assert len(nsecs) == 4 - assert np.all(nsecs == [21, 21, 21, 21]) - - nsecs = features.get('number_of_sections_per_neurite', NEURON, neurite_type=NeuriteType.axon) - assert len(nsecs) == 1 - assert nsecs == [21] - - nsecs = features.get('number_of_sections_per_neurite', NEURON, - neurite_type=NeuriteType.basal_dendrite) - assert len(nsecs) == 2 - assert np.all(nsecs == [21, 21]) - - nsecs = features.get('number_of_sections_per_neurite', NEURON, - neurite_type=NeuriteType.apical_dendrite) - assert len(nsecs) == 1 - assert np.all(nsecs == [21]) + for use_subtrees in (True, False): + nsecs = features.get('number_of_sections_per_neurite', + NEURON, + use_subtrees=use_subtrees) + assert len(nsecs) == 4 + assert np.all(nsecs == [21, 21, 21, 21]) + + nsecs = features.get('number_of_sections_per_neurite', + NEURON, + neurite_type=NeuriteType.axon, + use_subtrees=use_subtrees) + assert len(nsecs) == 1 + assert nsecs == [21] + + nsecs = features.get('number_of_sections_per_neurite', + NEURON, + neurite_type=NeuriteType.basal_dendrite, + use_subtrees=use_subtrees) + assert len(nsecs) == 2 + assert np.all(nsecs == [21, 21]) + + nsecs = features.get('number_of_sections_per_neurite', + NEURON, + neurite_type=NeuriteType.apical_dendrite, + use_subtrees=use_subtrees) + assert len(nsecs) == 1 + assert np.all(nsecs == [21]) def test_trunk_origin_radii(): diff --git a/tests/features/test_section.py b/tests/features/test_section.py index d549410a..8d1f478e 100644 --- a/tests/features/test_section.py +++ b/tests/features/test_section.py @@ -37,7 +37,6 @@ import pytest import numpy as np from numpy import testing as npt -from mock import Mock from neurom import load_morphology, iter_sections from neurom import morphmath @@ -74,6 +73,21 @@ def test_segment_taper_rates(): sec = Mock(points=np.array([[0., 0., 0., 2.], [1., 0., 0., 1.], [2., 0., 0., 0.]])) npt.assert_almost_equal(section.segment_taper_rates(sec), [-2., -2.]) +def test_section_path_length(): + m = load_morphology( + """ + 1 1 0 0 0 0.5 -1 + 2 3 1 0 0 0.1 1 + 3 3 2 0 0 0.1 2 + 4 3 3 0 0 0.1 3 + 5 3 2 1 0 0.1 3 + """, + reader="swc", + ) + + sec = m.sections[1] + npt.assert_almost_equal(section.section_path_length(sec), 2.0) + def test_section_area(): sec = load_morphology(StringIO(u"""((CellBody) (0 0 0 2)) diff --git a/tests/test_mixed.py b/tests/test_mixed.py new file mode 100644 index 00000000..9cb11822 --- /dev/null +++ b/tests/test_mixed.py @@ -0,0 +1,2157 @@ +import sys +import warnings +import pytest +import neurom +import numpy as np +import numpy.testing as npt +from neurom import NeuriteType +from neurom.features import get +from neurom.core import Population +from neurom.features import _POPULATION_FEATURES, _MORPHOLOGY_FEATURES, _NEURITE_FEATURES +import collections.abc + +from neurom.core.types import tree_type_checker as is_type + +import neurom.core.morphology +import neurom.features.neurite + + +@pytest.fixture +def mixed_morph(): + """ + (1, 4, 1) + | + S7:B | + | + (1, 4, -1)-----(1, 4, 0) (2, 4, 0) (3, 3, 1) + S8:B | | | + | S10:A | S12:A | + | | S11:A | + S6:B | (2, 3, 0)-----(3, 3, 0) + | / | + | S9:A / S13:A | + | / | + (1, 2, 0) (3, 3, -1) + / + S5:B / + / Axon on basal dendrite + (-3, 0, 1) (-2, 1, 0) (0, 1, 0) + | | + S2 | S4 | + | S1 | S0 + (-3, 0, 0)-----(-2, 0, 0)-----(-1, 0, 0) (0, 0, 0) Soma + | + S3 | Basal Dendrite + | + (-3, 0, -1) (0, -1, 0) + | + S14 | + | S17 + Apical Dendrite (0, -2, 0)-----(1, -2, 0) + | + S15 | + S17 | S16 + (0, -3, -1)-----(0, -3, 0)-----(0, -3, 1) + + basal_dendrite: homogeneous + section ids: [0, 1, 2, 3, 4] + + axon_on_basal_dendrite: heterogeneous + section_ids: + - basal: [5, 6, 7, 8] + - axon : [9, 10, 11, 12, 13] + + apical_dendrite: homogeneous: + section_ids: [14, 15, 16, 17, 18] + """ + return neurom.load_morphology( + """ + 1 1 0 0 0 0.5 -1 + 2 3 -1 0 0 0.1 1 + 3 3 -2 0 0 0.1 2 + 4 3 -3 0 0 0.1 3 + 5 3 -3 0 1 0.1 4 + 6 3 -3 0 -1 0.1 4 + 7 3 -2 1 0 0.1 3 + 8 3 0 1 0 0.1 1 + 9 3 1 2 0 0.1 8 + 10 3 1 4 0 0.1 9 + 11 3 1 4 1 0.1 10 + 12 3 1 4 -1 0.1 10 + 13 2 2 3 0 0.1 9 + 14 2 2 4 0 0.1 13 + 15 2 3 3 0 0.1 13 + 16 2 3 3 1 0.1 15 + 17 2 3 3 -1 0.1 15 + 18 4 0 -1 0 0.1 1 + 19 4 0 -2 0 0.1 18 + 20 4 0 -3 0 0.1 19 + 21 4 0 -3 1 0.1 20 + 22 4 0 -3 -1 0.1 20 + 23 4 1 -2 0 0.1 19 + """, + reader="swc") + +@pytest.fixture +def three_types_neurite_morph(): + return neurom.load_morphology( + """ + 1 1 0 0 0 0.5 -1 + 2 3 0 1 0 0.1 1 + 3 3 1 2 0 0.1 2 + 4 3 1 4 0 0.1 3 + 5 3 1 4 1 0.1 4 + 6 3 1 4 -1 0.1 4 + 7 2 2 3 0 0.1 3 + 8 2 2 4 0 0.1 7 + 9 2 3 3 0 0.1 7 + 10 2 3 3 1 0.1 9 + 11 4 3 3 -1 0.1 9 + """, + reader="swc") + +def test_heterogeneous_neurites(mixed_morph): + + assert not mixed_morph.neurites[0].is_heterogeneous() + assert mixed_morph.neurites[1].is_heterogeneous() + assert not mixed_morph.neurites[2].is_heterogeneous() + + +def test_is_homogeneous_point(mixed_morph): + + heterogeneous_neurite = mixed_morph.neurites[1] + + sections = list(heterogeneous_neurite.iter_sections()) + + # first section has one axon and one basal children + assert not sections[0].is_homogeneous_point() + + # second section is pure basal + assert sections[1].is_homogeneous_point() + + +def test_homogeneous_subtrees(mixed_morph, three_types_neurite_morph): + + basal, axon_on_basal, apical = mixed_morph.neurites + + assert neurom.core.morphology._homogeneous_subtrees(basal) == [basal] + + sections = list(axon_on_basal.iter_sections()) + + subtrees = neurom.core.morphology._homogeneous_subtrees(axon_on_basal) + + assert subtrees[0].root_node.id == axon_on_basal.root_node.id + assert subtrees[0].root_node.type == NeuriteType.basal_dendrite + + assert subtrees[1].root_node.id == sections[4].id + assert subtrees[1].root_node.type == NeuriteType.axon + + with pytest.warns( + UserWarning, + match="Neurite is not an axon-carrying dendrite." + ): + three_types_neurite, = three_types_neurite_morph.neurites + neurom.core.morphology._homogeneous_subtrees(three_types_neurite) + + +def test_iter_neurites__heterogeneous(mixed_morph): + + subtrees = list(neurom.core.morphology.iter_neurites(mixed_morph, use_subtrees=False)) + + assert len(subtrees) == 3 + assert subtrees[0].type == NeuriteType.basal_dendrite + assert subtrees[1].type == NeuriteType.basal_dendrite + assert subtrees[2].type == NeuriteType.apical_dendrite + + subtrees = list(neurom.core.morphology.iter_neurites(mixed_morph, use_subtrees=True)) + + assert len(subtrees) == 4 + assert subtrees[0].type == NeuriteType.basal_dendrite + assert subtrees[1].type == NeuriteType.basal_dendrite + assert subtrees[2].type == NeuriteType.axon + assert subtrees[3].type == NeuriteType.apical_dendrite + + +def test_core_iter_sections__heterogeneous(mixed_morph): + + def assert_sections(neurite, section_type, expected_section_ids): + + it = neurom.core.morphology.iter_sections(neurite, section_filter=is_type(section_type)) + assert [s.id for s in it] == expected_section_ids + + basal, axon_on_basal, apical = mixed_morph.neurites + + assert_sections(basal, NeuriteType.all, [0, 1, 2, 3, 4]) + assert_sections(basal, NeuriteType.basal_dendrite, [0, 1, 2, 3, 4]) + assert_sections(basal, NeuriteType.axon, []) + + assert_sections(axon_on_basal, NeuriteType.all, [5, 6, 7, 8, 9, 10, 11, 12, 13]) + assert_sections(axon_on_basal, NeuriteType.basal_dendrite, [5, 6, 7, 8]) + assert_sections(axon_on_basal, NeuriteType.axon, [9, 10, 11, 12, 13]) + + assert_sections(apical, NeuriteType.all, [14, 15, 16, 17, 18]) + assert_sections(apical, NeuriteType.apical_dendrite, [14, 15, 16, 17, 18]) + + +def test_features_neurite_map_sections__heterogeneous(mixed_morph): + + def assert_sections(neurite, section_type, iterator_type, expected_section_ids): + function = lambda section: section.id + section_ids = neurom.features.neurite._map_sections( + function, neurite, iterator_type=iterator_type, section_type=section_type + ) + assert section_ids == expected_section_ids + + basal, axon_on_basal, apical = mixed_morph.neurites + + # homogeneous tree, no difference between all and basal_dendrite types. + assert_sections( + basal, NeuriteType.all, neurom.core.morphology.Section.ibifurcation_point, + [0, 1], + ) + assert_sections( + basal, NeuriteType.basal_dendrite, neurom.core.morphology.Section.ibifurcation_point, + [0, 1], + ) + # heterogeneous tree, forks cannot be heterogeneous if a type other than all is specified + # Section with id 5 is the transition section, which has a basal and axon children sections + assert_sections( + axon_on_basal, NeuriteType.all, neurom.core.morphology.Section.ibifurcation_point, + [5, 6, 9, 11], + ) + assert_sections( + axon_on_basal, NeuriteType.basal_dendrite, + neurom.core.morphology.Section.ibifurcation_point, + [6], + ) + assert_sections( + axon_on_basal, NeuriteType.axon, + neurom.core.morphology.Section.ibifurcation_point, + [9, 11], + ) + # homogeneous tree, no difference between all and basal_dendrite types. + assert_sections( + apical, NeuriteType.all, neurom.core.morphology.Section.ibifurcation_point, + [14, 15], + ) + assert_sections( + apical, NeuriteType.apical_dendrite, neurom.core.morphology.Section.ibifurcation_point, + [14, 15], + ) + + +@pytest.fixture +def population(mixed_morph): + return Population([mixed_morph, mixed_morph]) + + +def _assert_feature_equal(obj, feature_name, expected_values, kwargs, use_subtrees): + + def innermost_value(iterable): + while isinstance(iterable, collections.abc.Iterable): + try: + iterable = iterable[0] + except IndexError: + # empty list + return None + return iterable + + + assert_equal = lambda a, b: npt.assert_equal( + a, b, err_msg=f"ACTUAL: {a}\nDESIRED: {b}", verbose=False + ) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + + values = get(feature_name, obj, use_subtrees=use_subtrees, **kwargs) + # handle empty lists because allclose always passes in that case. + # See: https://github.com/numpy/numpy/issues/11071 + if isinstance(values, collections.abc.Iterable): + if isinstance(expected_values, collections.abc.Iterable): + if isinstance(innermost_value(values), (float, np.floating)): + npt.assert_allclose(values, expected_values, atol=1e-5) + else: + assert_equal(values, expected_values) + else: + assert_equal(values, expected_values) + else: + if isinstance(expected_values, collections.abc.Iterable): + assert_equal(values, expected_values) + else: + if isinstance(values, (float, np.floating)): + npt.assert_allclose(values, expected_values, atol=1e-5) + else: + assert_equal(values, expected_values) + + +def _dispatch_features(features, mode): + for feature_name, configurations in features.items(): + + for cfg in configurations: + kwargs = cfg["kwargs"] if "kwargs" in cfg else {} + + if mode == "with-subtrees": + expected = cfg["expected_with_subtrees"] + elif mode == "wout-subtrees": + expected = cfg["expected_wout_subtrees"] + else: + raise ValueError("Uknown mode") + + yield feature_name, kwargs, expected + + +def _population_features(mode): + + features = { + "sholl_frequency": [ + { + "kwargs": {"neurite_type": NeuriteType.all, "step_size": 3}, + "expected_wout_subtrees": [0, 4], + "expected_with_subtrees": [0, 4], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite, "step_size": 3}, + "expected_wout_subtrees": [0, 4], + "expected_with_subtrees": [0, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon, "step_size": 3}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite, "step_size": 2}, + "expected_wout_subtrees": [0, 2], + "expected_with_subtrees": [0, 2], + }, + + ], + } + + features_not_tested = list( + set(_POPULATION_FEATURES) - set(features.keys()) + ) + + assert not features_not_tested, ( + "The following morphology tests need to be included in the tests:\n\n" + + "\n".join(sorted(features_not_tested)) + "\n" + ) + + return _dispatch_features(features, mode) + + +@pytest.mark.parametrize("feature_name, kwargs, expected", _population_features(mode="wout-subtrees")) +def test_population__population_features_wout_subtrees(feature_name, kwargs, expected, population): + _assert_feature_equal(population, feature_name, expected, kwargs, use_subtrees=False) + + +@pytest.mark.parametrize("feature_name, kwargs, expected", _population_features(mode="with-subtrees")) +def test_population__population_features_with_subtrees(feature_name, kwargs, expected, population): + _assert_feature_equal(population, feature_name, expected, kwargs, use_subtrees=True) + + +def _morphology_features(mode): + + features = { + "soma_radius": [ + { + "expected_wout_subtrees": 0.5, + "expected_with_subtrees": 0.5, + } + ], + "soma_surface_area": [ + { + "expected_wout_subtrees": np.pi, + "expected_with_subtrees": np.pi, + } + ], + "soma_volume": [ + { + "expected_wout_subtrees": np.pi / 6., + "expected_with_subtrees": np.pi / 6., + } + ], + "number_of_sections_per_neurite": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [5, 9, 5], + "expected_with_subtrees": [5, 4, 5, 5], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [5, 9], + "expected_with_subtrees": [5, 4], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [5], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [5], + "expected_with_subtrees": [5], + } + ], + "max_radial_distance": [ + { + # without subtrees AoD is considered a single tree, with [3, 3] being the furthest + # with subtrees AoD subtrees are considered separately and the distance is calculated + # from their respective roots. [1, 4] is the furthest point in this case + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 3.741657, + "expected_with_subtrees": 3.316625, + }, + { + # with a global origin, AoD axon subtree [2, 4] is always furthest from soma + "kwargs": {"neurite_type": NeuriteType.all, "origin": np.array([0., 0., 0.])}, + "expected_wout_subtrees": 4.47213595499958, + "expected_with_subtrees": 4.47213595499958, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 3.741657, + "expected_with_subtrees": 3.316625, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite, "origin": np.array([0., 0., 0.])}, + "expected_wout_subtrees": 4.472136, + "expected_with_subtrees": 4.242641, + + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0.0, + "expected_with_subtrees": 2.44949, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon, "origin": np.array([0., 0., 0.])}, + "expected_wout_subtrees": 0.0, + "expected_with_subtrees": 4.47213595499958, + } + ], + "total_length_per_neurite": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [5., 10.828427, 5.], + "expected_with_subtrees": [5., 5.414213, 5.414213, 5.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [5., 10.828427], + "expected_with_subtrees": [5., 5.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [5.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [5.], + "expected_with_subtrees": [5.], + } + ], + "total_area_per_neurite" : [ + { + # total length x 2piR + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [3.141593, 6.803702, 3.141593], + "expected_with_subtrees": [3.141593, 3.401851, 3.401851, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [3.141593, 6.803702], + "expected_with_subtrees": [3.141593, 3.401851], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [3.401851], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [3.141593], + "expected_with_subtrees": [3.141593], + } + ], + "total_volume_per_neurite": [ + # total_length * piR^2 + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.15708 , 0.340185, 0.15708 ], + "expected_with_subtrees": [0.15708 , 0.170093, 0.170093, 0.15708], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.15708 , 0.340185], + "expected_with_subtrees": [0.15708 , 0.170093], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.170093], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.15708], + "expected_with_subtrees": [0.15708], + } + ], + "trunk_origin_azimuths": [ # Not applicable to distal subtrees + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [3.1415927, 0.0, 0.0], + "expected_with_subtrees": [3.1415927, 0.0, 0.0], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [3.1415927, 0.0], + "expected_with_subtrees": [3.1415927, 0.0], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.0], + "expected_with_subtrees": [0.0], + }, + ], + "trunk_origin_elevations": [ # Not applicable to distal subtrees + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.0, 1.5707964, -1.5707964], + "expected_with_subtrees": [0.0, 1.5707964, -1.5707964], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.0, 1.5707964], + "expected_with_subtrees": [0.0, 1.5707964], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [-1.570796], + "expected_with_subtrees": [-1.570796], + }, + ], + "trunk_vectors": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [[-1., 0., 0.], [0., 1., 0.], [0., -1., 0.]], + "expected_with_subtrees": [[-1., 0., 0.], [0., 1., 0.], [1., 2., 0.], [0., -1., 0.]], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [[-1., 0., 0.], [0., 1., 0.]], + "expected_with_subtrees": [[-1., 0., 0.], [0., 1., 0.]], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [[1., 2., 0.]], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [[0., -1., 0.]], + "expected_with_subtrees": [[0., -1., 0.]], + }, + + ], + "trunk_angles": [ # Not applicable to distal subtrees + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [1.570796, 3.141592, 1.570796], + "expected_with_subtrees": [1.570796, 3.141592, 1.570796], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [1.5707964, 1.570796], + "expected_with_subtrees": [1.5707964, 1.570796], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.], + "expected_with_subtrees": [0.], + }, + ], + "trunk_angles_from_vector": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [ + [np.pi / 2., - np.pi / 2, np.pi], + [0., 0., 0.], + [np.pi, np.pi, 0.], + ], + "expected_with_subtrees": [ + [np.pi / 2., - np.pi / 2, np.pi], + [0., 0., 0.], + [0.463648, -0.463648, 0.], + [np.pi, np.pi, 0.], + ], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [[np.pi / 2., - np.pi / 2, np.pi], [0., 0., 0.]], + "expected_with_subtrees": [[np.pi / 2., - np.pi / 2, np.pi], [0., 0., 0.]], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [[0.463648, -0.463648, 0.]], + }, + + ], + "trunk_angles_inter_types": [ + { + "kwargs": { + "source_neurite_type": NeuriteType.basal_dendrite, + "target_neurite_type": NeuriteType.axon, + }, + "expected_wout_subtrees": [], + "expected_with_subtrees": [ + [[ 2.034444, 1.107149, -3.141593]], + [[ 0.463648, -0.463648, 0. ]], + ], + }, + ], + "trunk_origin_radii": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.1, 0.1, 0.1], + "expected_with_subtrees": [0.1, 0.1, 0.1, 0.1], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.1, 0.1], + "expected_with_subtrees": [0.1, 0.1], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.1], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.1], + "expected_with_subtrees": [0.1], + }, + ], + "trunk_section_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [1., 1.414213, 1.], + "expected_with_subtrees": [1., 1.414213, 1.414213, 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [1., 1.414213], + "expected_with_subtrees": [1., 1.414213], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414213], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.], + "expected_with_subtrees": [1.], + }, + ], + "number_of_neurites": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 3, + "expected_with_subtrees": 4, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 2, + "expected_with_subtrees": 2, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0, + "expected_with_subtrees": 1, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 1, + "expected_with_subtrees": 1, + }, + ], + "neurite_volume_density": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.235619, 0.063785, 0.235619], + "expected_with_subtrees": [0.235619, 0.255139, 0.170093, 0.235619], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.235619, 0.063785], + "expected_with_subtrees": [0.235619, 0.255139], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.170093], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.235619], + "expected_with_subtrees": [0.235619], + }, + ], + "sholl_crossings": [ + { + "kwargs": {"neurite_type": NeuriteType.all, "radii": [1.5, 3.5]}, + "expected_wout_subtrees": [3, 2], + "expected_with_subtrees": [3, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite, "radii": [1.5, 3.5]}, + "expected_wout_subtrees": [2, 2], + "expected_with_subtrees": [2, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon, "radii": [1.5, 3.5]}, + "expected_wout_subtrees": [0, 0], + "expected_with_subtrees": [0, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite, "radii": [1.5, 3.5]}, + "expected_wout_subtrees": [1, 0], + "expected_with_subtrees": [1, 0], + }, + ], + "sholl_frequency": [ + { + "kwargs": {"neurite_type": NeuriteType.all, "step_size": 3}, + "expected_wout_subtrees": [0, 2], + "expected_with_subtrees": [0, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite, "step_size": 3}, + "expected_wout_subtrees": [0, 2], + "expected_with_subtrees": [0, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon, "step_size": 3}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite, "step_size": 2}, + "expected_wout_subtrees": [0, 1], + "expected_with_subtrees": [0, 1], + }, + + ], + "total_width": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 6.0, + "expected_with_subtrees": 6.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 6.0, + "expected_with_subtrees": 4.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0.0, + "expected_with_subtrees": 2.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 1.0, + "expected_with_subtrees": 1.0, + }, + ], + "total_height": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 7.0, + "expected_with_subtrees": 7.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 4.0, + "expected_with_subtrees": 4.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0.0, + "expected_with_subtrees": 2.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 2.0, + "expected_with_subtrees": 2.0, + }, + ], + "total_depth": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 2.0, + "expected_with_subtrees": 2.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 2.0, + "expected_with_subtrees": 2.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0.0, + "expected_with_subtrees": 2.0, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 2.0, + "expected_with_subtrees": 2.0, + }, + ], + "volume_density": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 0.01570426, + "expected_with_subtrees": 0.01570426, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 0.02983588, + "expected_with_subtrees": 0.04907583, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": np.nan, + "expected_with_subtrees": 0.17009254, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 0.23561945, + "expected_with_subtrees": 0.23561945, + }, + ], + "aspect_ratio":[ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 0.630311, + "expected_with_subtrees": 0.630311, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 0.305701, + "expected_with_subtrees": 0.284467, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": np.nan, + "expected_with_subtrees": 0.666667, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 0.5, + "expected_with_subtrees": 0.5, + }, + ], + "circularity": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 0.739583, + "expected_with_subtrees": 0.739583, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 0.525588, + "expected_with_subtrees": 0.483687, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": np.nan, + "expected_with_subtrees": 0.544013, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 0.539012, + "expected_with_subtrees": 0.539012, + }, + ], + "shape_factor": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 0.40566, + "expected_with_subtrees": 0.40566, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 0.21111, + "expected_with_subtrees": 0.18750, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": np.nan, + "expected_with_subtrees": 0.3, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 0.25, + "expected_with_subtrees": 0.25, + }, + ], + "length_fraction_above_soma": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 0.567898, + "expected_with_subtrees": 0.567898, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 0.61591, + "expected_with_subtrees": 0.74729, + }, + ], + } + + features_not_tested = set(_MORPHOLOGY_FEATURES) - set(features.keys()) + + assert not features_not_tested, ( + "The following morphology tests need to be included in the mixed morphology tests:\n" + f"{features_not_tested}" + ) + + return _dispatch_features(features, mode) + + +@pytest.mark.parametrize("feature_name, kwargs, expected", _morphology_features(mode="wout-subtrees")) +def test_morphology__morphology_features_wout_subtrees(feature_name, kwargs, expected, mixed_morph): + _assert_feature_equal(mixed_morph, feature_name, expected, kwargs, use_subtrees=False) + + +@pytest.mark.parametrize("feature_name, kwargs, expected", _morphology_features(mode="with-subtrees")) +def test_morphology__morphology_features_with_subtrees( + feature_name, kwargs, expected, mixed_morph +): + _assert_feature_equal(mixed_morph, feature_name, expected, kwargs, use_subtrees=True) + + +def _neurite_features(mode): + + features = { + "number_of_segments": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 19, + "expected_with_subtrees": 19, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 14, + "expected_with_subtrees": 9, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0, + "expected_with_subtrees": 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 5, + "expected_with_subtrees": 5, + }, + ], + "number_of_leaves": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 11, + "expected_with_subtrees": 11, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 8, + "expected_with_subtrees": 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0, + "expected_with_subtrees": 3, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 3, + "expected_with_subtrees": 3, + }, + ], + "total_length": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 20.828427, + "expected_with_subtrees": 20.828427, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 15.828427, + "expected_with_subtrees": 10.414214, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0.0, + "expected_with_subtrees": 5.414214, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 5., + "expected_with_subtrees": 5., + } + ], + "total_area": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 13.086887, + "expected_with_subtrees": 13.086887, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 9.945294, + "expected_with_subtrees": 6.543443, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0., + "expected_with_subtrees": 3.401851, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 3.141593, + "expected_with_subtrees": 3.141593, + } + ], + "total_volume": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 0.654344, + "expected_with_subtrees": 0.654344, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 0.497265, + "expected_with_subtrees": 0.327172, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0., + "expected_with_subtrees": 0.170093, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 0.15708, + "expected_with_subtrees": 0.15708, + } + ], + "section_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1.] * 5 + [1.414214, 2., 1., 1.] + [1.414214, 1., 1., 1., 1.] + [1.] * 5, + "expected_with_subtrees": + [1.] * 5 + [1.414214, 2., 1., 1.] + [1.414214, 1., 1., 1., 1.] + [1.] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1.] * 5 + [1.414214, 2., 1., 1.] + [1.414214, 1., 1., 1., 1], + "expected_with_subtrees": + [1.] * 5 + [1.414214, 2., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 1., 1., 1., 1.], + + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.] * 5, + "expected_with_subtrees": [1.] * 5, + } + ], + "section_areas": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [0.628318] * 5 + [0.888577, 1.256637, 0.628319, 0.628319] + + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319] + [0.628318] * 5, + "expected_with_subtrees": + [0.628318] * 5 + [0.888577, 1.256637, 0.628319, 0.628319] + + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319] + [0.628318] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [0.628318] * 5 + [0.888577, 1.256637, 0.628319, 0.628319] + + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319], + "expected_with_subtrees": + [0.628318] * 5 + [0.888577, 1.256637, 0.628319, 0.628319], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.888577, 0.628319, 0.628319, 0.628319, 0.628319], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.628318] * 5, + "expected_with_subtrees": [0.628318] * 5, + } + + ], + "section_volumes": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [0.031416] * 5 + [0.044429, 0.062832, 0.031416, 0.031416] + + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416] + + [0.031416] * 5, + "expected_with_subtrees": + [0.031416] * 5 + [0.044429, 0.062832, 0.031416, 0.031416] + + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416] + + [0.031416] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [0.031416] * 5 + [0.044429, 0.062832, 0.031416, 0.031416] + + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416], + "expected_with_subtrees": + [0.031416] * 5 + [0.044429, 0.062832, 0.031416, 0.031416], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.044429, 0.031416, 0.031416, 0.031416, 0.031416], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.031415] * 5, + "expected_with_subtrees": [0.031415] * 5, + } + ], + "section_tortuosity": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [1.0] * 19, + "expected_with_subtrees": [1.0] * 19, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [1.0] * 14, + "expected_with_subtrees": [1.0] * 9, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.0] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.0] * 5, + "expected_with_subtrees": [1.0] * 5, + } + ], + "section_radial_distances": [ + { + # radial distances change when the mixed subtrees are processed because + # the root of the subtree is considered + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1. , 2. , 2.236068, 2.236068, 1.414214] + + [1.414214, 3.162278, 3.316625, 3.316625] + + [2.828427, 3.605551, 3.605551, 3.741657, 3.741657] + + [1., 2., 2.236068, 2.236068, 1.414214], + "expected_with_subtrees": + [1. , 2. , 2.236068, 2.236068, 1.414214] + + [1.414214, 3.162278, 3.316625, 3.316625] + + [1.414214, 2.236068, 2.236068, 2.44949 , 2.44949] + + [1., 2., 2.236068, 2.236068, 1.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1. , 2. , 2.236068, 2.236068, 1.414214] + + [1.414214, 3.162278, 3.316625, 3.316625] + + [2.828427, 3.605551, 3.605551, 3.741657, 3.741657], + "expected_with_subtrees": + [1. , 2. , 2.236068, 2.236068, 1.414214] + + [1.414214, 3.162278, 3.316625, 3.316625], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 2.236068, 2.236068, 2.44949 , 2.44949], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1., 2., 2.236068, 2.236068, 1.414214], + "expected_with_subtrees": [1., 2., 2.236068, 2.236068, 1.414214], + } + + ], + "section_term_radial_distances": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [2.236068, 2.236068, 1.414214] + + [3.316625, 3.316625] + + [3.605551, 3.741657, 3.741657] + + [2.236068, 2.236068, 1.414214], + "expected_with_subtrees": + [2.236068, 2.236068, 1.414214] + + [3.316625, 3.316625] + + [2.236068, 2.44949 , 2.44949] + + [2.236068, 2.236068, 1.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [2.236068, 2.236068, 1.414214] + + [3.316625, 3.316625] + + [3.605551, 3.741657, 3.741657], + "expected_with_subtrees": + [2.236068, 2.236068, 1.414214] + + [3.316625, 3.316625], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [2.236068, 2.44949 , 2.44949], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [2.236068, 2.236068, 1.414214], + "expected_with_subtrees": [2.236068, 2.236068, 1.414214], + } + + ], + "section_bif_radial_distances": [ + { + # radial distances change when the mixed subtrees are processed because + # the root of the subtree is considered instead of the tree root + # heterogeneous forks are not valid forking points + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1., 2., 1.414214, 3.162278, 2.828427, 3.605551, 1., 2.], + "expected_with_subtrees": + [1., 2., 3.162278, 1.414214, 2.236068, 1., 2.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1., 2., 1.414214, 3.162278, 2.828427, 3.605551], + "expected_with_subtrees": [1., 2., 3.162278], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 2.236068], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1., 2.], + "expected_with_subtrees": [1., 2.], + } + ], + "section_end_distances": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.] + + [1.414214, 1., 1., 1., 1.] + + [1.] * 5, + "expected_with_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.] + + [1.414214, 1., 1., 1., 1.] + + [1.] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.] + + [1.414214, 1., 1., 1., 1.], + "expected_with_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 1., 1., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.] * 5, + "expected_with_subtrees": [1.] * 5, + } + ], + "section_term_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [1.] * 11, + "expected_with_subtrees": [1.] * 11, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [1.] * 8, + "expected_with_subtrees": [1.] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.] * 3, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.] * 3, + "expected_with_subtrees": [1.] * 3, + } + ], + "section_taper_rates": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.0] * 19, + "expected_with_subtrees": [0.0] * 19, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.0] * 14, + "expected_with_subtrees": [0.0] * 9, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.0] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.0] * 5, + "expected_with_subtrees": [0.0] * 5, + } + ], + "section_bif_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1., 1., 1.414214, 2., 1.414214, 1., 1., 1.], + "expected_with_subtrees": + [1., 1., 2., 1.414214, 1., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [1., 1., 1.414214, 2., 1.414214, 1.], + "expected_with_subtrees": [1., 1., 2.], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1., 1.], + "expected_with_subtrees": [1., 1.], + }, + ], + "section_branch_orders": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [0, 1, 2, 2, 1, 0, 1, 2, 2, 1, 2, 2, 3, 3, 0, 1, 2, 2, 1], + "expected_with_subtrees": + [0, 1, 2, 2, 1, 0, 1, 2, 2, 1, 2, 2, 3, 3, 0, 1, 2, 2, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0, 1, 2, 2, 1, 0, 1, 2, 2, 1, 2, 2, 3, 3], + "expected_with_subtrees": [0, 1, 2, 2, 1, 0, 1, 2, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1, 2, 2, 3, 3], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0, 1, 2, 2, 1], + "expected_with_subtrees": [0, 1, 2, 2, 1], + }, + ], + "section_bif_branch_orders": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0, 1, 0, 1, 1, 2, 0, 1], + "expected_with_subtrees": [0, 1, 1, 1, 2, 0, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0, 1, 0, 1, 1, 2], + "expected_with_subtrees": [0, 1, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0, 1], + "expected_with_subtrees": [0, 1], + }, + ], + "section_term_branch_orders": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [2, 2, 1, 2, 2, 2, 3, 3, 2, 2, 1], + "expected_with_subtrees": [2, 2, 1, 2, 2, 2, 3, 3, 2, 2, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [2, 2, 1, 2, 2, 2, 3, 3], + "expected_with_subtrees": [2, 2, 1, 2, 2], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [2, 3, 3], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [2, 2, 1], + "expected_with_subtrees": [2, 2, 1], + }, + ], + "section_strahler_orders": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [2, 2, 1, 1, 1, 3, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 1, 1, 1], + "expected_with_subtrees": + [2, 2, 1, 1, 1, 3, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 1, 1, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [2, 2, 1, 1, 1, 3, 2, 1, 1, 2, 1, 2, 1, 1], + "expected_with_subtrees": [2, 2, 1, 1, 1, 3, 2, 1, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [2, 1, 2, 1, 1], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [2, 2, 1, 1, 1], + "expected_with_subtrees": [2, 2, 1, 1, 1], + }, + ], + "segment_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.] + + [1.414214, 1., 1., 1., 1.] + + [1.] * 5, + "expected_with_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.] + + [1.414214, 1., 1., 1., 1.] + + [1.] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.] + + [1.414214, 1., 1., 1., 1.], + "expected_with_subtrees": + [1.] * 5 + + [1.414214, 2., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 1., 1., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.] * 5, + "expected_with_subtrees": [1.] * 5, + } + ], + "segment_areas": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [0.628319] * 5 + + [0.888577, 1.256637, 0.628319, 0.628319] + + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319] + + [0.628319] * 5, + "expected_with_subtrees": + [0.628319] * 5 + + [0.888577, 1.256637, 0.628319, 0.628319] + + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319] + + [0.628319] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [0.628319] * 5 + + [0.888577, 1.256637, 0.628319, 0.628319] + + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319], + "expected_with_subtrees": + [0.628319] * 5 + + [0.888577, 1.256637, 0.628319, 0.628319], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": + [0.888577, 0.628319, 0.628319, 0.628319, 0.628319], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.628318] * 5, + "expected_with_subtrees": [0.628318] * 5, + } + ], + "segment_volumes": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [0.031415] * 5 + + [0.044429, 0.062832, 0.031416, 0.031416] + + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416] + + [0.031416] * 5, + "expected_with_subtrees": + [0.031415] * 5 + + [0.044429, 0.062832, 0.031416, 0.031416] + + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416] + + [0.031416] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [0.031415] * 5 + + [0.044429, 0.062832, 0.031416, 0.031416] + + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416], + "expected_with_subtrees": + [0.031415] * 5 + + [0.044429, 0.062832, 0.031416, 0.031416], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": + [0.044429, 0.031416, 0.031416, 0.031416, 0.031416], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.031415] * 5, + "expected_with_subtrees": [0.031415] * 5, + } + ], + "segment_radii": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.1] * 19, + "expected_with_subtrees": [0.1] * 19, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.1] * 14, + "expected_with_subtrees": [0.1] * 9, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.1] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.1] * 5, + "expected_with_subtrees": [0.1] * 5, + } + ], + "segment_taper_rates": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [0.0] * 19, + "expected_with_subtrees": [0.0] * 19, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [0.0] * 14, + "expected_with_subtrees": [0.0] * 9, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.0] * 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.0] * 5, + "expected_with_subtrees": [0.0] * 5, + }, + ], + "segment_radial_distances": [ + { + # radial distances change when the mixed subtrees are processed because + # the root of the subtree is considered + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [0.5, 1.5, 2.061553, 2.061553, 1.118034] + + [0.707107, 2.236068, 3.201562, 3.201562] + + [2.12132 , 3.201562, 3.201562, 3.640055, 3.640055] + + [0.5, 1.5, 2.061553, 2.061553, 1.118034], + "expected_with_subtrees": + [0.5, 1.5, 2.061553, 2.061553, 1.118034] + + [0.707107, 2.236068, 3.201562, 3.201562] + + [0.707107, 1.802776, 1.802776, 2.291288, 2.291288] + + [0.5, 1.5, 2.061553, 2.061553, 1.118034], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [0.5, 1.5, 2.061553, 2.061553, 1.118034] + + [0.707107, 2.236068, 3.201562, 3.201562] + + [2.12132 , 3.201562, 3.201562, 3.640055, 3.640055], + "expected_with_subtrees": + [0.5, 1.5, 2.061553, 2.061553, 1.118034] + + [0.707107, 2.236068, 3.201562, 3.201562], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.707107, 1.802776, 1.802776, 2.291288, 2.291288], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [0.5, 1.5, 2.061553, 2.061553, 1.118034], + "expected_with_subtrees": [0.5, 1.5, 2.061553, 2.061553, 1.118034], + }, + ], + "segment_midpoints": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [ + [-1.5, 0.0, 0.0], [-2.5, 0.0, 0.0], [-3.0, 0.0, 0.5], [-3.0, 0.0, -0.5], + [-2.0, 0.5, 0.0], [0.5, 1.5, 0.0], [1.0, 3.0, 0.0], [1.0, 4.0, 0.5], + [1.0, 4.0, -0.5], [1.5, 2.5, 0.0], [2.0, 3.5, 0.0], [2.5, 3.0, 0.0], + [3.0, 3.0, 0.5], [3.0, 3.0, -0.5], [0.0, -1.5, 0.0], [0.0, -2.5, 0.0], + [0.0, -3.0, 0.5], [0.0, -3.0, -0.5], [0.5, -2.0, 0.0]], + "expected_with_subtrees": [ + [-1.5, 0.0, 0.0], [-2.5, 0.0, 0.0], [-3.0, 0.0, 0.5], [-3.0, 0.0, -0.5], + [-2.0, 0.5, 0.0], [0.5, 1.5, 0.0], [1.0, 3.0, 0.0], [1.0, 4.0, 0.5], + [1.0, 4.0, -0.5], [1.5, 2.5, 0.0], [2.0, 3.5, 0.0], [2.5, 3.0, 0.0], + [3.0, 3.0, 0.5], [3.0, 3.0, -0.5], [0.0, -1.5, 0.0], [0.0, -2.5, 0.0], + [0.0, -3.0, 0.5], [0.0, -3.0, -0.5], [0.5, -2.0, 0.0]], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [ + [-1.5, 0.0, 0.0], [-2.5, 0.0, 0.0], [-3.0, 0.0, 0.5], [-3.0, 0.0, -0.5], + [-2.0, 0.5, 0.0], [0.5, 1.5, 0.0], [1.0, 3.0, 0.0], [1.0, 4.0, 0.5], + [1.0, 4.0, -0.5], [1.5, 2.5, 0.0], [2.0, 3.5, 0.0], [2.5, 3.0, 0.0], + [3.0, 3.0, 0.5], [3.0, 3.0, -0.5]], + "expected_with_subtrees": [ + [-1.5, 0.0, 0.0], [-2.5, 0.0, 0.0], [-3.0, 0.0, 0.5], [-3.0, 0.0, -0.5], + [-2.0, 0.5, 0.0], [0.5, 1.5, 0.0], [1.0, 3.0, 0.0], [1.0, 4.0, 0.5], + [1.0, 4.0, -0.5]], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [ + [1.5, 2.5, 0.0], [2.0, 3.5, 0.0], [2.5, 3.0, 0.0], + [3.0, 3.0, 0.5], [3.0, 3.0, -0.5]], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [ + [0.0, -1.5, 0.0], [0.0, -2.5, 0.0], [0.0, -3.0, 0.5], + [0.0, -3.0, -0.5], [0.5, -2.0, 0.0]], + "expected_with_subtrees": [ + [0.0, -1.5, 0.0], [0.0, -2.5, 0.0], [0.0, -3.0, 0.5], + [0.0, -3.0, -0.5], [0.5, -2.0, 0.0]], + }, + ], + "segment_meander_angles": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [], + }, + ], + "number_of_sections": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 19, + "expected_with_subtrees": 19, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 14, + "expected_with_subtrees": 9, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0, + "expected_with_subtrees": 5, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 5, + "expected_with_subtrees": 5, + }, + ], + "number_of_bifurcations": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 8, + "expected_with_subtrees": 7, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 6, + "expected_with_subtrees": 3, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0, + "expected_with_subtrees": 2, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 2, + "expected_with_subtrees": 2, + }, + ], + "number_of_forking_points": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": 8, + "expected_with_subtrees": 7, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": 6, + "expected_with_subtrees": 3, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": 0, + "expected_with_subtrees": 2, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": 2, + "expected_with_subtrees": 2, + }, + ], + "local_bifurcation_angles": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1.570796, 3.141593, 0.785398, 3.141593, + 1.570796, 3.141593, 1.570796, 3.141593], + "expected_with_subtrees": + [1.570796, 3.141593, 3.141593, 1.570796, 3.141593, 1.570796, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1.570796, 3.141593, 0.785398, 3.141593, 1.570796, 3.141593], + "expected_with_subtrees": [1.570796, 3.141593, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.570796, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.570796, 3.141593], + "expected_with_subtrees": [1.570796, 3.141593], + }, + ], + "remote_bifurcation_angles": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1.570796, 3.141593, 0.785398, 3.141593, + 1.570796, 3.141593, 1.570796, 3.141593], + "expected_with_subtrees": + [1.570796, 3.141593, 3.141593, 1.570796, 3.141593, 1.570796, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [1.570796, 3.141593, 0.785398, 3.141593, 1.570796, 3.141593], + "expected_with_subtrees": [1.570796, 3.141593, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.570796, 3.141593], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.570796, 3.141593], + "expected_with_subtrees": [1.570796, 3.141593], + }, + ], + "sibling_ratios": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [1.0] * 8, + "expected_with_subtrees": [1.0] * 7, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [1.0] * 6, + "expected_with_subtrees": [1.0] * 3, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.0] * 2, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [1.0] * 2, + "expected_with_subtrees": [1.0] * 2, + }, + ], + "partition_pairs": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [[3.0, 1.0], [1.0, 1.0], [3.0, 5.0], + [1.0, 1.0], [1.0, 3.0], [1.0, 1.0], [3.0, 1.0], [1.0, 1.0]], + "expected_with_subtrees": + [[3.0, 1.0], [1.0, 1.0], [1.0, 1.0], + [1.0, 3.0], [1.0, 1.0], [3.0, 1.0], [1.0, 1.0]], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [[3.0, 1.0], [1.0, 1.0], [3.0, 5.0], [1.0, 1.0], [1.0, 3.0], [1.0, 1.0]], + "expected_with_subtrees": [[3.0, 1.0], [1.0, 1.0], [1.0, 1.0]], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [[1.0, 3.0], [1.0, 1.0]], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [[3.0, 1.0], [1.0, 1.0]], + "expected_with_subtrees": [[3.0, 1.0], [1.0, 1.0]], + }, + ], + "diameter_power_relations": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [2.0] * 8, + "expected_with_subtrees": [2.0] * 7, + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [2.0] * 6, + "expected_with_subtrees": [2.0] * 3, + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [2.0] * 2, + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [2.0] * 2, + "expected_with_subtrees": [2.0] * 2, + }, + ], + "bifurcation_partitions": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [3., 1., 1.666667, 1., 3., 1., 3., 1.], + "expected_with_subtrees": [3., 1., 1., 3., 1., 3., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [3., 1., 1.666667, 1., 3., 1. ], + "expected_with_subtrees": [3., 1., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [3., 1.], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [3., 1.], + "expected_with_subtrees": [3., 1.], + }, + ], + "section_path_distances": [ + { + # subtree path distances are calculated to the root of the subtree + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [ + 1.0, 2.0, 3.0, 3.0, 2.0, 1.414213, 3.414213, 4.414213, + 4.414213, 2.828427, 3.828427, 3.828427, 4.828427, 4.828427, + 1.0, 2.0, 3.0, 3.0, 2.0 + ], + "expected_with_subtrees": [ + 1., 2., 3., 3., 2., 1.414214, 3.414214, 4.414214, 4.414214, 1.414214, + 2.414214, 2.414214, 3.414214, 3.414214, 1., 2., 3., 3., 2. + ], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [ + 1., 2., 3., 3., 2., 1.414214, 3.414214, 4.414214, 4.414214, + 2.828427, 3.828427, 3.828427, 4.828427, 4.828427 + ], + "expected_with_subtrees": + [1., 2., 3., 3., 2., 1.414214, 3.414214, 4.414214, 4.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [1.414214, 2.414214, 2.414214, 3.414214, 3.414214], + }, + ], + "terminal_path_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [3., 3., 2., 4.414214, 4.414214, 3.828427, 4.828427, 4.828427, 3., 3., 2.], + "expected_with_subtrees": + [3., 3., 2., 4.414214, 4.414214, 2.414214, 3.414214, 3.414214, 3., 3., 2.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": + [3., 3., 2., 4.414214, 4.414214, 3.828427, 4.828427, 4.828427], + "expected_with_subtrees": [3., 3., 2., 4.414214, 4.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [2.414214, 3.414214, 3.414214], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [3., 3., 2.], + "expected_with_subtrees": [3., 3., 2.], + }, + ], + "principal_direction_extents": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": [2., 3.596771, 2.], + "expected_with_subtrees": [2., 3.154926, 2.235207, 2.], + }, + { + "kwargs": {"neurite_type": NeuriteType.basal_dendrite}, + "expected_wout_subtrees": [2., 3.596771], + "expected_with_subtrees": [2., 3.154926], + }, + { + "kwargs": {"neurite_type": NeuriteType.axon}, + "expected_wout_subtrees": [], + "expected_with_subtrees": [2.235207], + }, + { + "kwargs": {"neurite_type": NeuriteType.apical_dendrite}, + "expected_wout_subtrees": [2.], + "expected_with_subtrees": [2.], + }, + ], + "partition_asymmetry": [ + { + "kwargs": { + "neurite_type": NeuriteType.all, + "variant": "branch-order", + "method": "petilla", + }, + "expected_wout_subtrees": [0.5, 0.0, 0.25, 0.0, 0.5, 0.0, 0.5, 0.0], + "expected_with_subtrees": [0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.basal_dendrite, + "variant": "branch-order", + "method": "petilla", + }, + "expected_wout_subtrees": [0.5, 0.0, 0.25, 0.0, 0.5, 0.0], + "expected_with_subtrees": [0.5, 0.0, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.axon, + "variant": "branch-order", + "method": "petilla", + }, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.5, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.apical_dendrite, + "variant": "branch-order", + "method": "petilla", + }, + "expected_wout_subtrees": [0.5, 0.0], + "expected_with_subtrees": [0.5, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.all, + "variant": "length", + }, + "expected_wout_subtrees": [0.4, 0.0, 0.130601, 0.0, 0.184699, 0.0, 0.4, 0.0], + "expected_with_subtrees": [0.4, 0.0, 0.0, 0.369398, 0.0, 0.4, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.basal_dendrite, + "variant": "length", + }, + "expected_wout_subtrees": [0.4, 0.0, 0.130601, 0.0, 0.184699, 0.0], + "expected_with_subtrees": [0.4, 0.0, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.axon, + "variant": "length", + }, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.369398, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.apical_dendrite, + "variant": "length", + }, + "expected_wout_subtrees": [0.4, 0.0], + "expected_with_subtrees": [0.4, 0.0], + }, + ], + "partition_asymmetry_length": [ + { + "kwargs": { + "neurite_type": NeuriteType.all, + }, + "expected_wout_subtrees": [0.4, 0.0, 0.130601, 0.0, 0.184699, 0.0, 0.4, 0.0], + "expected_with_subtrees": [0.4, 0.0, 0.0, 0.369398, 0.0, 0.4, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.basal_dendrite, + }, + "expected_wout_subtrees": [0.4, 0.0, 0.130601, 0.0, 0.184699, 0.0], + "expected_with_subtrees": [0.4, 0.0, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.axon, + }, + "expected_wout_subtrees": [], + "expected_with_subtrees": [0.369398, 0.0], + }, + { + "kwargs": { + "neurite_type": NeuriteType.apical_dendrite, + }, + "expected_wout_subtrees": [0.4, 0.0], + "expected_with_subtrees": [0.4, 0.0], + }, + ], + "segment_path_lengths": [ + { + "kwargs": {"neurite_type": NeuriteType.all}, + "expected_wout_subtrees": + [1.0, 2.0, 3.0, 3.0, 2.0] + + [1.414213, 3.414213, 4.414213, 4.414213] + + [2.828427, 3.828427, 3.828427, 4.828427, 4.828427] + + [1.0, 2.0, 3.0, 3.0, 2.0], + "expected_with_subtrees": + [1.0, 2.0, 3.0, 3.0, 2.0] + + [1.414213, 3.414213, 4.414213, 4.414213] + + [1.414214, 2.414214, 2.414214, 3.414214, 3.414214] + + [1.0, 2.0, 3.0, 3.0, 2.0], + }, + ], + } + + features_not_tested = list( + set(_NEURITE_FEATURES) - set(features.keys()) - set(_MORPHOLOGY_FEATURES) + ) + + assert not features_not_tested, ( + "The following morphology tests need to be included in the tests:\n\n" + + "\n".join(sorted(features_not_tested)) + "\n" + ) + + return _dispatch_features(features, mode) + + +@pytest.mark.parametrize( + "feature_name, kwargs, expected", _neurite_features(mode="wout-subtrees") +) +def test_morphology__neurite_features_wout_subtrees(feature_name, kwargs, expected, mixed_morph): + _assert_feature_equal(mixed_morph, feature_name, expected, kwargs, use_subtrees=False) + + +@pytest.mark.parametrize( + "feature_name, kwargs, expected", _neurite_features(mode="with-subtrees") +) +def test_morphology__neurite_features_with_subtrees(feature_name, kwargs, expected, mixed_morph): + _assert_feature_equal(mixed_morph, feature_name, expected, kwargs, use_subtrees=True) +