From c5e8027d9c9d3b6fe52900c08b7d9dc5b92df92e Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Wed, 20 Dec 2023 10:53:24 -0800 Subject: [PATCH 1/6] updating plots to use subcell-analysis for compression metrics --- examples/actin/docker/src/actin.py | 31 +- examples/actin/template.xlsx | Bin 210581 -> 209775 bytes .../actin/actin_analyzer.py | 67 +- .../visualization/actin_visualization.py | 594 +++--------------- 4 files changed, 147 insertions(+), 545 deletions(-) diff --git a/examples/actin/docker/src/actin.py b/examples/actin/docker/src/actin.py index 97b7619..0571677 100644 --- a/examples/actin/docker/src/actin.py +++ b/examples/actin/docker/src/actin.py @@ -3,6 +3,7 @@ import os import argparse +import time import numpy as np import pandas @@ -116,6 +117,7 @@ def report_hardware_usage(): def analyze_results(parameters, save_pickle=False): # get analysis parameters + plot_actin_structure = parameters.get("plot_actin_structure", False) plot_actin_compression = parameters.get("plot_actin_compression", False) visualize_edges = parameters.get("visualize_edges", False) visualize_normals = parameters.get("visualize_normals", False) @@ -135,7 +137,7 @@ def analyze_results(parameters, save_pickle=False): fiber_chain_ids = None axis_positions = None new_chain_ids = None - if visualize_normals or visualize_control_pts or plot_actin_compression or visualize_edges: + if visualize_normals or visualize_control_pts or visualize_edges or plot_actin_structure or plot_actin_compression: periodic_boundary = parameters.get("periodic_boundary", False) post_processor = ReaddyPostProcessor( trajectory=ReaddyLoader( @@ -168,20 +170,27 @@ def analyze_results(parameters, save_pickle=False): ], polymer_number_range=5, ) - if visualize_normals or visualize_control_pts: - axis_positions, new_chain_ids = post_processor.linear_fiber_axis_positions( - fiber_chain_ids=fiber_chain_ids, - ideal_positions=ActinStructure.mother_positions[2:5], - ideal_vector_to_axis=ActinStructure.vector_to_axis(), - ) + axis_positions, new_chain_ids = post_processor.linear_fiber_axis_positions( + fiber_chain_ids=fiber_chain_ids, + ideal_positions=ActinStructure.mother_positions[2:5], + ideal_vector_to_axis=ActinStructure.vector_to_axis(), + ) # create plots + if plot_actin_structure: + print("plot actin structure metrics") + traj_data.plots = ActinVisualization.generate_filament_structure_plots( + post_processor.trajectory, + parameters["box_size"], + periodic_boundary=True, + plots=traj_data.plots, + ) + if plot_actin_compression: - print("plot actin compression") + print("plot actin compression metrics") traj_data.plots = ActinVisualization.generate_actin_compression_plots( - post_processor, - fiber_chain_ids, - temperature_c=parameters["temperature_C"], + axis_positions, + plots=traj_data.plots, ) # add annotation objects to the spatial data diff --git a/examples/actin/template.xlsx b/examples/actin/template.xlsx index 7cfcace7e7f33695a69ef37b2f54ff15f3f084bc..c215d429d324253d930a5f5b5d2fa85c62b3640f 100644 GIT binary patch delta 23688 zcmdpec|25o`+w(@6h)}WzGRmnWoMM7WJ@Y7GA(y$Axqhp(h8Xu`}uyK=exXKzu)hlaL)C)miP6(uIqEAGoPf9CDi^U z`Rf)7I9t|s7{VrtWg%DJRJ~9yu8 z-9?h?o!TYux4TFxVABZE5qpHjzs|FIC{BD|5-YYX)L#2eqmA_mXxHw=yTz7TGcRpd z(bU{-pZmyIC<422#et`#8!q)}N*#Q9(8S5MbLW$fGRy2u-dW+29=KZUNy3BfV$Ao3 zp1SuJDkfhs;2hkrFXv5fUCo}&kH2;W#l$$PcS}|(u--1<*yo=^Tw_LgEMwA z&6~J~6b?izJ2z^)Fh?nRGt}&LB;eU0)w{uuF9+>GwyJ&Xu zY`A~L!L&}FtfA}K%#B*G=BR0!M!Ck#Y}=k|`!zn&_irp}cyL&Dr@qJ;$&ydCOuHBJ z0n_XYMf{8i0#k@DjM z<)M&_oFb;z^PJHh)3K-&&kfQQ30|9pvMN4xes++Uth7$Alftso8tvZsB|4Y6uDcZ$ zBAwlNTI_=7t%Y~}xuej*F{IJ1p{q3yS2+=~QE>^XE^L1~dunQyblHEE0SQuNkg$ zR`-Uv!Qu&{C*q7V&Lx>}Wp7IV8Ido&hlDl89+$iAeb<&{ooV>=!@^I2zjy{fTnnM;JPVcbUFn@7%fFUp}+>lq9f;cq@OxRe*X^P2wS zQq6a3JWNit4I3N!6s31KJu;MN>ymFbNS;1%+~e~oH?`Yn`a(e19si)6Me|+c0cW1nBe#{n_#M9%*$XmVnD++Hdc`7KRQ$LtGZ={hZnVVJvL2SQ^=H04k}<8OCttzsxh@B;Sa_s&Zy$Zb}{8*iiQL$Rl;#h6wh% z&WLlqo=XYt_&k4P>snAz`NM%zpG5uc_g8h?FE$_WxJqjEaPhbw8b{6flbqUjIXmn= z@<3rf>fm}jTQGMh!NJhia;C|ZTqf;Fu=IlIQtC%XJ{8(*lFFD`5oj^~;)z5= z@1BtIDUU1G!`F?oL{h9rw;JCNNwOY!Yy3nc-I}Urd_@GV8qPOsV!So21;avz+jp`;MFY6*4w`@nGap7-hs zA@YyBGCZS)V;*;r&1aN#gsT0;R@HSTP&PqNn9#gb$0(NsoHq(6-7rd1QxGgnLZ zI>K93l@+t=%0hB(J#y&FR8xrNtEqT9Hb8%v*G((Og-;Z$gEz>oR&qfG_?q17tzAN( zl8jSgqL`y1rH9GB2gm)^NY?6rAU&7 zDo~4KfrhZZ&T0iWC*M({D_Dy++$D`sEt51_ih~okdkzFjlUkz6YT4G)q{Hl`h`%=( zj~6Y!+2Sw7Ae2zdNc|Uf%fh14DrzPjzAK@s-R$o5Q6FaNGOsVA=OV;#?#<@Rh1R@zm7+xjWXU`r)EU=w^> zNe>ijg8?uxKP0YhXGm|kSUf%~Ssisi~j@Q%@!i+v3W5CqGa0DLL z)D#-;$)4dT=i@7azu(P4@9rJ=d3PuOZX9}-q@p3>CZ^!oZcEzxvc)H7>2mm*islk! z&vq*kZoyQes4jqrnS)9&PJW@VZclYNjb*#PjotM_$`F+bs{gB$9ABz%PHK_z$u=GA z{=y6JG9^vHOKpaAqMvkNez0BMRCu33fX5G#QBPUU7{hJdnumi zlk7rF|5ftRuaXCUQ6B%Ra+Uv)B^tsx%9b5PyDd9{zWn03#IMRP{OZqvU;J_T)t~?R zlI52}LL%FaGW#l8j@t`uR?<|i&bol$)&Wv0vZ9|P8a78M$MqFi(7M+zbV@up?goz= zn8C}|FA!{d&m5E9s}66)SP8v8SnMMZ8FF>i4WwWTh!qz3t)SiTx@Fm`jav{uXY_*L z<6S4)oE1C|nXAH1YpfSG+!rMjZqo8&4cWnJD5OK)=++-z}7@oHA=I(TLCmM@!Gu{z0h9t?C>F~T5~g-z3B%D=F8 zK$*|aM>TV&$ZN;Oy6VVYvnCE;lm$$`|F(Pq8iX!9s_y{w!_c}dKs_vKv;|0k@g`g5 zuWju>e71u&e8+*q@R-r&`KS^2jHw>Hf6YpPc#>JPY4ZDX|)uMkp`sNyft|EpHE0Dtf1CehocbmEP4pSdv7B|Wx|&4W*@*I0rS zk6v_1-h|^?V$UE}_8@ivnvj2V{hulQo7i2aC=gtXk-F$hozNXT>tSdKQp>6i} zH?PlCdYX0ZJ-o1a^R3MW$CMssU44Q_$x|O_6!ExhvX4rezZrJLU>ESO{rm9F`M6e- z8!us2^Jsp}V-3@L;sz?pE%$5++^wH}WeaW-)DncSd0nHEUFib7TM)5Npxc6yD}4mm z$;wUusUuI865l-qH-bPUe9YEvK@hVvst})O1&XXqhD0RWu~c2j_6<}re0_`A0<137 zN)%&(8$je+Ft`+Uw=xG53ePs!fdCab;0f9Z5Pm4cg3BZeJUC6C&#@q??#gz1syb5X z0~#&n@tOiAXBA~~P(@&9FkHUf7NmL(x2imd_djeeH zGNV7>)7(Uygup`}y#JUHG5#F56NGJlGeTse!G{8MSqI1O!wKi5i2c#vRbW=+ehgn^ z;G7ZB=rZUI!r>>4iLb8kuPydBBD%$+zx`Xr-zV^m^R@Wi;A`1`Y#ouA4DOk;aXE!= zqcPH$xGD|Y0m4Di#>DIlzMRJyBg7yR91|Ftb#jBF6k>G;bcgkJb;HL_imdOrooY{n zN7Xu9!fw)eX`f7zSG;*@`%c#$_cc6Gc=h!aJCfh>F0xm=jx~$F=SMN_ekk^;0-nrU z0|xXJp1Fz=%@p2|XuE*q!R;i*?*^?wq$vOt0!69tw>&d|2os`5iB?{qI{+|YOF$vG zF54Uor4apY@cko+&0oU5LKGCAZvY5nsIQ1~hr3;+z>*9EHGzS*Zg@FYk*H#nxU;KZ zEO`lN=>ZpAHY0W?fSW*o0ZZd;pK8wE3<4(b)ihfGDjdmN3lQ~gg5tokOML_=ZOt=Y zWN<;3^NijAd?3-#h|9pEm-(@5*YV_F)r0br1QPGL(vwWJ3-gZ)ox6Nn`BD!HyBVLa zK29|oyg3?pC}`)BibDTdBcSkkqac7>A48)NuCgG;`hj-f64K9xVB`*uzBdDKM6^=v zK$sLcD+UA+pND~(0MZo&euB4-+AP7w2#?}}5#iC>sd=z%gBejj1Pumu$z2Rn%Pm0O zUl#x$;QD8_fCusV3osB6Ai$Ta^I-q&QecM@oY}kv-fwe$f$kx+XHd$&_Y4~l83+c2 z;gUiNzVfXgUoo*F4AlPV#>Qo0b)`I$y?#9EY^skm!=IiR3NU!Q-JOM&o9)1*TvHD} zURkM!pI}Fs46OOuNI;IQ0~Hr`x7vbm-kVJVDBm8D3&%5cX91T7VCx18;*C7i-29SI zBY0_*8BzHWXag=lyQTuXv)K+Pg}YnLH0dOi<2&WH6rmh1XYam5QC_zB_l2={cuSeu zH=Nr;A1n#JCebNu@T;)BT%s>_uTq zmzlK#Ls`CaU4oKTWA6e-In7uMgf`b@wc7hi+UT5kWH!AUrFg!rH3&bXE!8A1ElD}e z=Xa?+1GaRYS%?k<8a6Q0`3IRS_-(hHAd@wV4E%1(Nni$E&ej3k ziDF!^836haz4?IMB7bIh;QHTpMHKpvTY?}-9;pxjY!?x9A;P1DE8+^TlP$pcSNe#E z!y;6;p^y6E!|j9PM61gn76hV@&q9DLNGc^N3Zp+=uxP8m^Z0uQ0%74kmy#6vk>-&s zBHZZ|hDEs;4_i9TjzU@%29{T|a7O3fv1Xz@S>3`VEQQ{C?$Y6gKtZHW6d}30pTAsP zsGo!dT`1k>STu6qIkRt*A8rrX8^N89CF%5_=2CQb5qWdAMuOo{^~enR;pQr|IB@C~ zPBVrBtSpK>29EMrR`8&BwCgc@V1hAtd~3lpQ7H|y1OYDcNeHl8%sc%@_@RIVSD_Be zW#LKO8P*IL;<0eN`NR#Pz(PPs?P9&L{Yqfo^FaPbe%jEcIO5BTeTuPS>#RiWngzaH z&9*ZX)77#MJ@yu&UQaxn)2MXmfE(SfzaiDg+KCXn%|=Eq`9u<|<7~g+^1!*g9nKX{ zDe0MVLg4ll+m8yEd_J$*TG)io-s_VwBK_R<(rs_&kBvLv)dk*56Tj?ozde@q-fmkA z)Sy+61njw>2x%sI-nu4kU#z_J#mh25Yedycuk}IlZIAf}t+vLnM{c=Kol*N7-nDvt zQ($1_lLyNV)(Ql(=50}~RTI^|w{)M-gK%W2MqJQBk_CHG^w>F$4*a$LXr64_xiA=cWt-t}e%Pjbu3#l(CwCL6b3E?(^{c86R&u(Y zohXk~ea5X#y=b!g=I8CL&VB;I5})@RtePS!(hb#(??3ELAbf0;N&LKSphxJMO3&fg zB{xipYt=}vdJ`Cnb+N|YBej+RrLT&>5=rjv zsJZy^^INxrd-kjNUtP_>=$DnAZhew(PgbaL=Xic;3X>Z>yGw>_bH`apye67CF*k7y zLh$$d<)&lJ)%8jNi^5_oWruo96XOpByngh zI4%>x@&>Rx4v#yA=S}oB#ADfP1&CSGG(*64R&d7A&qLgJliV5wE~}r5fjGS+F3AOA zF$oYWW~_ng3$bZ!<(rwr}S&gICb&V9%-GMBQIGP{f@HGLl%+dk<%V;2Ch-I z8HR+5^!lG~bI=rV5ebi}@9u)xQmgb%EU?!MCO^E* zWllgzpUH}tnx%KP%l_)p&P9paM@f6NYRo=c_cf$fKTdbmsUcYNXvLTv3Qe0DE7ML! zh^U?W+ImM&i4^8JdB(g`lk>EUh^q}!{*4-6LlbNwU2Iv}L|Vw#{q0ZTYAs6KN`EtG z79CLLCjCZG_^{ddMh%K+hP6q@nC#`M-i59jzXqU%2A~lZOHbcu(pihYw)0KU*J#vr zb0naow|F^U{kA{3zIIXKHoDW7V6DCobxbwxw2xOJ;jYnCM?3Si$<3ClO_j#xCzlBn zR(3FLSup;S*>56($fFg&O5M)~{CBBrKF1vEIEPv~S6oc=stfARt4xbeJ8#l4EqjIU zn%=Lj>7uUb)tsfIY^pbcXuLkl&N$Lm(U^*vEu4UUk9uY&L}YR*RPd9Az0xX@sz|*0|xl7`f0HeZH*pZ7kf=zUlMQz?pBvHJE3#zV9nt{mU_|1 zKKp~*H{YgTHoI539Sl}z=8Z&ok@1n{IKR+K>6Q)m^pG0A9UeCcr(BE+V}83nAKsuQ zWpc<~*i2*(nUwFN<{(ssyc((B`Duh2KNLd8TD;V$jW@~>){XvUb$r+_tM$>VYl-EN@)D8DY+8_dt%S#zwRb4au z%q5EF^{FUKttnp+QYuO}v8WmyUs5P+IVdVh#N;6_Jc7C%>0VLj@`$^UFH(L$X(a8i z=ac^@HaFVm)UjA43_?qK5^vUx9)XSZ$!D-J2L1=5j>H?R`Z^Ce zs|_qi&X|A;^mjGCye2WoaPgKo0>5O7hQ|Arn?eRmZ`Yjp$Q{&YHGfx2aY)3g=MKAG z-fQh^04Xp4r7#6qwsJ-vskrw=?P|76{6ey~S4DSfO^d*6pxOQhDzx~EUi`XfL&56k z4omDuD73J1l58GToL4SaM;=_&=p9|oDIdWxO36$I7EOTiLa|hFyp4H{^c-U9haFRB zw29vElOg&jLkLs=y4POxz}-_i4~H$8d$#eqR#d8$xrg#SB?9c^pIyRDKwVpxtNcM zDmS4VqV0wx8(&5YSHHkt-o`xhUEDT zdRnuI?bk2k8FSzzmfW-LQ;?eO5|t=)!4ehAAAHe7`4XNa#!FV=K8&C4t~V+Zsk#VN zRNvF{p6+%m+c|xZcflZh^N(P3e+2V+%(u9kJabW@F`{Zn&wGZl*F&hS@nsYGYxn)> zRTrHUeiB2oAS=53O*T~L#QSK(nRBS+r$ec;4yFEdNav8fP-%FQ6R&PxixZhFz;QFC z4dG(aF~$&P@n?ZpAz|^zqzmSMBwuqb`2%GC781<~i>97X)yS{TYy}k911+Lncf6#2 z6RlMgg#DGJbY3N{oYG3=${BIY*lCxhuQEyfO{!H+cS$n| zM6W_=g@)_sA!Qhu_Fcc$FZ#)-gzvOgtWiIjv8a0rt96&8Mj1;;{tC{z@8Hxp@fuNZ zj%W{$G8hAG!JSxI_D9+ZBu1(;1&t&2z#+h#d3paYv|T$(+YhAP z`w*3SSH$ZY1aPYHJ3>rOg`V^wemsmkaK?`$uE=Kl4C&YX`YIve}qm zryO>;!k(&CljHSG$G9QPZzdYeqwLVe4Fa(h;u|v*|DdGa91g#`?;fQVx06|qq#MhY zzYmJ}=1(68sAMh<QFI|{-Z5j^wykY)iR#Ob2G3fkLfgyMo83Pe zg$6ng-{;kjG19uJk*k+_c?umbwHE&_oL|Bh-qrC#*n}@!m{-ch8P5n8HnFWWX1djL zpv8R|!mRQpe0ROWxOed89$*BWlj5RNbU(Fgg5jn=Y8V|+y`e!ayzG!?>10-ZC@Gye zsNsiV9euiDcp{NS|E#;Xd*4{Mv;W2rN%J48&_a_43&QZSFFd0VLy$iiRl2&uU$re zHl?dGV!VjMzOj@#x#o_U%HSct{%G~4zPl84(PJp&lveW%cS3&If#{+jHWijm8dXJI>bMiMzH_v?ek*_TU8z+m#_|suet?Bm+f|AYp+VFm zRUyJh#LF}6bNy3hK&5I%3#GWs`_&g&-MHS$FIWCNuu)!yMA5#mpFd?*w)Xq=-wV;x zf=WF6Vg=1k880V!zhXW1Kk`chd__wQ?iU=b7o8ToOe!&mF?$I$+I~sLHucr^=1uT- z6cf{b>?rscF&oTi^U*crl&TaY!@NBOBk7T09%2+|yc`=Eu_Nk_y^9(@3BG90dS*bH zOgBjj)s!EtM}&39_j`%N^n~Oa{Lxy_eiTa?L1jtpmAqqzc-Jgsg(jzQ*V!S=F=lEY zTJ9;0xM$~Vm{z}*`}(|K|Axfs|HVdt>Cw4$mXv_9bBy%2!wIFiE?-GhK1q?NxBUdfI2} z9~&F(f~?xICtAFXYts^aJEw!U`$E;p=zqg6j@cVvydviO8gz2w*MCa$4UZ+A(haYj<5olfFvl^^c^Z72?!?PkbzG&=&A~K zaZ2{PA`FQ9Dv2agc_w|n$P`BC{dtAdH z=_CQsl^#(nWb7N2L>AN;;w&-ogaGc9~0^Ks4bRkQ7A#4W?iq!%$MjT z3z8zHo0r7&(7a04ozFB*izzfH(+7crOo z0}5)pnl&EnbtR|od|%fp;Eh*6vl99AdrZ1r6r--6qg)bgEDtwh_eY}9NbZC{9kR^6$Gx(lQT>9k zcNHU*NBr_?6kV5JSAYFu$Dzu%Sdu%=Id;6v=4{f~A+#GUQ=FCaUp}3=@DGj3hwSe@ zJLwP+zea!aMAma~IyNWh*@r{djfKh{F8||EPzBwLEpK&-rRu9Lw7>J5WsMzwFG37k zBEG5HtWV;v;VAO4X_F)$!YxU(sZ`K@K8LoffSTPe<*iVV&YIW?nyDBc< z{gaaWqiHI>O0@IuSi+Py`0e^?p~mY!QK^~t2fm84(XN`XIPzN2etE=QwI7;pD*d5J zc~+CmdNOMX(+$;Cu^9dlsm`rG@KBh=<6T6&t-)q-0ZFHfFWoy)HfiJ^?fIEcCA5F` zR0YQMBXxZ|t$p{0H{t; z^g68g(Wdk-~0<1 z8Y9trzZp>8$B~`Or>B`If0k6@OZsFP`RXO|B@G6CNb27A)JsZyABRUO4FR25v4$K} zA)YTLJorQGxBH&?cm6E4CJ)v2kT2HzQ|uP1E&r#OCg14QM>C;JEO)Fh61s48vUgzD z`|B0h7bcTEL^q7(llzk1gw>Nh$hZ-pqwHYwu`V$DLDVkG@NQfnR?Z>3a!#(I^1fvK z;~#P=vvO_^qqbRwSFX*Y2)BOmaxCwnzQ=2P#^rySxDht)hP*RDG0I;LEL2t=>ML2{ zjt#$AzW>P%jBEg2AGO^=t48)(;Pm=`Itz*2D{uEGrbZR(Msmb8LdzXkS@@l7)u-Ru zM}K?g5@+Gcj#%vO`7JhQD(lK=#lEa_+Y`06IUV1;CRb%L>+HHjwPwdXn@?=~Sm$5M ziYKwH3HyAxA$f=i#LMA2kfR`u*Xz!9Z{huPHyP&}Q;+2wJV_TETF}JC=d`JoPlMvQ zh${xr<|kG_DfV@kiEd~s{qdSJ_86Ia<Te&P!b|UQa z%I8_g)n+u6Xg<+?Q|+}0*`S9q@~tkgNai0H`FG^f5O@7j*O9~;MTVLqjkY!mkF+fH z34R17_2_#vJDfp1?koU zG!_2ZD_Uqjs8iGQyeKT5$Ih8HxAk30JwmX=*xD`0BfYn3RsM|Tm(#rir;{5fS7urMYuOFo=mHG zLM;^~u^~kGIA92zKqjMsWHDbJjl`q!CLnGN#KF06rHKa?fLj3Mb`fwHDH{gF5Na+U zEYy+UY033t@s=1M`L6-hBwtWi6n;8yl#Ilj4%^*S^s4SCBDN5;LGm{Mi-;NsAORwH zH(&>`J{j;u-^FW3GTi_-poZ9x44ec3GD3R-YzK;nzGu*&=F=w_NOvxJfk1S=27mz3 zst-5@EHEu+_aZ|e=tlfV1$@8-&KO=9qW2ipLz&Y9LY;GZXaDvFJ*w(RyteYlV9w~L zV1<^^KHQTW-W`inN(MIDhmbnZ=SVmbX3>(G(%DFdOo~s*IwcRc-*HniK3j9YYlWm)lWHwaluj^LD}4 zPYoB0utuv2Bgt|D)xDxbUXxx+Ua~8y;p$glu|(jSGb~1*Qo85I!iVE)r7;!^3emG2 z*b8dbHbgQQ)W%4LlAPy}@R{F=yyAVx!3r58;T@YX^3XGCY43+T=HR3wg0ljY)l9dQ zc{p@lrPlqD=WCI{6}|T;+z7`QL^q*qU3nUEAOWxmsa zUD)_^6&fDFv%#ILy06dx7n1&BU9N7= zJO;JS+O+(ES51W(dA#`)WSL#%fNb`><=KdcZ3nDSCi(dQVxYv|z>K851~$sI{=xG9 zEP7rLKmUkGAR2v13t%kgdje=DVMGV$k2B<0qY*U@O&sP;5zId2&|_gH{H$9cAIkv? z;@f+GJ(}=Cz~x2ho?K3U1&f7HK~Qq3WIVQ=(J?;wA!0Ht8?icrrbSNE^N}j`6{C;b z3O<8y^34N>EzG7mnAeaD7ENEyv&3einI3@X+W}UfC7$TY0z^TDGWNL3FI|JwlEF=% z^5!q_3GZ;kAp0i(Z55U3x~E?`WD*8qx9pimq}Cnc1Bd1$xD!S0Z&?V}x4vb7Xg z41OkIK3e6defPwl8q9+P$OJd!GHZPmD;)@rAcY znSraQnaWf4Igd)kc$tNM8=@8xh|d>+s_46`_J9)-U;&5&jR^KJU<0(xHD56@TES4% z&k80r64Bc)&{xFfZ5h_I)$^3KAj(`C9NY6d= z#md8gwBQr3_>pObD>7~kSc9}kMClw*g3d!Ce;!KA=Q0pPArP=6F1i6kfyg5_z>esj zgHBZRok6)AU^~&(nt#2*6R-nXiC;FO)>m_3Z*q;LfM!)1Y^JyxQ(T#e2jovdJ)Ar{Z&O(QJ6FHqVYHeCpKlPP?Zaa!1$&n#N)IgQ;5Oy~Tjys`04-q>KtVVo41B>NI z1OY#E`;ygC7Zx;JN|Qr(2!Dp>$J{_ab*N}~>3~BrF9R}*`9eI)LGfCo;2^LOjKm?8 zmw{6n{E5$F|17mA{K&weHnh=tvB;bL)&U{e2mhjdF;Knj-;ECBeaCenb(Ss&sS|Jk zSVZz`bf3ARs!o9y=1bm3LI-BelW->X3(UKVoc;({$WqZ;C4Z^}1$1v=;%RDTiAT^O zYWC-S^*_=Zs{c3WeRsaXf5vYtx<&h6;5Rz^{#*S1m)HLje*Y)y|KK-*^4m+F38l50 zTs!QP?D;NtS=zk961EDED0G*~a_(DV44oaiSbVpB%KqM?k^|IdkNu)u#)#ebf#-9Z zXGG^QAV##I07ujjtzCe%MGdd%P5JPipd>G;iWOgl&q(Rzt5{wBomLqbAQ%3+V57AC zYQLO4_oJAkugx~=;9S-*g{LGNFS(n#CV0RZy8GItfzN4`g(;KwXjkw z>*@a4E_*K{hiXg5I1`h5<~5Nq#AA%d*L+x`<~0t4-g||;Rp{u~^t=^vtmj@VqBi>r zLXZMAk#sT8MvV6m#~x>pdBX%~deFrP7f&QBgXnWZ?sCu;Ldo4%9m$|k&<<(n<9Vbs z*r1J58q^u=8=ubTNRaP{9HaB_2eQ=-EYb%)J2rXE1Z7Jal&1}Ra3pwz8MFi&l&24L zI8wc4^jd0CEPXURuc=aBFzPU^UMG*N=+gC>+jNGgd4^q{^Akyt z1=l32gNqi<1%YHrfrZ2dIdD0Eyj4d1YaZ=$mt=G_;n7~v_9QfAq_e@BLH!a+F1HBk zL>~Z7Gb(8nX3iEVIT?*~?`(=eq*+TGe=~4AyVamPWuTP*2zN4@VGwz9pw|&!v>a2D z7~IY3b3Xx?L!@?e@R=5kl+t7|<)fARzLcGg{d#X;+_8+l9OG2BZoC<3lLcjwW%8hy zXawoXt6$*$J5O8xH#}|q-|#dwjx|bYLLU!(zLeDmQI8a{f=9}zo$_YrG_WXPb%m1U zax&`O{ylOw`~x||{Sd>srxHyP`+(2#M zI5K10lC!3|uHfsw=_uy^xR)dU8(xn5Ph;8fzu{%af7?q0&A0&M8^Y|T)Zn_vJO(+6 zP2h#ltI*G9a8c|k_t%`JRuY%S;GvhpsWGHTUKMnlN2d`;%)xjPkLv|-s@!=2w9=S^(|?N ze~pqXDGx!mP_E=7(rcpRqbDT4br^GE+1PaHh(oy_XX&$Yzk8hMyP@61ebF^IA>3f} z5aHTLwDt%_wn4aQ4aS`7Gjq@L88^N3r{D;tj9Ck*;8@Nx5v^Qmf7_FCzoKxXy>-Iz zWHN?M9y`KG#`9!4vpLa2pQbQRp%D+ykceKs54ia6s-db=b;3;~^NRIY-YXBbLMV5v zo&8d@$y+|e5)sUR4s$C78x>Mg@zV(@4fC5Q6VGXwH(2ykl=A^J39Lq8Xcv!nKHwC0 z15;4J7;>>BkCIg(p72)}BlH&Af>}-GPE?6vJ6W`NXaTpE2Ovu{K}4u@g7dL!{=9kq z%aDhfpe(eH&|TdIWqj! z(vT8jUNaRjIX&KNSc>Z^TRUSC&g#EDVkMcqy(FwLDuTzDjPa(p4sJ8Ft2z}VfWA>( zJ~Q;qV)W?4sixgmkMcfKBVDbvIWmVfNCX%@O4^Kt_B~~MJgXDA^o&m-WGdT9-J5&5 zNZTvZzmcBrSdcE$bcb~JjA-a%hG|9yC_5ox z?Em?<{3j-|HR)H$$7Gat6Gq0Sp5X^W4tk9*t7Tm<7o3# ziAEUW;L3r_Jc|nj8VxhHOVvcY1+Gy6Zp>F;`=w=fP)1}k0H#M|0GeBfdg zeLBg0|M}Vq^rh3g*NsMlpKQlim5sbQgXirN6_jUJG z6(UQevMzxP6_31AI(-R-Rw-@icuw6>Ym(GP+{T>l<rC-G@`_ z5t2TruKeL!sTlUH2YuflDz`p(%99q*v%R6CJhFoMvAV0Tr+5{zOdnK}DLq0jW^8Fl z8-3#2e||kG^uB(61vZj)d@KODwiYxH=jQhFErgebO#ix25nZT0q+Q-JyAp=N!%QaV znXf9$B766Y1&pL!K@!%2>I!S|gyi}G=j1C;Fe-4PKAxML zCwH(v)bstW)ABq^W==XSC5Nc=03uh&reR1cCnru+tn3r--T|MLl~>>USu${sFk$=q zy0K@PgD^He#DqrQ6u1`>F(Q{prx7ZqzP5Q?GC&V*QcwxhW@R{Q_auddty*nl(aq|p zOL22q+1A}JBYB=y-FhoRY$-x80awek9NSm6%?V@Omz^!hkBoN2iTq5eeOfo{0_>-* zq*0O->6FtkzGZ#66+8{!Ik|Gba_LpaT~55q)nHeS@&nMuotSVne3^nCi==RL*272| z-@|*Ux7Z1_4YZUN(G1Q|O96VE8xhB2kH(NVGXyAt5kh0;RJFuHElt)Htd`ym7H={n zhxiPN%&5Z-3Vtx;6`OEHETa0RNERJii$!N|8RO4QD9tO7*IYj2H; zrKDPZV2QS?8S_qh4pp}Zq#>17pt4zOczKR6F`~qe>LQ7rMzyD$c%N}v@dW@%2skU|)!?V8i;-k%8mMYI&|mBQ>__j@v9i&aMMhu=_lzv9H(QFSWBI z5-?sBANvFYLrdaH-jlV$qBfPcMaYNndUN)^WU}ijn7T9L%$6Y;jBq1vc zHfZwFTf<)j)T!$*xEK2BQV1WIu7Rizrr7fL%mch*lU0-v=*jX=%f=(l=QIhklzN!i z82nCBaB3I3H1wgQ%V|F(-vm@aUKoSRwYHV#l>XkyI#rIVppbau{ZZsJ4(S!mfz&+0 zsVf~0J}I6CNfJKH6lTp9q0t9VYVErF7P04329xZ|L0>uU@aZ%c zPm3>0o(v#MB9Blx^o)~AEU2lvPc^tTd@cQS9tX}8gf?>f>S@?){~70_(28NR*7G!i ze-z~`zvi@2P~2L>Jx|q8){An}(D&CV*`rbIP+1rO+ceXHKvM#h%)F3?B2ByJVpr{$ z%qVOhH>o?1?d1ma-!o-xC#ZikRFUU?mRvv~^gJ~+SE-yWmer{0z}9fPfEusWqSKPm z8BONyjE@RnWGO441)84Ii||ZAZN>KS#jnp}-B5jf<&8UutX^tF<+K-WSsLkb>)_aL z^Jh7*`olhI?JGIFxA`=!3*`WNclXt%zzPPp5VNQeSA8NdM0t*Et5LEIK4RGPj9yh2 z=RRRX^bb&PvcK`N$=&?eCtI_}%FwU4p|o1;w2vBEneICnkDT2EYAvpvEMl7R+GOYc z6s(;Lnlc+6=HJs?>RuZ&G&1Z${K8GOnY)j?H3t_g5Utoc3%ZWr-V*He9rhx`Y9olA zn&#n~!tD|sQ$Z{&K-~rBF-3EW7S&EBlSM1?E%}K;*Br-AADNk3gE(1%`pBducu+X+ zYwu)HYtG)0BJ`Xbb5IA7+yEXB81f+~6d=#BpdN4#8ApH2e>o*IYi=QGE}U}-d1DFc zAw(qD504GpHK)#Z~X^3Dd-5=cbv;SmxR)CUgCX?pQplWtO3arxlB zZj!Ghn_cGuaq(TPEMX~kw!5HbpA|K`=pqRyLxnCG5u2%!XZO)zM9dn6 zSG$MXu#aA9YhLb6Jj=&MYq6N`$=m2TY~d7SnKh`hxMA)vwX$G-#M+CCeWFNHB5B&A(j@P2C~iuJg|7m)}v^5=b^$2vu+>?cA&s~YQgSV{Q_c$Pcj2( zvO$$Sx)f7qSeJv05cO@rtpGA(gC450Vx0N7VPfELDPn2^VnFtl|JKh1>?;l8$PyGT zk;~6Qq0dK|?2?Y-_NLd7Ee@c@qPnU2REE^LpNEBs9=2Xtci4^65bmA+ZLCjEZ~ zJ7jG8dOxNZCbpt#;MoperaOBx8S<7Nv_8F7zrNobb{5!nIUgo|rs+F^n7kzy*6Ngt z15fl?e_MaYg2TKT3bLL^XNm~|&lG!iBFXY>sra7FcYc2QUa>vjk7+Iw>%&99Gsvb3 z1~Qql@$dU%dk=i_m(#RRCfn%^flRZtru{$PIXC%wV{)X(+@Q{xf)if`ExNIKVN}sS z=9SD-oA=KTE%;vb;apSU((N0&m^>M`ZwX^ERh_)RhE?PLX^_r^(a|DrFI6x5le=6g zbmFRr;8*p||C;7JI43OmQilDe6vwvBianx?+qG+%a#SZP*s^GU0-b+UCTlUx^`d5K zSeK-rp6V*am2GofvM}tNoVfJ%^nHTNyxZ3;W}2otIl-1y<7d3z8`;%-PrDXnOuu*j zgr6f*t)K$ugkSds)=!$e>$A(clD^=ZQldGNKVDn9UF;Oo9o5MOEv%YQv&#(QTULNeZzf`$=*DCAVrOnFI z4a=EKwzG*d-WuYwy$>6R9KFx_zkSkoGhgmkPPY$EvCUPmDEuKd zX}8{unh7qu|9wxpyXNy_!#TzE-|tM$stbty#CXnoakk4Q8Px?l1%tMmR^8LNY(3@2 zfh`GhF9s}Tw3%goJzRAA6lwPUa}foOP5w%r)u(cMzbQs^ZG2uGa3=lY!o@p`{Wc2B zUkBWj^U?eolisI4Y}yy8Y(*KUMqS@rIDO?vr**o+ni2si^<-@8Pu0SxWk+rif_H%;j=A zsjHi|Ui2St#xYa-`jtOfrt2Fs-+-LZJ)OsdIf$_ncswDq17qa$1t!c2GKeFGIT$2> z=R+WXoyg=^5vl3lOqd0v5vP=(tJ0sYXv(ZDgBU1cU=Tr=0c5RHo1S6HtR@5QRsLgO zIGuoQ!~w18t4x_)m>z3R|6s~&!?a#!x}F)ck_;%8)fpKWJ^`^HikTaAr{|b4E3@$% I0LP{n08%@D1ONa4 delta 23196 zcmdqJc|27A`#0_!Bq<6Z+mN-yBudsPg;1ouEYl_}_NB3pw0O6nPB}?PL{dhSu@6Zq zF(Fizu}{d1eT>=eGg{xD_viEbe((Et|8d{vdhv>v=t|>-CzIKDhYy=;Ay> zQQ?Z+W!}b+v8W>6)u&HeIs$u&cfO|ds^7Qp13@;-{K8$Pk58D49<%b?m438OxT8qE z9#0Xw+?;XnX}*#UAP83%Gn8}Dq`=!_Zi!7Ut*chb0uu+4>!RDptLFStnaBE1Z5fGP zbg)3WT>pB73x`=dVEq`reG9ED@X^HTq5h8OXASoIoFY^e@cXLY9eR94Uhc8S6NQ?S zY4)$TYxv3uZEhbb+wb_ML)LL2ng0)AurseV6Y=d-YNC zPYWEa5jOG>mk!yR;(}@0)p%QgXet?&98FouRgcpA<{14DN?)C#lRfaoMZmD@rp0)w z(4F}7rGZ)c8zolZUtTHqPpUT8S@ti_>&QzNZP`$eAa+1w^Y)b&2s8c_yVpi_I9@#R zIP+4yyVs`Ero$v1ofvbp#$TUUuesjEM}2XdEW5>@8#;cqqd4Jrhdz~D^%7$^8?!&! zU&9w2c3t+{@|Byf_SO@QS6z7hBz>>io!5p&Wv@kgGqQ|&@Nb~O4JL9c^WDoMcydZA zgqv>BIO~wH$nv6dhT&O+c z=v5qiP{z@wDzM}l9)Hf@^qu~jg_#`#?Uq^E(fuiFI-7E_s+*f8dXrz14>M*9c^{_E zke*Z*=9~K#HqzhaZ6Sp92dNr0q!YQx7moS*_U%&+Ir2_=Nw44YRMypTs^nM8w*`<2 z-d@MbCQ!!NWQAFUSM1krL5tid5ZxB`LaN-}MEbO|vuo6iEJ=LWg%Z8Lrl}u3xl^9! z$n^W?i7fgO={W0079%7Uz6iy?>aGm$Hf)lr%YH+@QyMmpt$0LXp!F0+5QRbOF5hg^1|M*QQ%vvIGQ(1$(7Wf_RT;_n!UBiqP&G zJv48DS2mC*i*{#Uc;^GYNZc*Fnei40UO)(|j)3R3@s{9O zoHhzqGkkV@I!8YhKQ|pp#$bIy%l2|hr$~{JvqO%&=}AZA@sZ19v+*H+3&}Fy{%Pj6 z*}#5BEDG|H)x}KPXP2?Yt4sG|*mJ2QULS?$GuO^z5^E#f-cu)N#i=6@cfyVrPOI&m zWe>bcsx8!$kVx1+#;NORo=h!;^j2so4-eLNmK0mOEeV-DNicKtP{!Vl3=j4om{3TH z7Os$j92!5D5nw#yHpv`%^&-hAxQa<@858 z2E9@Ak9nCE_(<_|%p_Y`+1gPdO0kTE zYna|^NLfxx`%YbT=gpAtkHsrk&sj|7n?BQ9aS?62TO%B_%i>fzyKvH634c@01G=$t zrJla4$L2LPMw|md%>1>d(V<m+r5o=1@W>mXbKRRUDQuDa zy5>aKLOp_hbT)k4d56X2*Bc>7Ut3izcgO_Ym7SUOaAGHf@2xOP4hbh!=QJ~NFOeCM zBX4V)zQfcscg?Cbot}DrR5|n|H4DL)!W~*-fr()FKceJtCj2-EKK%yhh_;$qNyAOt z>fuv8X3_(sXQ3e`le5gXm^eHY2~8#54L^F_PZL%(0Tm%0YJpO7P;S}o9gcO>v`j~? za9Wl@93AJ@jMKX$i#ntREmKj~XbW`Jr|)TQRtcRgjHd&s1nj}EE?AHIFpe(}LS$IT zgwqBjVReO_Z8s(zHk*$=q(sj}J?kENf7RzKrC81_JmosnjZ#PLX&zze_Ql0bpda>( zT&rcaW0O--eZA`SE|FzZvG@2k%J-T%h2w1Zvz5b0`cokxI zX}nF5Z>MyQ5K|G^K35ypi0RQ7t~H*w~YuXZp*=Kf0e?1W*`betBA_DvgWkpQ66bonilxm$271F;m}f?6_loWKl`A!+I(!Tx!#@TYMsi&y_KM%$WZAJC^3~B0y~*;b@(I7o zSqyor=!hf*Q)V0HRzOFzRxPrl5!))pedO?$Ag5KUAd;#uq@=c1(CCvJKQl;Wr7)lX z;nmgxOQ1aDx&Q)H*9x8=U5+f)UYupu)#Vcw5&rR1*hT*G0+Oma0J*3_i}t)NWqA?e z;!L~>qSPD#h>rN(p%yE)Is{J|ogR0}n)k2@O3+v*!f1D+9Q*D7N>vvIo2A^2!TY<(UAT zQ(Lu2_8$v{#T=}{%e)H0%SNC7wzv7ye%ZeX^y<$*&;Iu<6k9$%R9HWOu=PiP)IT!O z`V*Pxzd`o$f8WAY;(J~V^b|Hu)J0^q+?0-owYGt(;u>BVoUJl-USWd+zf)`81MFLnlYyz2`0-~Q zg~@KpYg2B*fD^Ll+?s^;6D@UuXZJ!~Xiez8?g|Ox_!|HA-q4WBe!E?3BJChQz4egq zIuX(Q1KtY9&=LoZLC^HIE-)P(0Ji8NAT}NNs#rJhT9`3jF&QnkNbi2U^o^stEbVMp ztkPPZadXd&)4R5uwpGg3T9R?Y@WzQI%R{ybr!Y5a{F6txeK3iGNFK(JT5Xa@!)WQE zuUeyYo~|H#6^&0t;%EvFZ@W_C2cfr zv_HY9f*nzj-~?_G`u{-ur?Y>TFtmE3*upY;S$A}_m^R$F6I7Ai*JS<0@%S?rYKel| z!LGRp1{3_Vh<_vGzjFTvpKDifFb>}oOT`uI9(=uJ_x{>n+hNHXwb^f;R3lAYOO{rUJ08*C|twJzCHcyalvF^4aL( z%PT8&@fC5!E-6V2G#&D8;$+0wqBJA!O`U3@+F~`$9YXI;nGKC(_&{!&HB2UoJFo^+IMT3R~o0yA1!DCDov^C0F_$)Jw#}~UL2yBF~7p;Wx zMeK5XJ^~ufI|qtEl@X>uG4a_bAD-fuMbMXMbK&GV3JV8~PYw%Dned5K4VMF(>Y>AL z))C|T!G(ZDNKrR81XtWGUJ_~a^+eNXMyV|WQs$;gnbzi878h!N(HOcK?&yC`DF4Ab z$vH|q2{iVjK9i;?YL%^j`3c&bffjZt{sCh#NnExAL)Vk-*o^9esxwW6`I*oQe0N!b zP+6SKg3$>ExAT6%kvE$`_;5QYv~WK1+?Y2Md(ZS&e%nDP=ZcvykDpy(4$_s-(Ol~V z!zevGU8}%v*JcnPLpmAOK=b!fAeA|P&l;fqs;LR|F~0K9t`Jl={c`Iyo!8);=OqE9qS4V4CAZIqwFiBoTdhC#!(Fu>5?n$ER zXAlF5nfz?JP8X>EqbZW%vr0sm1Zei9H2@c2s;xm7wGR{pzuyN*u0%XC6Zmv>mLt;b zD_DzF2cZuc$`W7KCE!s%h3=gTnbw;D^VVBI&_4*;S&H-{(XSB<1Qs-nWE>rT45b~D z18{s%>cV+`;HAY~zxrf;x#^FQ2!DHy(6{oL1<~*YLil2y^53GEg82C`TR@3d6Er-( zU&3(Dx*bXqp&K=*O zL2)L_w~8}K^RH40ob>(xsl3@BiW@6nkN(&ekj_BHrzq6>(Sq3XkdKEFBSuTb!D6)g zNjNCKz*LyR_YV2)GFru}yPSNiCGEJg^k`FE>cwC^@Z5#d@!EXNTMoTSu6^fA-PrSV z+3T)?s&B6WsM~4l5}xAE?7Bz_5N5^%83&MRiz)Fi8GH@EYh%H!i%YA88RXJBVcZ9# z5!J^{fd>B0Lj)EKh3tA3ycg=RmIF&|q0M(!!#l%4D}F040_ROl7Wn^Tf8z)F+xE5( z5d3Rk7;ZKRnvU?oY_!a2156a`N5qEDQlVUm|H)@$e`MyAvdjZ*QL*%6w|sr zm|63>gCYY&^+WON=O>5+MAPhaoZC3mjBE+1S`mp;K}|$kiWcY2R6tW-MwTr}laqq% z4emtKli*1=atZ!3tNx-^r-CUQfmcbV2Q!PM80-{+m zp0C-q{pWMFI=6=!44rLix20>Ar%oN2FIeY}lGLe?CfHU2sZ>Prl;Kc&1|3S#D&Icm zxA|56V9EKqiJUr0PR{0#kl?k91KXO^J*vJ zcw8!h#I4{mC_LmvMFFCb1h4>*prqluHCF%vXJZ!%$O#C*9DTtR_9E`ln&h+k4~kF!0LV}LpN>H(^O1-`%jN^JHwBCyXq`BE41vn z&YLf+&;<8b!e)K=uDoFie$VUGgX08O_uRtDcV9D(Z+Jg>b-R9=?bWpf)WdfkODfV9 z`$(UpL?4KMuZAbeyl78f6)Ar2=*+`xk&Z|5stJ~hNLTNxJB82WHNDP%_^ON;_&zxQ z#jH}V&vOZ>QLP)dx?4ByW=h{ENSHh|*X(&$ze_E75HGV>$wS7?O1m4tpPt=d7ua-c ziNxV_qiENu!<+R?&ul#yj}2k0=!Pz7VAik-D982VA|_7C1~dn+GBo@0YW*Uy!cySW z2H~oGql+DnzL=?s*Bv}a-@nN7hKl{swAF9cCr4z3Sct>D$pml%?PO|jVfX^k4#9~G zby3|4!$A#~@J+j8pPbHA{cM46`=mYIK+BuX5WtI;ugV6!>lby)h! z#*I6by>73^Y__q`_q;DZCKkQ3vqdJodaDD?X8YhJF8S#a$!*`x6Jy)9Wd^+Rdc#W% z#y>L^4{@w5Kl6zJitLy{Oa zHLgW99*HZ}YGo2r9_^lbm0_eb{-%-X8Z5f~7+3!bN$D6-A?`-#1_?(#8) zP=BAv1Do@{)HKW!Mg_K?MwK#qqUJICZR~3vZkwVwJgx`I<3+|q zUiga-&Sf8_Me=AgB=!V8l0(I_>9r&di#NBc5nATjnXF^dHAIPwl$ANUOjn9d5;{9Q zczylM&_E0>D~8;;X6)?9iE|@Q6=G~+%D43T`fgKN7t@jEtF5#;=1ZF2NhOn*t|Pv3 zN}4goQIC?6=&kkR4t_>6B@=%5#3GZ5!<-REpT_z+33`J>YPGR>A5+-loQmRq!j*x#__k@ zk8%jz?<%G$rn+e|>~1D0s<|MUCw4$(7G2rXc!VCmQ4vy=lh+ehxnfO8h{!q_i@TB@ z;=w^noiQZJ)Z0tY8YEKEjm_JQ?+QnsR8&@)5x-nr1#@)46?%gL4^B=+wBW7`24`)& z%aEyG_sqJF@~|8U#qKt~iv&tPsPfnAKf>VoVKB2l!Y~gvu8|^dooqT^I(eRY(t&(7 z_QW|A`Tq_TcO@vgcl(mr3Whl+R)l2Ug_S#!QD*xCLQ4YHiwU)CG&A|8BpGvl%W%AL zzxc&1LpzN>ieKL{0>`4kMFzmu7Db{BW$49V-P0!utVl+m$>jEpXp53|<<5OOYe}rW z-QK2;tfGTk>o$a3_@*k_F8~4BvY5f0z4v1kbe^n3&wU;~eK2w;=XhS(lrXhK#W`W| zey%Gy^lGysc4o{wEZBZ{D#Ll#)eE-8X1PXRPRzs$`A5O}>p)qP0ByfrweGDryVIYq zTG^IAOC?w+w~5!*gce1Np*l!MP7v!ogLH}F38^DH545cd2#9nFd-Y!V<|c9I$Vxdw zagE%OqVYM0AlTD5jRbat(CKjM1z0>M=;Zsht)!Zi!% z@ZfyP8<3PsVoU8R&X&(D-dTHY_6Gf(k?_ceIoz5S- zSm!Zw$5urI@>7;`ijm9;>NsgJ={}CUwbn-Bj-s$#FUHNM=vqy%xIHOrTeDMeo`^*$mo6zkaw{|mjBGQepNNOt4{DkBi03O!|<&4u_i7GSF zDVQqv5qXWjQUwwcckeVp^^$0K55o1-T(^)aT6mRAxBu0(o3L_tWVzs<%Q0!M0A<`$ z3Z~Niyn^swD9MCGjmi-UiV`Hne%Q9DEN0}k?6+9DUb93o{XCs6BPuvGFiD+$ymnLS z=VO}Z6d`?0xpj(mirs!*ERLW{LX}e#H}_(W``~gy6rt=eF7ndoa^GSrx%>2I7{dOG zz}$OfgCFuS8r=CnfPmeoT*NM|MB`pSAYSB^NkFJ@S&h#9^5c=Aqr#<=!c>9%qLjOT z^<4i@N7>5NvPd^Fssci&lH;{Tx8}U zry`DO&3P#haUw49r{W9|#YLeNE!dX?Qa8qoMbu5ND-fFP!i!fF^OPJ)-6Cpe;WY^= z4%)4WhwDWC(7P^Dd?D<{wAm%16M0~J|A_ku+p)L373=a$DM!Q>$4HNCm%D?qH9LAj z;z!Nb#7K(Fyi~36S$K`Qyp0vj&F#V&&)Y$W9UX5O?SEmY*8ZLqx!hFNi&p z66C@{3?{GauN-^Qr{K2kKd)K+uWLqg27mDWSHV0m$3KDT&I9vMYBL)tnnCUlCR)`g zIfh}w=GCdvVA7hnb`$1v@bAuEBd)P%Y%u0tP|e zJUN!koS_3h0XgsiRdxsLsR;ktr1r+m#VcS$4|=$bn|JHPIh>s@hV2`q7=J-a ze}h`{9hB{xC61QfW2amR=As9E=ISYiv<6vE#PjiFQNdTX3;zJ5_Z^UW|0OL|qeuim zS#Mb7Y@0>JoF1o>x#Hh}TK!K@yT{HrTF$EqibqtnR7>>_eS9YtNg(Ay$1H+*@(Mx} z%Ng|_a#;Om4r~N+NDbWG#%-d4-_ovZM>v`o$Q`$Rb3OIK69fAwBNFL4C2pE~ znPEnnyVA+L5)hpu*^C|#)(MoKf)*Z)So2J z&kRhzVShrFDv`U&a^$@ZxEL?dKYW?|E$`@jmmn*>fZpny)|^6h)Gvf<`Bcj~Y;j{& zp^0W-@w8fe&|r;(+}x$a&J6`ML7{WHdf{hVLJ!Kp z5m6#+>P(yL6*-7(FD9nuoYq7XHrx#@W;*%DMzA*RRWc`D@Vc||$&R^muT9b~rOc)m zTa?d!^$rU>Gct9{*?IQDt5+u{zWO{GcU?lln=^JeIVAI`EUMEhE(J?$a18N2bNjH`f8l6tm?Y z?rP7*M`>~2xJ0}e9U{Lk?-Q+jnkw+bfE@I^qbHGoZ-3!zlaEwZEqMnk|=2L;XqI$FW8Bq)FGg&JaiL$Jogl8(e$3v(2Ye!j>a8vlj(^ zPY)f0g4mRfRMKNL+Bm#1pq-nFIWsg)5>zorv12#V6Ra07(l zzrw~%DOV}w&zan*Uc|f+F|V07H~&|2MlWJ+I%Y~ahnNriHm{sFuT%cz9nn7DdpBav zqM1@~lw_kXWZo(K4?-OrNu#&(Scr}pcAp}C9wyc*+I{@ku7#rzP7T*7^_nx1>BVh{ z9u@mT#~V9mz80rX?oKg1q!H{DZD#rV6se62K^syCnUfHPL<%4vwVvut2`fE2vM`Jm zr2TE7&ym_+XGR@y4cb>DdNIYH96?2H4;IHZIFD&<5c z;YFJ6oI!fT22Cp@@@w4ghSWH({q-;3THk#aQv8XWR$TkY`v9f#qV}%bPLr43VM>ck zUwSM2;%60~A4*w!WL@WsGsSBjbkkrM(t1G`kTB9`FLZABW6-0&=Sbuvl#b>Y&7RXe zgk{xQgdd_mJ3A;rFN`7H%IiYlAlKh;%)~VuhKi?W>gYD1b;&a=3gPFKbYF~q?Qagf zy4!-p?!h9@Dkbm-DPzx(v;J_^E(GI9)vn(dk^QSEg_`gMn;nz%O^3I;$D$&_?>6VV zM3kT8BU~ZZ-&`RrF%sf^n{zv~h*VMX(x{n|R~gBg<`28u7le_d68cAmZGI~j*@Pe5 z>@p`-MiytF;m@eAVFMHg!WMLe(uBf^v+eJTGbSff_kYR}(%Afa?DIioQ1^YpuDTpH zqQlNUM|V541OHth`y^!AkS7Yr^^d;ONqaSIMl9mmnRo2Nx-y-&c9lPm|EXN%Cpq(3 zXH+a^11&go418bn;gjC|`1 zVY1onn^W#Sqrbzozuza>`Vj$zT>k*2oAye>oS4bAv+dY{fKvF_Wj+pQ?0cK#bLX>{ zfwS#_JxCbj`a6tXj3hkO@nB9z!~9ut#RTzDxv5VUu}%Hf6hiRX)Tip-vXM-X>+ejE zVQpXH>Bym7Heb)x{v4z|CyLTY#(gHY7Pixnk|a+fCApBi{tun;PdK*l=o54g)DnwhgUj+`kWXN(Kt7Bh-kq;85r+L5>7ucy z(x{{~NPBUmhca@G6Y#v`mue1GDXq{zaua|@&rw(oGkvR_1gPd?L-sYXyXPbKvVw%a!CNy}BaxYZPc2-XMS$r^@?8(<=)lf0?eFk5_n`9^pKT#9QRK0(bAHS?{u6W*YNm zGI8_>VZPaw$)uyi{fj=Bdr#-t$q5%;{zHhWXPAsyMu3?Toq3&Pg(Al?7BIklrx9j6SJg!W8<>XR$m>ngIkStBbWqUhC z(fo2^_P*54JNo;GuY#O&dh;LliL6oDfB8Xo(~a+SI3j((Vs>0jBZ%xBvgkIK8`qH!hq%(kDn4llX2jFhlfys~ZB zE=!a6UAwH`*$QW?@6EXTG5*jln_*kQQ|dc1j8~2fV{ZJ=E!+-XwgFIs4{rb#stO%1cPU(09(sPq*UP1+k2^|FM@ML9 zeZ1-j_ZR`n@TQG`xN6iBxbAGmexE#bbEgimwhYJj*|O{j&gBPgd8-F*g^sU+fenDp zO3C9s#`E3Y~| zz10uxl53%##`3;68h%$FA1%hPbjWMOOVHICd^D4wh~euj&0sD7jl=rZeVn3@;Uu5215S)NNgt> zSgJmc`0>o}lCv+gu3)i6PZ1Uf=-C#AlobK}>65=*`RpzDoDra-Ci*6{*z;E7M2lL4TE*SKE<=ZG;%vAI`C%1@=sD?Xf^yEBB& z77k!&wMV9Ld$t}ykK9%;_do4uw@qZ1-MekZ*($d)uIz#2ci9hbGe4zrJ>#;W&GB+Y zqYxgGivJpRBM0L)(teM`+s$*tM3UHkyXFEFx5L9VK3Q(zqmTGiTdRb;!wnon*Y4=jC|{#c`3Sz)i&*^4mBQa#Y15FK|7ZGG z&CjK#+yz@i;rP3jQCwFIUo0$r(3!fXgPwb=CsKAQxWqT<3nbZe4ys=xw{qEh-ogZQ zTit@S;SmQ{(cD*vZ7;9m%M|vbSsmwm~;D#tC zv{(%i3o)FpGfYOR_J5xjZSYE1aj`FnHBb^vD&6%_%nN5#iaA(BajktBq4ae;@JJD* zX3xuPL0ESiD7RW|K7TPctarka!?SbwyE^LD$&BKRvk>V@Ad_NzM@FdMm%k}qj6c)< z&i`+%$EMo;94zoALqJ#_X=r2b3B4ir&jQxvYork9C%Q$W6o|(sa7_Ou`1h@c-=Pr=xn^ywM%>$YQ-Ndam@eWxH!5%~6JiYkTI4 z6$v$Sm80SGZ~G2!`VXIj9Q5+DLHifMW4sb9~*?;|4X7))Su!R7YC@<-_?8iV0YCWfoY9t+=(d5nFvbu z`zWT^lO*LM0~L;|T!NBFTDRQy2Vdxh({2K1vg*tLv85w}6)Y-_TEU7JM;DlL&cIPT zU@Jf(x_$#Z!R1PBDeU(02pos@35RN?4PZ_=@*p_P4sZk=yI`6numz602v`!k-GI*k z%n$;lrMc%RleN$3r4l4|+?0Ps==$c7M7VJ;pa*Xk2A9E4h`?re#1asNhh_m|seT_0 zkvEOuB5y5{c#|Y9T?)H!yXN)4rsOAc~1$)8=_T2$&kQg~6z@XP{F)4rIQCne3`BkG`f~o8v zqW)>16fEML=n>2wVmnhw0#As<*?mZa@R(@&Qa^<3;{=54rCLUeP7EGLE0q`;Ti4#Kg}4B=5;r{owgbea0H_WK`~3{z2*iEo zZW+F!n9gPX8){YQh}IT5pn3kk7QvZQfF#lT8{h^CPxB=TKVt(sfn1`+98xji_WB40 z88&i8bclEUUPWO>5{s4!Z+ecHRSKbtYw z`d(q3+tD7s#JPGoheAT*+3vbg?Ti8aWVhCW2nxrPIBtg!uhrFoqBl>B>if@3wQog;m}z8!5?kncl=PGThMT zcKnqKx$&%4?-~<9g?^7Qm~4X(k9>THUrQkeG`amdh;^8-wz zxledAByKH_iRLl*-!tU&Q*4M;Z;;l4O`ikuFoppv6gGH26`Bd7zX7^SdCz1jIjo6z z2lUA2*5DP_Vb5>K;Ow473OQwov=cyw&7LDfMEu}|02??1d-(GRnPVv({K2&gG&@LC z-3HhJaO8VH42I2tDz_t|�AxOj-9%?A1cgQ! zd-T6o8s{KX4-VM7Uj)@J7NF>3F-?08U>Ew&+5Y0bthrM zOQgqCT#!+Drx{oVE^n^G{u3$uEd;QH%`YM|&Y>A_p11P(M?UjbzzEU2nZFk|*HQkH zBPR%uB04q$hxk_W<&XS1n<2`HusQg?A7BX-64PD+H-H6m8>G+?HSlITKwh+?pf-@r zdSs7$H^v^Y1g9e4@rOX50>5$keoZN+S7FwVDTPPDYjJ=HGGP&4-k{Wf3I_iT66v2@ z!g=9$!xb%n zxfEl7fp-<_o{0$Ao&`^PQgFN~xbYKki5Twy+yvlm7hvluq%|D*t?>sp9$|BzxDE8r z$G#wy#d6Kd?-1M?rr`1WNMP>|fGCJR1TQWJGC;gMakDSK`+RQ#*R)dqrB-A_d&mH< zMm7zlp_)DkA}{qw=6-{OP(V*=^7n?!Keo7q(moS?Jdv8io7w?kAQB#J0vzE^5+Jtl z7$un=JXGvKoc04;zy*YPu;1^yc>eSLckbu!@AlG-Vlf6(`5V_yFL6U-O&T}hOqf69 z%M3a5CzXTH!~KD)g|sT4=<&88CzY?PirvP74cP9IBQ))Th+Ww^-k(f!*SPdZ-Gev# z04Br;I>I~=nHD%LV>W&t&Alx@!z}%4&izn<`%z#MsD8DKLMJM(hOa3UYG~Wfw1+*p z7I3(ReM&dqbh+s|p)Fw^k(TfcHE<2Qb2)g07_SH(0*t?uhllx5I1|MYlyci7_%p35 zakR^eWlCZZhsx+(;S=(bQe0t}nRN@Wi)DD3eDBHS$tR`F3R z2TVuZBK`h}!jxWWRNH`T#}Yo6;5nYhouTpVjweD*nhO$B$WuoVKFjY8!Q4nRJX=?wt`sM48b``#$!C48Ps(=n=Z1c{0)M{gnTAtxK=& z+Eu1#-*u4=W1k9-u{r+GS@g)KLV80$*;EOOcxrN9k2Z)Ncf;qlATR&j7YCP0Omf!A zqswFAhFFPts?HrnjuCii&>7%ky1Kp;rt8s4u@*$6E~h5fwRf5pbe4KiBN&69LI?^& zJ)qOJkNf@OFYQh95aG6*I*1GRu8-LTUf*oa$N0)8}go$yA;C+D4k`?|G zV`Ahgiw2XCcHi-(?EP6rgD1ngJ%mzje70-EIy9J7HCbWmV4jkd@_s5X)5E|#b*Sd& zNwe$y)Y(j&L3!#>^>&oUoMAInx2En{zw2xybNU}bF2yHpK)bFCITdOsg*LCTmqPC- zGlc)z1}MVaDCCJ$yE3>4I7~#a%Yue=ROESrrkC4^(lC$0FybiELyD&KAT{$cPIm}J zDNC8W`@fR#iQf_~jKHNJG<1@7B+ZXv^SgZAK2rn@_&v2qID9RF2tKzB&?X`kM5wV? z0ZPF4WWf#Vzm&6QdMj810j!EJXal;Sv)rJ)y!~t!&L=jRQF24MD{`!i9=kh^nyqDE zo;>tzR+t-wZ>%2ZA=>UMl7e{!t7!W21zjnF8W_4e>Kx4X*q ze9<~B0}UkV8FeSircB&SO#f|Lts8Ztxpd%2;ocOSPt%}#VJYrJDXtyq>-uMtZbnVj z?Gh=BPQT%Cv+3i=m;B9>hM4f?Xh0F~v&t=2!zY{+&S>t#1LGVXd%Bj#CC6acwIm+1 zmV)o==0&n{FcCa7Zw8OU&yFxKY_)S>Kv(4mPkCb847~--OT@deB|_`f+LQKw zp_tt@N_AEaW~v2varhSz4@GkE0^SlXr|{h;m4m%Qy2r^E z2ve2Q=-AP7I-?Sb^@3zw2{W$S+->C|?v+w=--n!Y=^@?41Ltb77b=a9d*5_OU&vQ69_$MC(}CNQvYJ|LI{2A33~x-VIu6I1>PzL0@dRKEO%TRc zvQ9`!3ke8pT_zwPjeG<;*b^Oa@`%6t@!f&`r@TD^^n$%SuN~M-Xy357LgEvCVYtua zoQZUGJbd(Iqg&Y5hGDIAcvN*4N`f_JRzCS9;O#w;zFm)6JqBN)=R_ATSi?vco>|5Dn&o@_EdC{S zmq6#qsP%^~YR^hr5#{63jy+u%yhUd36a8azduJB_*TV-IS4yi$pGlz*cFTLWii+QN zzmmN5Ix9(UY1Kk~(Z`b!b|-5DWUAkM*wov){DN=rgt||gy7!0WOVTe+e6{wN##eaY z@w|i=S`Q=KvWmMFO&{CsmA3rjY3~uQX1wIS&Gy%?3oQ<8%+xv3d3QmQ=mb4 zu$$X-NEhv}aYN@qXZ35RVkP4)%{&lZK+U8+KEJ$bP3iq>#ja;;^6{@ONDw`ro8XIu zH+JB?4}T7UHtTh+QZ$k$y_4^GjCEG))zR5{?<=X~;9yB=%X8M5!rosk6EowSVH&l9w99U-Qccl)8~(pai^DFZ;5?`O&!QgB!vrqOLo>?wlIxuWa6a zS@YqVe364Pp9eG+BoL#&9P>?`vfbuMy@g#p*S5=weA7ON6Kf_*ti2`czuDMRFez@6 z?O1F}Hh*^cj)ijoCZsc~=HDaF>TvwGTl~&lDW~D-svpkRQz+Z*!e;-AYqzB=zOzu~Oo*O+ z=n)s-bkS%akfeT0$G_h4OA+ORUHO!!;$h5{pvWEDohc&@IcqQXJlVY9*(XAeYl5Jk zwA98E8xE|txu>Ul`9ZGVgO`~_-9B<+o(Den#V^ATRt8S5oxSHbzFFrwbgf$VdipC( zpJBD3js4g9H!jjkUvN0=ni@6+_jyZ4Ub{=?$008{RC;LH3Mc7K@3gQBW=bnFFWtH( zHTuT=WbFOXDCrwIPhQ^F69uRohuLL3p~LZKo*sI6;M!~bj-!*dPm(>3O(gG|c&VOg zObkesuRUXO6)$9_N>+`$zBEAd+Kax*h07bPPoEZ1PCq0$lFRMuRc85~_}cr3irt|& zrgce-XPX;K8%Mtx$It&1k??vJCTNz}yc+Ze@liLeu1IEEe%_>`))#JJ)6sIPuCGE4 zS3NLP-A;DcwcYxgmF?)(xtD7zW(~FiV%h!Tc3dWd#GPr!L}R(*b*YhD7KuZZ!1MZ& zl8i8HE}g_4#c>F!-k1tz4}tfA7aBZSjqmK>OpG`p&nCNx6g+3t=Oxa^hr%ABuE#nO zm~it|W$bF6>l_Ln%9uoP`y4U6o;I|W#pJ7621_kwe{9ESoDkWDJ2p^g6Z1jUA40EC zu?Z5y@6@2~@y6N@RJXGvd?IJ3c^qyWo<$n>;mzTdSt09FlO(8-4!F5K1rDP)RbhlW z^omtm=-pLbP$o{PP&A|PI`2fmx}06neY2bo1906$Le#cCR+lG z;9)&bOQo~J^=VRb@=RyPq4i`fhn9};F=CRB+maMQN{-(`rIpZOIn|X$c{$d2`EkK| zE4*DY?K$DK#mn$W3%|#MobreYW@}A16JCo34W(aOG)E$jpB@D877>%D$HJ#wz9WHC z(4eRk>Ki$WDwTDbHrm_cEPI}Mm3ji+ycX0_zT1dv-5;7|+Ls}xyp0%j42;p6W8rV#^KSlU6$dk_Vsvex-?ZQvIneJSjeYfF7|*`sG)&A1U; zcaS}lDVOrN$v)%51Kr`whmWYNa_Ird9^`wEC}j=j&$2;>#0YynPL36=dbhtSiE0o> zdE7o1Fgw`N9X|a~i(;{5-aTT;2t2fm@UeMPrFduv!camj0*2pOY*q5E@YCOj4bOfn z+qe|2G6dBW)+e{t=(J3k?Wy<}7?<^s?BVuwNI~|z0>vs08e6$igP2bkdT#2LokClu7M(df-GXx!>b?W-n zD&}qUdS)T)w+_@mWj5~$d3ZO)hP-13_tm#S=KJhS6rH_!goqZQ)}gInpLL)vD#fUA zhZEOqbgKGv+?mE~?eJrapKi1jdbh)>7;xp$%;pNdr>rK#)0_Umq4z@~WbQ=pHv);x zL4K@@*+VCB2I71yDB*KMW5K=H>B)E>?i7<#%A01ETSRgOPU-RPqxQ73se`j>q3m~~ zczju!$w$Z3G7ilv8jT&OgE0o+8dH>yyI>in82POnoq?_j7UIh?t@B>p-tv49uI7;4dsatMH@&-aOij-FkdH=X*c@N zOWjnx4p#F;v`580`qN?aPOItRH*hF9RsKl_gHT4dh!8E}G%q2ApPof>n56eVbTO-~ zgzp=HYDUFvVX=O7R5=D$v9B(XaGT{efS~ZFr_{9&ryNspK9o6$2N9#*js*Oo<%C5I zSgJi7yfPYA(FfJc^azZt{Rxg+mD!iyASsMYIJhE<^`?neQo5WdJ{A&XVXSc(KMPJi zq2g*}@Z8kgySnJf{j6PCm;m_322f)Ybvh-lV?wk}kbCw^b0VRIxj87s;ZZ{zIjCzC z^Rz4C3)jZ}9B=du8v8ouMP&Hny{Y?1xTJ1&X-KxLW1z1l>|_k8!jIO2%XRmZXO%wc zVg;0Ax$p6enb&sbwT89B;;2Hv^7fY zy0G-OANrooI9^g_i!xyYA;U>BelHehS)gMu5rmqT+nmvL{vJhi}V)_{igNXaR1oDMze zC2>NW_e~6;eSDWv%8l;teWF$}Vy+fK_E4ALdd}4)W(~WsG(Qka@iW3z%$(qCamP(R zxJ6?-r&FTJ@B^X?xk@>3?ee)=_y5t!xd%0MZE@J_BS}d>K?3+FjaVZj5E7&+At)G% zK+y{56)Chp9i>{(S`n02j9&XdpuyZ^tHDaqilTs2Bo{+?TxbCqc?dROnuB3@4Ac-2 zF#=7}lThj%r!&2C&p&I=nYGq$?X~xuIcLwD4;_c-b{)1a8(x0iJ+h%$Tf5SmXU}}7 zZpF_vhH^uO3UhVkzI{(e8ZZjg^hat@$I)7y;g6ExVQ=8Ra^t1ojAzAi>}kbo9@KDZ z?r01bEXxQZin@1Bk63uy3F}=cR?;61jiq&v(AleoWsnjGXhgSabiENZu5V6&*x57f zjq<@X&TRUVC$N(XXx5Rfs;QJ0&CyB5ZABfXTXf6TDe)vH_L^S+tX~gYiJ4SI^j&R6 zxWNjR`T%$Q`PxZ~mtJr&2OCqT+Z~m!!umO}yZ<}vfHuQRZw_NDlsCR0niBlg$ki%E zl<{x$#muoAsyVM30#>H-S3)9;-w2qv#V~&(U?G%Sa#+E}bd?|Z#igG%Z4VlxYyuWY z{an}wB>MC*Ifk*?CDYiOrw^Nafe^NHfR}ruG*a&No8rC-`=|2xh;dPtnK*crcPQK_ zxM`TuUgGK1T3*H^lTvZpB!l6l%Exc^Gl66*2*NI`HN+c|G2NN^teLJMXu|=nB*%%d ztQq~PTcWIFXl-(+xhu-5!-!Rs4^7RZS=;K$$BX(I0~32?)?{-*q%+5kBKj=bUXIr6 z<7$GY zE!?1Goif9mxP&(R14-)cV2qpON7@_U(6(1OdY_Y1Y{Egki;tQbc6vp5BK2XEr@n!- z$wfjp4C3TRN^JSKIx$Cn{L6J>Sh@uSkgB>d?Xb{!ZE+}eKypNVa`q)$911+0b*5KD zonLMB^yn^1lA**Hz+aJBtk)PJU$P?But{j*OHo0Oz+@ezc_H994weN2bx`Lfv66&Q zm=qevm`*kRe{Dv=@Ldubbue>!HUpV23U-ErP~lSsUc3zqNFRk@3x$7X_z%XSf)2qXL2;of1a#A;jRn2z6en+SU>8kvZkg?Q4{El2c?osIPv8^AemNUuky=jB zx-X5q9ls+6CoR|kT&T2Xb>^VEcLKXS+8WpM&Re68RX#c~1!Bm-_iXpkD|M(%j+64#6%Wfr_DvV7&dxloWboGiOWx)1 z*|%2}yH6Guq-6)iCGK10p$p_qfdhL;&-{_Njrng=D+~oxG%e21m_{&p; zBDaay>cK=}@ZkQ{L{6)x8df?p#wi!@9Dm9sct6^ zKuZ?VJI8N^!z?5Mmjl_Jh(9g{rg diff --git a/simularium_readdy_models/actin/actin_analyzer.py b/simularium_readdy_models/actin/actin_analyzer.py index 4c6a4c6..6d83ff1 100644 --- a/simularium_readdy_models/actin/actin_analyzer.py +++ b/simularium_readdy_models/actin/actin_analyzer.py @@ -1149,7 +1149,12 @@ def analyze_filament_length( return np.array(filament_length) @staticmethod - def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1): + def analyze_bond_stretch( + trajectory, + box_size, + periodic_boundary, + stride=1 + ): """ Get the difference in bond length from ideal for lateral and longitudinal actin bonds. @@ -1163,22 +1168,18 @@ def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1): ActinStructure.mother_positions[2] - ActinStructure.mother_positions[0] ) print("Analyzing bond stretch...") - for time_index in range(0, len(monomer_data), stride): + for time_index in range(0, len(trajectory), stride): stretch_lat.append([]) stretch_long.append([]) new_time_index = math.floor(time_index / stride) - filament = monomer_data[time_index]["topologies"][0]["particle_ids"] + filament = trajectory[time_index].topologies["topologies"][0].particle_ids for index in range(len(filament) - 2): - particle = monomer_data[time_index]["particles"][filament[index]] - particle_lat = monomer_data[time_index]["particles"][ - filament[index + 1] - ] - particle_long = monomer_data[time_index]["particles"][ - filament[index + 2] - ] - type_name = particle["type_name"] - type_name_lat = particle_lat["type_name"] - type_name_long = particle_long["type_name"] + particle = trajectory[time_index].particles[filament[index]] + particle_lat = trajectory[time_index].particles[filament[index + 1]] + particle_long = trajectory[time_index].particles[filament[index + 2]] + type_name = particle.type_name + type_name_lat = particle_lat.type_name + type_name_long = particle_long.type_name if ( "fixed" in type_name or "fixed" in type_name_lat @@ -1187,9 +1188,9 @@ def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1): stretch_lat[new_time_index].append(0.0) stretch_long[new_time_index].append(0.0) continue - pos = particle["position"] - pos_lat = particle_lat["position"] - pos_long = particle_long["position"] + pos = particle.position + pos_lat = particle_lat.position + pos_long = particle_long.position if periodic_boundary: pos_lat = ReaddyUtil.get_non_periodic_boundary_position( pos, pos_lat, box_size @@ -1208,7 +1209,12 @@ def analyze_bond_stretch(monomer_data, box_size, periodic_boundary, stride=1): return np.array(stretch_lat), np.array(stretch_long) @staticmethod - def analyze_angle_stretch(monomer_data, box_size, periodic_boundary, stride=1): + def analyze_angle_stretch( + trajectory, + box_size, + periodic_boundary, + stride=1 + ): """ Get the difference in angles (degrees) from ideal for actin angles between: @@ -1223,24 +1229,24 @@ def analyze_angle_stretch(monomer_data, box_size, periodic_boundary, stride=1): ideal_angle_lat_long = ActinStructure.actin_to_actin_angle(True, False, True) ideal_angle_long_long = ActinStructure.actin_to_actin_angle(False, False, True) print("Analyzing angle stretch...") - for time_index in range(0, len(monomer_data), stride): + for time_index in range(0, len(trajectory), stride): result_lat_lat.append([]) result_lat_long.append([]) result_long_long.append([]) new_time_index = math.floor(time_index / stride) - filament = monomer_data[time_index]["topologies"][0]["particle_ids"] + filament = trajectory[time_index].topologies[0].particle_ids for index in range(len(filament) - 4): particles = [] positions = [] fixed = False for d in range(5): particles.append( - monomer_data[time_index]["particles"][filament[index + d]] + trajectory[time_index].particles[filament[index + d]] ) - if "fixed" in particles[d]["type_name"]: + if "fixed" in particles[d].type_name: fixed = True break - positions.append(particles[d]["position"]) + positions.append(particles[d].position) if d > 0 and periodic_boundary: positions[d] = ReaddyUtil.get_non_periodic_boundary_position( positions[d - 1], positions[d], box_size @@ -1287,7 +1293,12 @@ def analyze_angle_stretch(monomer_data, box_size, periodic_boundary, stride=1): ) @staticmethod - def analyze_dihedral_stretch(monomer_data, box_size, periodic_boundary, stride=1): + def analyze_dihedral_stretch( + trajectory, + box_size, + periodic_boundary, + stride=1 + ): """ Get the difference in dihedral angles (degrees) from ideal for actin angles between: @@ -1305,23 +1316,23 @@ def analyze_dihedral_stretch(monomer_data, box_size, periodic_boundary, stride=1 False, False, False, True ) print("Analyzing dihedral stretch...") - for time_index in range(0, len(monomer_data), stride): + for time_index in range(0, len(trajectory), stride): result_lat_lat_lat.append([]) result_long_long_long.append([]) new_time_index = math.floor(time_index / stride) - filament = monomer_data[time_index]["topologies"][0]["particle_ids"] + filament = trajectory[time_index].topologies[0].particle_ids for index in range(len(filament) - 6): particles = [] positions = [] fixed = False for d in range(7): particles.append( - monomer_data[time_index]["particles"][filament[index + d]] + trajectory[time_index].particles[filament[index + d]] ) - if "fixed" in particles[d]["type_name"]: + if "fixed" in particles[d].type_name: fixed = True break - positions.append(particles[d]["position"]) + positions.append(particles[d].position) if d > 0 and periodic_boundary: positions[d] = ReaddyUtil.get_non_periodic_boundary_position( positions[d - 1], positions[d], box_size diff --git a/simularium_readdy_models/visualization/actin_visualization.py b/simularium_readdy_models/visualization/actin_visualization.py index 4ed19d6..7c4c5eb 100644 --- a/simularium_readdy_models/visualization/actin_visualization.py +++ b/simularium_readdy_models/visualization/actin_visualization.py @@ -9,18 +9,15 @@ DisplayData, MetaData, ScatterPlotData, - TrajectoryConverter, TrajectoryData, UnitData, ) -from simulariumio.filters import AddAgentsFilter, MultiplyTimeFilter -from simulariumio.plot_readers import ScatterPlotReader +from simulariumio.filters import MultiplyTimeFilter from simulariumio.readdy import ReaddyConverter, ReaddyData -from subcell_analysis import SpatialAnnotator -from subcell_analysis.readdy import ReaddyPostProcessor +from subcell_analysis import SpatialAnnotator, CompressionAnalyzer +from subcell_analysis.readdy import ReaddyPostProcessor, FrameData -from ..actin import ACTIN_REACTIONS, ActinAnalyzer -from ..common import ReaddyUtil +from ..actin import ActinAnalyzer class ActinVisualization: @@ -262,392 +259,13 @@ def ACTIN_DISPLAY_DATA(longitudinal_bonds: bool) -> Dict[str, DisplayData]: return display_data @staticmethod - def get_bound_monomers_plot(monomer_data, times): - """ - Add a plot of percent actin in filaments. - """ - return ScatterPlotData( - title="Monomers over time", - xaxis_title="Time (µs)", - yaxis_title="Monomers (%)", - xtrace=times, - ytraces={ - "Actin in filaments": 100.0 - * ActinAnalyzer.analyze_ratio_of_filamentous_to_total_actin( - monomer_data - ), - "Arp2/3 in filaments": 100.0 - * ActinAnalyzer.analyze_ratio_of_bound_to_total_arp23(monomer_data), - "Actin in daughter filaments": 100.0 - * ActinAnalyzer.analyze_ratio_of_daughter_to_total_actin(monomer_data), - }, - render_mode="lines", - ) - - @staticmethod - def get_avg_length_plot(monomer_data, times): - """ - Add a plot of average mother and daughter filament length. - """ - return ScatterPlotData( - title="Average length of filaments", - xaxis_title="Time (µs)", - yaxis_title="Average length (monomers)", - xtrace=times, - ytraces={ - "Mother filaments": ActinAnalyzer.analyze_average_for_series( - ActinAnalyzer.analyze_mother_filament_lengths(monomer_data) - ), - "Daughter filaments": ActinAnalyzer.analyze_average_for_series( - ActinAnalyzer.analyze_daughter_filament_lengths(monomer_data) - ), - }, - render_mode="lines", - ) - - @staticmethod - def get_growth_reactions_plot(reactions, times): - """ - Add a plot of reaction events over time - for each total growth reaction. - """ - ytraces = {} - for total_rxn_name in ACTIN_REACTIONS.GROWTH_RXNS: - rxn_events = ReaddyUtil.analyze_reaction_count_over_time( - reactions, total_rxn_name - ) - if rxn_events is not None: - ytraces[total_rxn_name] = rxn_events - return ScatterPlotData( - title="Growth reactions", - xaxis_title="Time (µs)", - yaxis_title="Reaction events", - xtrace=times, - ytraces=ytraces, - render_mode="lines", - ) - - @staticmethod - def get_structural_reactions_plot(reactions, times): - """ - Add a plot of the number of times a structural reaction - was triggered over time - Note: triggered != completed, the reaction may have failed - to find the required reactants. - """ - ytraces = {} - for total_rxn_name in ACTIN_REACTIONS.STRUCTURAL_RXNS: - rxn_events = ReaddyUtil.analyze_reaction_count_over_time( - reactions, total_rxn_name - ) - if rxn_events is not None: - ytraces[total_rxn_name] = rxn_events - return ScatterPlotData( - title="Structural reaction triggers", - xaxis_title="Time (µs)", - yaxis_title="Reactions triggered", - xtrace=times, - ytraces=ytraces, - render_mode="lines", - ) - - @staticmethod - def get_growth_reactions_vs_actin_plot(reactions, monomer_data, box_size): - """ - Add a plot of average reaction events over time - for each total growth reaction. - """ - ytraces = {} - for rxn_group_name in ACTIN_REACTIONS.GROUPED_GROWTH_RXNS: - group_reaction_events = [] - for total_rxn_name in ACTIN_REACTIONS.GROUPED_GROWTH_RXNS[rxn_group_name]: - group_reaction_events.append( - ReaddyUtil.analyze_reaction_count_over_time( - reactions, total_rxn_name - ) - ) - if len(group_reaction_events) > 0: - ytraces[rxn_group_name] = np.sum( - np.array(group_reaction_events), axis=0 - ) - return ScatterPlotData( - title="Growth vs [actin]", - xaxis_title="[Actin] (µM)", - yaxis_title="Reaction events", - xtrace=ActinAnalyzer.analyze_free_actin_concentration_over_time( - monomer_data, box_size - ), - ytraces=ytraces, - render_mode="lines", - ) - - @staticmethod - def get_capped_ends_plot(monomer_data, times): - """ - Add a plot of percent barbed ends that are capped. - """ - return ScatterPlotData( - title="Capped barbed ends", - xaxis_title="Time (µs)", - yaxis_title="Capped ends (%)", - xtrace=times, - ytraces={ - "Capped ends": 100.0 - * ActinAnalyzer.analyze_ratio_of_capped_ends_to_total_ends( - monomer_data - ), - }, - render_mode="lines", - ) - - @staticmethod - def get_branch_angle_plot(monomer_data, box_size, periodic_boundary, times): - """ - Add a plot of branch angle mean and std dev. - """ - angles = ActinAnalyzer.analyze_branch_angles( - monomer_data, box_size, periodic_boundary - ) - mean = ActinAnalyzer.analyze_average_for_series(angles) - stddev = ActinAnalyzer.analyze_stddev_for_series(angles) - return ScatterPlotData( - title="Average branch angle", - xaxis_title="Time (µs)", - yaxis_title="Branch angle (°)", - xtrace=times, - ytraces={ - "Ideal": np.array(times.shape[0] * [70.9]), - "Mean": mean, - "Mean - std": mean - stddev, - "Mean + std": mean + stddev, - }, - render_mode="lines", - ) - - @staticmethod - def get_helix_pitch_plot(monomer_data, box_size, periodic_boundary, times): - """ - Add a plot of average helix pitch - for both the short and long helices - ideal Ref: http://www.jbc.org/content/266/1/1.full.pdf. - """ - return ScatterPlotData( - title="Average helix pitch", - xaxis_title="Time (µs)", - yaxis_title="Pitch (nm)", - xtrace=times, - ytraces={ - "Ideal short pitch": np.array(times.shape[0] * [5.9]), - "Mean short pitch": ActinAnalyzer.analyze_average_for_series( - ActinAnalyzer.analyze_short_helix_pitches( - monomer_data, box_size, periodic_boundary - ) - ), - "Ideal long pitch": np.array(times.shape[0] * [72]), - "Mean long pitch": ActinAnalyzer.analyze_average_for_series( - ActinAnalyzer.analyze_long_helix_pitches( - monomer_data, box_size, periodic_boundary - ) - ), - }, - render_mode="lines", - ) - - @staticmethod - def get_filament_straightness_plot( - monomer_data, box_size, periodic_boundary, times - ): - """ - Add a plot of how many nm each monomer is away - from ideal position in a straight filament. - """ - return ScatterPlotData( - title="Filament bending", - xaxis_title="Time (µs)", - yaxis_title="Filament bending", - xtrace=times, - ytraces={ - "Filament bending": ( - ActinAnalyzer.analyze_average_for_series( - ActinAnalyzer.analyze_filament_straightness( - monomer_data, - box_size, - periodic_boundary, - ) - ) - ), - }, - render_mode="lines", - ) - - @staticmethod - def generate_polymerization_plots( - monomer_data, - times, - reactions, - box_size, - periodic_boundary=True, - plots=None, - ): - """ - Use an ActinAnalyzer to generate plots of observables - for polymerizing actin. - """ - if plots is None: - plots = { - "scatter": [], - "histogram": [], - } - plots["scatter"] += [ - ActinVisualization.get_bound_monomers_plot(monomer_data, times), - ActinVisualization.get_avg_length_plot(monomer_data, times), - ActinVisualization.get_growth_reactions_plot(reactions, times), - ActinVisualization.get_growth_reactions_vs_actin_plot( - reactions, monomer_data, box_size - ), - # ActinVisualization.get_capped_ends_plot(monomer_data, times), - ActinVisualization.get_branch_angle_plot( - monomer_data, box_size, periodic_boundary, times - ), - ActinVisualization.get_helix_pitch_plot( - monomer_data, box_size, periodic_boundary, times - ), - ActinVisualization.get_filament_straightness_plot( - monomer_data, box_size, periodic_boundary, times - ), - ActinVisualization.get_structural_reactions_plot(reactions, times), - ] - return plots - - @staticmethod - def get_total_axis_twist_plot(axis_twist, times): - """ - Add a plot of total axis twist vs time. - """ - axis_sum = np.sum(axis_twist, axis=1) - return ScatterPlotData( - title="Total filament axis twist", - xaxis_title="T (μs)", - yaxis_title="Twist (rotations)", - xtrace=times, - ytraces={ - "<<<": 9.5 * np.ones(axis_sum.shape), - ">>>": 15.5 * np.ones(axis_sum.shape), - "Ideal": axis_sum[0] * np.ones(axis_sum.shape), - "Axis angle": axis_sum, - }, - render_mode="lines", - ) - - @staticmethod - def get_total_plane_twist_plot(plane_twist, times): - """ - Add a plot of total plane twist vs time. - """ - plane_sum = np.sum(plane_twist, axis=1) - plane_sum /= plane_sum[0] - return ScatterPlotData( - title="Total filament plane twist", - xaxis_title="T (μs)", - yaxis_title="Twist (normalized)", - xtrace=times, - ytraces={ - "Total": plane_sum, - }, - render_mode="lines", - ) - - @staticmethod - def get_twist_per_monomer_plot(twist_angles, filament_positions): - """ - Add a plot of twist vs position of the monomer in filament. - """ - end_time = len(twist_angles) - 1 - return ScatterPlotData( - title="Final twist along filament", - xaxis_title="Filament position (index)", - yaxis_title="Twist (rotations)", - xtrace=filament_positions[end_time], - ytraces={ - "End": twist_angles[end_time], - }, - render_mode="lines", - ) - - @staticmethod - def get_total_bond_energy_plot( - lateral_bond_energies, longitudinal_bond_energies, times - ): - """ - Add a plot of bond energies (lat and long) vs time. - """ - sum_lat = np.sum(lateral_bond_energies, axis=1) - sum_long = np.sum(longitudinal_bond_energies, axis=1) - return ScatterPlotData( - title="Total bond energy", - xaxis_title="T (μs)", - yaxis_title="Strain energy (KT)", - xtrace=times, - ytraces={ - "Lateral": sum_lat, - "Longitudinal": sum_long, - }, - render_mode="lines", - ) - - @staticmethod - def get_bond_energy_per_monomer_plot( - lateral_bond_energies, longitudinal_bond_energies, filament_positions - ): - """ - Add a plot of bond energies (lat and long) vs index of monomer in filament. - """ - end_time = len(lateral_bond_energies) - 1 - return ScatterPlotData( - title="Final bond energy along filament", - xaxis_title="Filament position (index)", - yaxis_title="Strain Energy (KT)", - xtrace=filament_positions[end_time], - ytraces={ - "Lateral": lateral_bond_energies[end_time], - "Longitudinal": longitudinal_bond_energies[end_time], - }, - render_mode="lines", - ) - - @staticmethod - def get_filament_length_plot( - normals, axis_positions, box_size, periodic_boundary, stride, times - ): - """ - Add a plot of distance from first to last particle - in the first filament vs time. - """ - filament_length = ActinAnalyzer.analyze_filament_length( - normals, axis_positions, box_size, periodic_boundary, stride - ) - return ScatterPlotData( - title="Filament length", - xaxis_title="T (μs)", - yaxis_title="Length of filament (nm)", - xtrace=times[::stride], - ytraces={ - "<<<": 450.0 * np.ones(filament_length.shape), - ">>>": 550.0 * np.ones(filament_length.shape), - "Ideal": filament_length[0] * np.ones(filament_length.shape), - "Filament length": filament_length, - }, - render_mode="lines", - ) - - @staticmethod - def get_bond_stretch_plot(monomer_data, box_size, periodic_boundary, stride, times): + def get_bond_stretch_plot(trajectory, box_size, periodic_boundary, stride, times): """ Add a scatter plot of difference in bond length from ideal for lateral and longitudinal actin bonds vs time. """ (stretch_lat, stretch_long,) = ActinAnalyzer.analyze_bond_stretch( - monomer_data, box_size, periodic_boundary, stride + trajectory, box_size, periodic_boundary, stride ) mean_lat = np.mean(stretch_lat, axis=1) # stddev_lat = np.std(stretch_lat, axis=1) @@ -669,7 +287,7 @@ def get_bond_stretch_plot(monomer_data, box_size, periodic_boundary, stride, tim @staticmethod def get_angle_stretch_plot( - monomer_data, box_size, periodic_boundary, stride, times + trajectory, box_size, periodic_boundary, stride, times ): """ Add a scatter plot of difference in angles from ideal @@ -680,7 +298,7 @@ def get_angle_stretch_plot( stretch_lat_long, stretch_long_long, ) = ActinAnalyzer.analyze_angle_stretch( - monomer_data, box_size, periodic_boundary, stride + trajectory, box_size, periodic_boundary, stride ) mean_lat_lat = np.mean(stretch_lat_lat, axis=1) mean_lat_long = np.mean(stretch_lat_long, axis=1) @@ -702,7 +320,7 @@ def get_angle_stretch_plot( @staticmethod def get_dihedral_stretch_plot( - monomer_data, box_size, periodic_boundary, stride, times + trajectory, box_size, periodic_boundary, stride, times ): """ Add a scatter plot of difference in angles from ideal @@ -712,7 +330,7 @@ def get_dihedral_stretch_plot( stretch_lat_lat_lat, stretch_long_long_long, ) = ActinAnalyzer.analyze_dihedral_stretch( - monomer_data, box_size, periodic_boundary, stride + trajectory, box_size, periodic_boundary, stride ) mean_lat_lat_lat = np.mean(stretch_lat_lat_lat, axis=1) mean_long_long_long = np.mean(stretch_long_long_long, axis=1) @@ -732,11 +350,8 @@ def get_dihedral_stretch_plot( @staticmethod def generate_filament_structure_plots( - monomer_data, - times, + trajectory: List[FrameData], box_size, - normals, - axis_positions, periodic_boundary=True, plots=None, ): @@ -750,118 +365,110 @@ def generate_filament_structure_plots( "scatter": [], "histogram": [], } - axis_twist, _, _ = ActinAnalyzer.analyze_twist_axis( - normals, axis_positions, STRIDE - ) + times = [frame.time for frame in trajectory] plots["scatter"] += [ - ActinVisualization.get_total_axis_twist_plot(axis_twist, times[::STRIDE]), - ActinVisualization.get_filament_length_plot( - normals, axis_positions, box_size, periodic_boundary, STRIDE, times - ), ActinVisualization.get_bond_stretch_plot( - monomer_data, box_size, periodic_boundary, STRIDE, times + trajectory, box_size, periodic_boundary, STRIDE, times ), ActinVisualization.get_angle_stretch_plot( - monomer_data, box_size, periodic_boundary, STRIDE, times + trajectory, box_size, periodic_boundary, STRIDE, times ), ActinVisualization.get_dihedral_stretch_plot( - monomer_data, box_size, periodic_boundary, STRIDE, times + trajectory, box_size, periodic_boundary, STRIDE, times ), ] return plots @staticmethod - def generate_bend_twist_plots( - monomer_data, - times, - box_size, - normals, - axis_positions, - periodic_boundary=True, - plots=None, - ): + def generate_actin_compression_plots( + axis_positions: List[List[np.ndarray]], + plots: List[ScatterPlotData] = None, + ) -> List[Dict[str, Any]]: """ Use an ActinAnalyzer to generate plots of observables - for actin being bent or twisted. + for actin being compressed. """ if plots is None: plots = { "scatter": [], "histogram": [], } - STRIDE = 10 - (axis_twist, _, _) = ActinAnalyzer.analyze_twist_axis( - normals, axis_positions, STRIDE - ) - (twist_angles, filament_positions1) = ActinAnalyzer.analyze_twist_planes( - monomer_data, box_size, periodic_boundary, STRIDE - ) - ( - lateral_bond_energies, - longitudinal_bond_energies, - filament_positions2, - ) = ActinAnalyzer.analyze_bond_energies( - monomer_data, box_size, periodic_boundary, STRIDE - ) + perp_dist = [] + bending_energy = [] + non_coplanarity = [] + peak_asym = [] + contour_length = [] + total_steps = len(axis_positions) + for time_ix in range(total_steps): + first_polymer_trace = axis_positions[time_ix][0] + perp_dist.append(CompressionAnalyzer.get_average_distance_from_end_to_end_axis( + polymer_trace=first_polymer_trace, + )) + bending_energy.append(1000. * CompressionAnalyzer.get_bending_energy_from_trace( + polymer_trace=first_polymer_trace, + )) + non_coplanarity.append(CompressionAnalyzer.get_third_component_variance( + polymer_trace=first_polymer_trace, + )) + peak_asym.append(CompressionAnalyzer.get_asymmetry_of_peak( + polymer_trace=first_polymer_trace, + )) + contour_length.append(CompressionAnalyzer.get_contour_length_from_trace( + polymer_trace=first_polymer_trace, + )) plots["scatter"] += [ - ActinVisualization.get_total_axis_twist_plot(axis_twist, times[::STRIDE]), - ActinVisualization.get_total_plane_twist_plot( - twist_angles, times[::STRIDE] + ScatterPlotData( + title="Average Perpendicular Distance", + xaxis_title="normalized time", + yaxis_title="distance (nm)", + xtrace=np.arange(total_steps), + ytraces={ + "filament" : perp_dist, + }, + render_mode="lines" ), - ActinVisualization.get_twist_per_monomer_plot( - twist_angles, filament_positions1 + ScatterPlotData( + title="Bending Energy", + xaxis_title="normalized time", + yaxis_title="energy", + xtrace=np.arange(total_steps), + ytraces={ + "filament" : bending_energy, + }, + render_mode="lines" ), - ActinVisualization.get_total_bond_energy_plot( - lateral_bond_energies, longitudinal_bond_energies, times[::STRIDE] + ScatterPlotData( + title="Non-coplanarity", + xaxis_title="normalized time", + yaxis_title="3rd component variance from PCA", + xtrace=np.arange(total_steps), + ytraces={ + "filament" : non_coplanarity, + }, + render_mode="lines" ), - ActinVisualization.get_bond_energy_per_monomer_plot( - lateral_bond_energies, longitudinal_bond_energies, filament_positions2 + ScatterPlotData( + title="Peak Asymmetry", + xaxis_title="normalized time", + yaxis_title="normalized peak distance", + xtrace=np.arange(total_steps), + ytraces={ + "filament" : peak_asym, + }, + render_mode="lines" ), - ] - return plots - - @staticmethod - def generate_actin_compression_plots( - post_processor: ReaddyPostProcessor, - fiber_chain_ids: List[List[List[int]]], - temperature_c: float, - plots: List[ScatterPlotData] = None, - ) -> List[Dict[str, Any]]: - """ - Use an ActinAnalyzer to generate plots of observables - for actin being compressed. - """ - if plots is None: - plots = [] - times = [] - for frame in post_processor.trajectory: - times.append(frame.time) - bond_energies, _ = post_processor.fiber_bond_energies( - fiber_chain_ids=fiber_chain_ids, - ideal_lengths=ActinAnalyzer.ideal_bond_lengths(), - ks=ActinAnalyzer.bond_energy_constants(temp_c=temperature_c), - stride=10, - ) - sum_lat = np.sum(bond_energies[1], axis=1) - sum_long = np.sum(bond_energies[2], axis=1) - plots.append( ScatterPlotData( - title="Total bond energy", - xaxis_title="T (μs)", - yaxis_title="Strain energy (KT)", - xtrace=np.array(times[::10]), + title="Contour Length", + xaxis_title="normalized time", + yaxis_title="filament contour length (nm)", + xtrace=np.arange(total_steps), ytraces={ - "Lateral": sum_lat, - "Longitudinal": sum_long, + "filament" : bending_energy, }, - render_mode="lines", - ) - ) - formatted_plots = [] - reader = ScatterPlotReader() - for plot in plots: - formatted_plots.append(reader.read(plot)) - return formatted_plots + render_mode="lines" + ), + ] + return plots @staticmethod def add_spatial_annotations( @@ -968,28 +575,3 @@ def simularium_trajectory( ] ) return traj_data - - @staticmethod - def save_actin( - trajectory_datas: List[TrajectoryData], - output_path: str, - plots: List[Dict[str, Any]] = None, - ): - """ - save a simularium file with actin trajector(ies). - """ - traj_data = trajectory_datas[0] - for index in range(1, len(trajectory_datas)): - traj_data = TrajectoryConverter(traj_data).filter_data( - [ - AddAgentsFilter( - new_agent_data=trajectory_datas[index].agent_data, - ) - ] - ) - converter = TrajectoryConverter(traj_data) - if plots is not None: - for plot_type in plots: - for plot in plots[plot_type]: - converter.add_plot(plot, plot_type) - converter.save(output_path) From 283785d1bab99a8c34a84e0748e4bcc398c89e2c Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Wed, 20 Dec 2023 16:53:20 -0800 Subject: [PATCH 2/6] move visualization out of main actin script for reusability, debugged actin compression visualization --- examples/actin/docker/src/actin.py | 111 +-------- examples/actin/template.xlsx | Bin 209775 -> 209772 bytes .../actin/actin_analyzer.py | 2 +- .../visualization/actin_visualization.py | 215 +++++++++++++++--- 4 files changed, 190 insertions(+), 138 deletions(-) diff --git a/examples/actin/docker/src/actin.py b/examples/actin/docker/src/actin.py index 0571677..8e9b300 100644 --- a/examples/actin/docker/src/actin.py +++ b/examples/actin/docker/src/actin.py @@ -9,18 +9,11 @@ import pandas import psutil -from subcell_analysis.readdy import ( - ReaddyLoader, - ReaddyPostProcessor, -) -from simulariumio import BinaryWriter - from simularium_readdy_models.actin import ( FiberData, ActinSimulation, ActinGenerator, ActinTestData, - ActinStructure, ) from simularium_readdy_models.visualization import ActinVisualization from simularium_readdy_models import ReaddyUtil @@ -115,103 +108,6 @@ def report_hardware_usage(): ) -def analyze_results(parameters, save_pickle=False): - # get analysis parameters - plot_actin_structure = parameters.get("plot_actin_structure", False) - plot_actin_compression = parameters.get("plot_actin_compression", False) - visualize_edges = parameters.get("visualize_edges", False) - visualize_normals = parameters.get("visualize_normals", False) - visualize_control_pts = parameters.get("visualize_control_pts", False) - - # convert to simularium - traj_data = ActinVisualization.simularium_trajectory( - path_to_readdy_h5=parameters["name"] + ".h5", - box_size=parameters["box_size"], - total_steps=parameters["total_steps"], - time_multiplier=1e-3, # assume 1e3 recorded steps - longitudinal_bonds=bool(parameters.get("longitudinal_bonds", True)), - ) - - # load different views of ReaDDy data - post_processor = None - fiber_chain_ids = None - axis_positions = None - new_chain_ids = None - if visualize_normals or visualize_control_pts or visualize_edges or plot_actin_structure or plot_actin_compression: - periodic_boundary = parameters.get("periodic_boundary", False) - post_processor = ReaddyPostProcessor( - trajectory=ReaddyLoader( - h5_file_path=parameters["name"] + ".h5", - min_time_ix=0, - max_time_ix=-1, - time_inc=1, - timestep=100.0, - save_pickle_file=save_pickle, - ).trajectory(), - box_size=parameters["box_size"], - periodic_boundary=periodic_boundary, - ) - if visualize_normals or visualize_control_pts or plot_actin_compression: - fiber_chain_ids = post_processor.linear_fiber_chain_ids( - start_particle_phrases=["pointed"], - other_particle_types=[ - "actin#", - "actin#ATP_", - "actin#mid_", - "actin#mid_ATP_", - "actin#fixed_", - "actin#fixed_ATP_", - "actin#mid_fixed_", - "actin#mid_fixed_ATP_", - "actin#barbed_", - "actin#barbed_ATP_", - "actin#fixed_barbed_", - "actin#fixed_barbed_ATP_", - ], - polymer_number_range=5, - ) - axis_positions, new_chain_ids = post_processor.linear_fiber_axis_positions( - fiber_chain_ids=fiber_chain_ids, - ideal_positions=ActinStructure.mother_positions[2:5], - ideal_vector_to_axis=ActinStructure.vector_to_axis(), - ) - - # create plots - if plot_actin_structure: - print("plot actin structure metrics") - traj_data.plots = ActinVisualization.generate_filament_structure_plots( - post_processor.trajectory, - parameters["box_size"], - periodic_boundary=True, - plots=traj_data.plots, - ) - - if plot_actin_compression: - print("plot actin compression metrics") - traj_data.plots = ActinVisualization.generate_actin_compression_plots( - axis_positions, - plots=traj_data.plots, - ) - - # add annotation objects to the spatial data - traj_data = ActinVisualization.add_spatial_annotations( - traj_data, - post_processor, - visualize_edges, - visualize_normals, - visualize_control_pts, - new_chain_ids, - axis_positions, - ) - - # save simularium file - BinaryWriter.save( - trajectory_data=traj_data, - output_path=parameters["name"] + ".h5", - validate_ids=False, - ) - - def main(): args = parse_args() parameters = setup_parameters(args) @@ -227,7 +123,12 @@ def main(): show_summary=False, ) report_hardware_usage() - analyze_results(actin_simulation.parameters, args.save_pickle) + ActinVisualization.analyze_and_visualize_trajectory( + actin_simulation.parameters["name"], + actin_simulation.parameters["total_steps"], + actin_simulation.parameters, + args.save_pickle, + ) if __name__ == "__main__": diff --git a/examples/actin/template.xlsx b/examples/actin/template.xlsx index c215d429d324253d930a5f5b5d2fa85c62b3640f..94c77a934ee8360d69145de9a547887feddf7283 100644 GIT binary patch delta 1917 zcmV-@2ZH$T<_zrS46rl_1>}8?;Wd*v2_1h=+b|G*pS1rFEyy1QP-V zkH}42V&T}C?a&~!|Gu-^l$K6Y$D-Db?epjF?mk~mzH~MBo~hDIipQwl2P5CphKQ1IDQM=uQ0;#?{ru@?l$|ny1gwFiRm?lL?7Zt6EK9UU;u(E%U zswM`8>JjNiQBvxPQd4sjg<*nf!i2vjog^j5G7%Kz*++((^|l6`6BzyTe1TQ}9dRk3(pp{v=r1JE~@dl3Fgz z5kR9A={^5(h%nv>LAM`+!;l6@(Pw{VbGTupk`J;}aGDFMEKZc=P(wzy_Hj|cdSj9L3Ug&fq=PLyR z$o*Gts1n4`sT2mJwx<5SvU1qdN`i>JYx=Wgio)Fi34uG%k|HH@eML;=wTgddKKh}- zKQ!lrDUN#CYc!+J5D8M2FF`F!$~B@LsT)o;B86cBZHbfY--LgVCn;=n5JY6fdiDN3 zIzV2fy#%j}^5Fe^b^&R)A{97bbnK|6Dyjq5?)3<5A2 zjwWF|j)Fms^I(9-VKA8_gJ3wzP(29qEO6oZ!xx0eh7rUMVng+&cV76AahbqD|$005W$F#{QYoK(qf z+b|Hl7wA6-yjzR6#)jn}3xWpdAxOFaz0uUjCd@@8DK8EB?;TRI?AAFLF>Hz+&b)c^ zW~kZqPKsy)#&WH4lEo>BKrvl#wam%4uMfQ|5?M!8L4{V3lRa2+J^%P&w$)~BR{$`v6KJ$)=A>r43yOtb3OYDT{Un zzkU1-4rzk@ER*ESXsb&XBXZJkJ}-ZgOcKgEr1SbOg2{%7ferT!JVNz<-$yAMby)o) z*x;RD;||QrWu`SRa`KwCpB_Flzof_IaQ27HXJJmB%zUP6C%A$q6WLl~!}hlz^fo6r zs;%LfFIUbRlKCv@;tGx<;$?{plsUP{ravc>vG)ls-#KjU$v$#)@d5&v1EiUe$WQn} z>$P`$L{o~SDyksby;L|lC&6cb+&W!-fzsW95IMO;`6zRI15X%{lZAFp%fK_-Cx^Du z=-;3MpODrw_4qkKsMfW?XfctK{S}uT5&;)~KkQCZ-c7nl+;-}EMGi`z%J>j2J2q-T zE!=7GgBNavI$@^@C}}O+iE$W@I`gyMeuk!r*(m;m2duUzSZBof@fUDf?C-%cy!qIU zI2(CNObcuyu94*a6^1vL?wSTR=?9beYS2#yS--*IL(vci8vf7np=Mq9vnn%-auo^t9b9)Yo0h-gM~wjj;B?Ir+c7&?%k zyu(?N;cu@8E8V{XlXBsd!Vf0gG{?PUb9NRN_sNOH*!O+X7gLOoKf^Ek-Q#p>EY#Q) z15}UwbkzBMsPZ6yZ-#o)QDLLqpp`4}U%BLL3P|KU5`~@QobXF}AujIjn=PvIG%w;T zkMp$23zFTEV)a_$-++2r^YfSLPXLz*Fay{HIuZ4VC3i0FsK&=wD)VhfgQh7-$qRvD^Ygm6- z30V;hL-~kQBPl6VNvW2;*MY_*7RUPw#}?%99A(7fET*h$@xaY z0CGR&n#!7JIv2GDscotMudFQgyc8f}=PUiuGD+d?fP}yuXh|NErMe-ybXtGO6Bm6~ z;2&CY!X!t%>=jDsGem+E#fw*of^vmuN2-QXg-EVhjrPPzj&H(0$dlwIItU`NVZC{O z8yz4o<6eR{S~~E4nO;E}zL9503Z%Z+-q}kN^}$1(%NRfWoX$qkaFR~EQ8b!*cp9eO zbPXz7}C{R;NumKCowe)`+TGyR*U zU;AYoHkIv;(Y1P2h8a5DJ*x|`wjRCo?>HEPgKN8T&k9{8?qE0$cXOX9dn|!2lsq#D z@;9^+_kd2ip_9SmndX2F+T&vjDDbkqQ(Mv7q(h-LwbcgLb!2Rgke7c2Hva4z+wALh z*@9g5_CP_~?FruQd15tHrdIEc80?u3tVVibIpO(@bWCTtOt2qJj0NrVmD08cc{2#W zY&f3zgGu0p8P2>APke7Si$ZUhW`p4fr}M=yGwC;M-;+@g6_*<^0}z9z47a8X0hR*^CWg#7NfrSBz?Vf30W5zG(ZnV(go;^rbadqE-FcRY0!V~kdkGm&LN0lQ}l4=&6_tvO|N%SL>n-cYn78M zPDuoc>4K|OPQHJ;zqlfibyO8pXazahgC*Cq&!47SZJzBKz(t6lY);nBRg)yKYz>my zSXZF%q|`=Ihkwl~v6TT@gdl$@l3tpQ5=prt4Z_5{2Sb-7XE4{SmY`e%Vn9%b{Oy`o z_JAhYd*CEB&vkXdw5$+f$pv@&07@dsCJ!s64P6SXdzTF;i*^RTJ^b|#X@dPMljO{3 zt4kLna?)@Klk{*}C*KZc*3KJ}*)-|m z3XUV;Es`D|V+T=a*jOCa2$b`g%rNqfp| zjrJi}im~T=P|Ro1p`*AnPO%W5u7Aar@i*22F7o^aD+o-m!0-{6+KGtHNWmwheRtgm z5FA4X3Uu(8WCVYE132aX9T=4hr;=f?@uqX!OSa%AadIELoQ!?nXMH)wi1`WLu2)av zsj1Y*t~j9P*iXlNSjQ@N62xL?79Eo=Sq(<}3jdWW!N-79okx; zmv=D(IsrbHvM~cI9VUj%I7t=(0Kjz!02lxO00000000000000l5CNA<5CIvNRWbt* zm$@ba5DNeR0000000000u9rSC13>}6mw++@MFHcN#xesv0k4+?GXo Dict[str, DisplayData]: return display_data @staticmethod - def get_bond_stretch_plot(trajectory, box_size, periodic_boundary, stride, times): + def get_bond_stretch_plot( + trajectory: List[FrameData], + box_size: np.ndarray, + periodic_boundary: bool, + stride: int, + times: np.ndarray, + ) -> ScatterPlotData: """ Add a scatter plot of difference in bond length from ideal for lateral and longitudinal actin bonds vs time. @@ -273,7 +282,7 @@ def get_bond_stretch_plot(trajectory, box_size, periodic_boundary, stride, times # stddev_long = np.std(stretch_long, axis=1) return ScatterPlotData( title="Bond stretch", - xaxis_title="T (μs)", + xaxis_title="Time (μs)", yaxis_title="Bond stretch (nm)", xtrace=times[::stride], ytraces={ @@ -287,8 +296,12 @@ def get_bond_stretch_plot(trajectory, box_size, periodic_boundary, stride, times @staticmethod def get_angle_stretch_plot( - trajectory, box_size, periodic_boundary, stride, times - ): + trajectory: List[FrameData], + box_size: np.ndarray, + periodic_boundary: bool, + stride: int, + times: np.ndarray, + ) -> ScatterPlotData: """ Add a scatter plot of difference in angles from ideal for lateral and longitudinal actin bonds vs time. @@ -305,7 +318,7 @@ def get_angle_stretch_plot( mean_long_long = np.mean(stretch_long_long, axis=1) return ScatterPlotData( title="Angle stretch", - xaxis_title="T (μs)", + xaxis_title="Time (μs)", yaxis_title="Angle stretch (degrees)", xtrace=times[::stride], ytraces={ @@ -320,8 +333,12 @@ def get_angle_stretch_plot( @staticmethod def get_dihedral_stretch_plot( - trajectory, box_size, periodic_boundary, stride, times - ): + trajectory: List[FrameData], + box_size: np.ndarray, + periodic_boundary: bool, + stride: int, + times: np.ndarray, + ) -> ScatterPlotData: """ Add a scatter plot of difference in angles from ideal for lateral and longitudinal actin bonds vs time. @@ -336,7 +353,7 @@ def get_dihedral_stretch_plot( mean_long_long_long = np.mean(stretch_long_long_long, axis=1) return ScatterPlotData( title="Dihedral stretch", - xaxis_title="T (μs)", + xaxis_title="Time (μs)", yaxis_title="Angle stretch (degrees)", xtrace=times[::stride], ytraces={ @@ -351,10 +368,10 @@ def get_dihedral_stretch_plot( @staticmethod def generate_filament_structure_plots( trajectory: List[FrameData], - box_size, - periodic_boundary=True, - plots=None, - ): + box_size: np.ndarray, + periodic_boundary: bool = True, + plots: Dict[str, List[ScatterPlotData or HistogramPlotData]] = None, + ) -> Dict[str, List[ScatterPlotData or HistogramPlotData]]: """ Use an ActinAnalyzer to generate plots of observables for a diffusing actin filament. @@ -365,7 +382,7 @@ def generate_filament_structure_plots( "scatter": [], "histogram": [], } - times = [frame.time for frame in trajectory] + times = np.array([frame.time for frame in trajectory]) plots["scatter"] += [ ActinVisualization.get_bond_stretch_plot( trajectory, box_size, periodic_boundary, STRIDE, times @@ -382,8 +399,9 @@ def generate_filament_structure_plots( @staticmethod def generate_actin_compression_plots( axis_positions: List[List[np.ndarray]], - plots: List[ScatterPlotData] = None, - ) -> List[Dict[str, Any]]: + timestep: float, + plots: Dict[str, List[ScatterPlotData or HistogramPlotData]] = None, + ) -> Dict[str, List[ScatterPlotData or HistogramPlotData]]: """ Use an ActinAnalyzer to generate plots of observables for actin being compressed. @@ -416,59 +434,83 @@ def generate_actin_compression_plots( contour_length.append(CompressionAnalyzer.get_contour_length_from_trace( polymer_trace=first_polymer_trace, )) + times = timestep * np.arange(total_steps) plots["scatter"] += [ ScatterPlotData( title="Average Perpendicular Distance", - xaxis_title="normalized time", + xaxis_title="Time (μs)", yaxis_title="distance (nm)", - xtrace=np.arange(total_steps), + xtrace=times, ytraces={ - "filament" : perp_dist, + "<<<": np.zeros(times.shape), + ">>>": 85.0 * np.ones(times.shape), + "filament" : np.array(perp_dist), }, render_mode="lines" ), ScatterPlotData( title="Bending Energy", - xaxis_title="normalized time", + xaxis_title="Time (μs)", yaxis_title="energy", - xtrace=np.arange(total_steps), + xtrace=times, ytraces={ - "filament" : bending_energy, + "<<<": np.zeros(times.shape), + ">>>": 10.0 * np.ones(times.shape), + "filament" : np.array(bending_energy), }, render_mode="lines" ), ScatterPlotData( title="Non-coplanarity", - xaxis_title="normalized time", + xaxis_title="Time (μs)", yaxis_title="3rd component variance from PCA", - xtrace=np.arange(total_steps), + xtrace=times, ytraces={ - "filament" : non_coplanarity, + "<<<": np.zeros(times.shape), + ">>>": 0.03 * np.ones(times.shape), + "filament" : np.array(non_coplanarity), }, render_mode="lines" ), ScatterPlotData( title="Peak Asymmetry", - xaxis_title="normalized time", + xaxis_title="Time (μs)", yaxis_title="normalized peak distance", - xtrace=np.arange(total_steps), + xtrace=times, ytraces={ - "filament" : peak_asym, + "<<<": np.zeros(times.shape), + ">>>": 0.5 * np.ones(times.shape), + "filament" : np.array(peak_asym), }, render_mode="lines" ), ScatterPlotData( title="Contour Length", - xaxis_title="normalized time", + xaxis_title="Time (μs)", yaxis_title="filament contour length (nm)", - xtrace=np.arange(total_steps), + xtrace=times, ytraces={ - "filament" : bending_energy, + "<<<": 480 * np.ones(times.shape), + ">>>": 505 * np.ones(times.shape), + "filament" : np.array(contour_length), }, render_mode="lines" ), ] return plots + + @staticmethod + def add_plots_to_trajectory( + trajectory: TrajectoryData, + plots: Dict[str,List[ScatterPlotData or HistogramPlotData]] + ) -> TrajectoryData: + reader = ScatterPlotReader() + for plot in plots["scatter"]: + trajectory.plots.append(reader.read(plot)) + reader = HistogramPlotReader() + for plot in plots["histogram"]: + trajectory.plots.append(reader.read(plot)) + return trajectory @staticmethod def add_spatial_annotations( @@ -560,6 +602,7 @@ def simularium_trajectory( up_vector=np.array([0.0, 1.0, 0.0]), fov_degrees=120.0, ), + scale_factor=1.0, ), display_data=ActinVisualization.ACTIN_DISPLAY_DATA(longitudinal_bonds), time_units=UnitData("µs"), @@ -575,3 +618,111 @@ def simularium_trajectory( ] ) return traj_data + + @staticmethod + def analyze_and_visualize_trajectory( + output_name: str, + total_steps: float, + parameters: Dict[str, Any], + save_pickle: bool = False + ) -> TrajectoryData: + """ + Do all visualization, annotation, and metric calculation on an actin trajectory + and save it as a simularium file. + """ + # get analysis parameters + plot_actin_structure = parameters.get("plot_actin_structure", False) + plot_actin_compression = parameters.get("plot_actin_compression", False) + visualize_edges = parameters.get("visualize_edges", False) + visualize_normals = parameters.get("visualize_normals", False) + visualize_control_pts = parameters.get("visualize_control_pts", False) + + # convert to simularium + traj_data = ActinVisualization.simularium_trajectory( + path_to_readdy_h5=output_name + ".h5", + box_size=parameters["box_size"], + total_steps=total_steps, + time_multiplier=1e-3, # ns to us + longitudinal_bonds=bool(parameters.get("longitudinal_bonds", True)), + ) + + # load different views of ReaDDy data + post_processor = None + fiber_chain_ids = None + axis_positions = None + new_chain_ids = None + if visualize_normals or visualize_control_pts or visualize_edges or plot_actin_structure or plot_actin_compression: + periodic_boundary = parameters.get("periodic_boundary", False) + post_processor = ReaddyPostProcessor( + trajectory=ReaddyLoader( + h5_file_path=output_name + ".h5", + min_time_ix=0, + max_time_ix=-1, + time_inc=1, + timestep=parameters.get("internal_timestep", 0.1) * 1E-3, # us + save_pickle_file=save_pickle, + ).trajectory(), + box_size=parameters["box_size"], + periodic_boundary=periodic_boundary, + ) + if visualize_normals or visualize_control_pts or plot_actin_compression: + fiber_chain_ids = post_processor.linear_fiber_chain_ids( + start_particle_phrases=["pointed"], + other_particle_types=[ + "actin#", + "actin#ATP_", + "actin#mid_", + "actin#mid_ATP_", + "actin#fixed_", + "actin#fixed_ATP_", + "actin#mid_fixed_", + "actin#mid_fixed_ATP_", + "actin#barbed_", + "actin#barbed_ATP_", + "actin#fixed_barbed_", + "actin#fixed_barbed_ATP_", + ], + polymer_number_range=5, + ) + axis_positions, new_chain_ids = post_processor.linear_fiber_axis_positions( + fiber_chain_ids=fiber_chain_ids, + ideal_positions=ActinStructure.mother_positions[2:5], + ideal_vector_to_axis=ActinStructure.vector_to_axis(), + ) + + # create plots + plots = None + if plot_actin_structure: + print("plot actin structure metrics") + plots = ActinVisualization.generate_filament_structure_plots( + post_processor.trajectory, + parameters["box_size"], + periodic_boundary=True, + ) + if plot_actin_compression: + print("plot actin compression metrics") + plots = ActinVisualization.generate_actin_compression_plots( + axis_positions, + parameters.get("internal_timestep", 0.1) * 1E-3, # us + plots=plots, + ) + traj_data = ActinVisualization.add_plots_to_trajectory(traj_data, plots) + + # add annotation objects to the spatial data + traj_data = ActinVisualization.add_spatial_annotations( + traj_data, + post_processor, + visualize_edges, + visualize_normals, + visualize_control_pts, + new_chain_ids, + axis_positions, + ) + + # save simularium file + BinaryWriter.save( + trajectory_data=traj_data, + output_path=output_name + ".h5", + validate_ids=False, + ) + return traj_data From fabad3f1be96df6072064d2b8ac34fca61d28f8a Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Wed, 20 Dec 2023 17:01:10 -0800 Subject: [PATCH 3/6] log clock time --- examples/actin/docker/src/actin.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/actin/docker/src/actin.py b/examples/actin/docker/src/actin.py index 8e9b300..1cfcad3 100644 --- a/examples/actin/docker/src/actin.py +++ b/examples/actin/docker/src/actin.py @@ -111,6 +111,7 @@ def report_hardware_usage(): def main(): args = parse_args() parameters = setup_parameters(args) + start_time = time.time() actin_simulation = ActinSimulation( parameters=parameters, record=True, @@ -122,6 +123,7 @@ def main(): timestep=actin_simulation.parameters.get("internal_timestep", 0.1), show_summary=False, ) + print("Run time: %s seconds " % (time.time() - start_time)) report_hardware_usage() ActinVisualization.analyze_and_visualize_trajectory( actin_simulation.parameters["name"], From bdc795712124c65c675d6ed7cc0164efebe34e49 Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 21 Dec 2023 12:12:12 -0800 Subject: [PATCH 4/6] a few fixes and tweaks --- .../visualization/actin_visualization.py | 67 +++++++++---------- 1 file changed, 30 insertions(+), 37 deletions(-) diff --git a/simularium_readdy_models/visualization/actin_visualization.py b/simularium_readdy_models/visualization/actin_visualization.py index 3a87189..f797dcd 100644 --- a/simularium_readdy_models/visualization/actin_visualization.py +++ b/simularium_readdy_models/visualization/actin_visualization.py @@ -282,7 +282,7 @@ def get_bond_stretch_plot( # stddev_long = np.std(stretch_long, axis=1) return ScatterPlotData( title="Bond stretch", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="Bond stretch (nm)", xtrace=times[::stride], ytraces={ @@ -318,7 +318,7 @@ def get_angle_stretch_plot( mean_long_long = np.mean(stretch_long_long, axis=1) return ScatterPlotData( title="Angle stretch", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="Angle stretch (degrees)", xtrace=times[::stride], ytraces={ @@ -353,7 +353,7 @@ def get_dihedral_stretch_plot( mean_long_long_long = np.mean(stretch_long_long_long, axis=1) return ScatterPlotData( title="Dihedral stretch", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="Angle stretch (degrees)", xtrace=times[::stride], ytraces={ @@ -366,7 +366,7 @@ def get_dihedral_stretch_plot( ) @staticmethod - def generate_filament_structure_plots( + def generate_actin_structure_plots( trajectory: List[FrameData], box_size: np.ndarray, periodic_boundary: bool = True, @@ -417,28 +417,28 @@ def generate_actin_compression_plots( peak_asym = [] contour_length = [] total_steps = len(axis_positions) + control_pts = ReaddyPostProcessor.linear_fiber_control_points(axis_positions, 10.) for time_ix in range(total_steps): - first_polymer_trace = axis_positions[time_ix][0] perp_dist.append(CompressionAnalyzer.get_average_distance_from_end_to_end_axis( - polymer_trace=first_polymer_trace, + polymer_trace=control_pts[time_ix][0], )) bending_energy.append(1000. * CompressionAnalyzer.get_bending_energy_from_trace( - polymer_trace=first_polymer_trace, + polymer_trace=control_pts[time_ix][0], )) non_coplanarity.append(CompressionAnalyzer.get_third_component_variance( - polymer_trace=first_polymer_trace, + polymer_trace=control_pts[time_ix][0], )) peak_asym.append(CompressionAnalyzer.get_asymmetry_of_peak( - polymer_trace=first_polymer_trace, + polymer_trace=control_pts[time_ix][0], )) contour_length.append(CompressionAnalyzer.get_contour_length_from_trace( - polymer_trace=first_polymer_trace, + polymer_trace=control_pts[time_ix][0], )) times = timestep * np.arange(total_steps) plots["scatter"] += [ ScatterPlotData( title="Average Perpendicular Distance", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="distance (nm)", xtrace=times, ytraces={ @@ -450,7 +450,7 @@ def generate_actin_compression_plots( ), ScatterPlotData( title="Bending Energy", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="energy", xtrace=times, ytraces={ @@ -462,7 +462,7 @@ def generate_actin_compression_plots( ), ScatterPlotData( title="Non-coplanarity", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="3rd component variance from PCA", xtrace=times, ytraces={ @@ -474,7 +474,7 @@ def generate_actin_compression_plots( ), ScatterPlotData( title="Peak Asymmetry", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="normalized peak distance", xtrace=times, ytraces={ @@ -486,7 +486,7 @@ def generate_actin_compression_plots( ), ScatterPlotData( title="Contour Length", - xaxis_title="Time (μs)", + xaxis_title="T (ms)", yaxis_title="filament contour length (nm)", xtrace=times, ytraces={ @@ -585,14 +585,14 @@ def simularium_trajectory( path_to_readdy_h5: str, box_size: np.ndarray, total_steps: int, - time_multiplier: float, + saved_frames: float, longitudinal_bonds: bool = True, ) -> TrajectoryData: """ Get a TrajectoryData to visualize an actin trajectory in Simularium. """ data = ReaddyData( - timestep=(ActinVisualization.TIMESTEP * total_steps * time_multiplier), + timestep=1e-6 * (ActinVisualization.TIMESTEP * total_steps / saved_frames), path_to_readdy_h5=path_to_readdy_h5, meta_data=MetaData( box_size=box_size, @@ -605,19 +605,10 @@ def simularium_trajectory( scale_factor=1.0, ), display_data=ActinVisualization.ACTIN_DISPLAY_DATA(longitudinal_bonds), - time_units=UnitData("µs"), + time_units=UnitData("ms"), spatial_units=UnitData("nm"), ) - converter = ReaddyConverter(data) - traj_data = converter.filter_data( - [ - MultiplyTimeFilter( - multiplier=time_multiplier, - apply_to_plots=False, - ) - ] - ) - return traj_data + return ReaddyConverter(data)._data @staticmethod def analyze_and_visualize_trajectory( @@ -642,7 +633,7 @@ def analyze_and_visualize_trajectory( path_to_readdy_h5=output_name + ".h5", box_size=parameters["box_size"], total_steps=total_steps, - time_multiplier=1e-3, # ns to us + saved_frames=1e3, longitudinal_bonds=bool(parameters.get("longitudinal_bonds", True)), ) @@ -651,6 +642,7 @@ def analyze_and_visualize_trajectory( fiber_chain_ids = None axis_positions = None new_chain_ids = None + time_step = max(1, total_steps * 1E-3) * parameters.get("internal_timestep", 0.1) * 1E-6 # ms if visualize_normals or visualize_control_pts or visualize_edges or plot_actin_structure or plot_actin_compression: periodic_boundary = parameters.get("periodic_boundary", False) post_processor = ReaddyPostProcessor( @@ -659,7 +651,7 @@ def analyze_and_visualize_trajectory( min_time_ix=0, max_time_ix=-1, time_inc=1, - timestep=parameters.get("internal_timestep", 0.1) * 1E-3, # us + timestep=time_step, save_pickle_file=save_pickle, ).trajectory(), box_size=parameters["box_size"], @@ -692,18 +684,19 @@ def analyze_and_visualize_trajectory( # create plots plots = None + if plot_actin_compression: + print("plot actin compression metrics") + plots = ActinVisualization.generate_actin_compression_plots( + axis_positions, + time_step, + plots=plots, + ) if plot_actin_structure: print("plot actin structure metrics") - plots = ActinVisualization.generate_filament_structure_plots( + plots = ActinVisualization.generate_actin_structure_plots( post_processor.trajectory, parameters["box_size"], periodic_boundary=True, - ) - if plot_actin_compression: - print("plot actin compression metrics") - plots = ActinVisualization.generate_actin_compression_plots( - axis_positions, - parameters.get("internal_timestep", 0.1) * 1E-3, # us plots=plots, ) traj_data = ActinVisualization.add_plots_to_trajectory(traj_data, plots) From 3d6e9a52c9ba5db9330e3c34d009be137716b797 Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 21 Dec 2023 12:20:26 -0800 Subject: [PATCH 5/6] lint --- .../actin/actin_analyzer.py | 21 +-- .../visualization/actin_visualization.py | 136 ++++++++++-------- 2 files changed, 83 insertions(+), 74 deletions(-) diff --git a/simularium_readdy_models/actin/actin_analyzer.py b/simularium_readdy_models/actin/actin_analyzer.py index e997f84..59af776 100644 --- a/simularium_readdy_models/actin/actin_analyzer.py +++ b/simularium_readdy_models/actin/actin_analyzer.py @@ -1149,12 +1149,7 @@ def analyze_filament_length( return np.array(filament_length) @staticmethod - def analyze_bond_stretch( - trajectory, - box_size, - periodic_boundary, - stride=1 - ): + def analyze_bond_stretch(trajectory, box_size, periodic_boundary, stride=1): """ Get the difference in bond length from ideal for lateral and longitudinal actin bonds. @@ -1209,12 +1204,7 @@ def analyze_bond_stretch( return np.array(stretch_lat), np.array(stretch_long) @staticmethod - def analyze_angle_stretch( - trajectory, - box_size, - periodic_boundary, - stride=1 - ): + def analyze_angle_stretch(trajectory, box_size, periodic_boundary, stride=1): """ Get the difference in angles (degrees) from ideal for actin angles between: @@ -1293,12 +1283,7 @@ def analyze_angle_stretch( ) @staticmethod - def analyze_dihedral_stretch( - trajectory, - box_size, - periodic_boundary, - stride=1 - ): + def analyze_dihedral_stretch(trajectory, box_size, periodic_boundary, stride=1): """ Get the difference in dihedral angles (degrees) from ideal for actin angles between: diff --git a/simularium_readdy_models/visualization/actin_visualization.py b/simularium_readdy_models/visualization/actin_visualization.py index f797dcd..1a4b0ba 100644 --- a/simularium_readdy_models/visualization/actin_visualization.py +++ b/simularium_readdy_models/visualization/actin_visualization.py @@ -263,10 +263,10 @@ def ACTIN_DISPLAY_DATA(longitudinal_bonds: bool) -> Dict[str, DisplayData]: @staticmethod def get_bond_stretch_plot( - trajectory: List[FrameData], - box_size: np.ndarray, - periodic_boundary: bool, - stride: int, + trajectory: List[FrameData], + box_size: np.ndarray, + periodic_boundary: bool, + stride: int, times: np.ndarray, ) -> ScatterPlotData: """ @@ -296,10 +296,10 @@ def get_bond_stretch_plot( @staticmethod def get_angle_stretch_plot( - trajectory: List[FrameData], - box_size: np.ndarray, - periodic_boundary: bool, - stride: int, + trajectory: List[FrameData], + box_size: np.ndarray, + periodic_boundary: bool, + stride: int, times: np.ndarray, ) -> ScatterPlotData: """ @@ -333,10 +333,10 @@ def get_angle_stretch_plot( @staticmethod def get_dihedral_stretch_plot( - trajectory: List[FrameData], - box_size: np.ndarray, - periodic_boundary: bool, - stride: int, + trajectory: List[FrameData], + box_size: np.ndarray, + periodic_boundary: bool, + stride: int, times: np.ndarray, ) -> ScatterPlotData: """ @@ -417,23 +417,36 @@ def generate_actin_compression_plots( peak_asym = [] contour_length = [] total_steps = len(axis_positions) - control_pts = ReaddyPostProcessor.linear_fiber_control_points(axis_positions, 10.) + control_pts = ReaddyPostProcessor.linear_fiber_control_points( + axis_positions, 10.0 + ) for time_ix in range(total_steps): - perp_dist.append(CompressionAnalyzer.get_average_distance_from_end_to_end_axis( - polymer_trace=control_pts[time_ix][0], - )) - bending_energy.append(1000. * CompressionAnalyzer.get_bending_energy_from_trace( - polymer_trace=control_pts[time_ix][0], - )) - non_coplanarity.append(CompressionAnalyzer.get_third_component_variance( - polymer_trace=control_pts[time_ix][0], - )) - peak_asym.append(CompressionAnalyzer.get_asymmetry_of_peak( - polymer_trace=control_pts[time_ix][0], - )) - contour_length.append(CompressionAnalyzer.get_contour_length_from_trace( - polymer_trace=control_pts[time_ix][0], - )) + perp_dist.append( + CompressionAnalyzer.get_average_distance_from_end_to_end_axis( + polymer_trace=control_pts[time_ix][0], + ) + ) + bending_energy.append( + 1000.0 + * CompressionAnalyzer.get_bending_energy_from_trace( + polymer_trace=control_pts[time_ix][0], + ) + ) + non_coplanarity.append( + CompressionAnalyzer.get_third_component_variance( + polymer_trace=control_pts[time_ix][0], + ) + ) + peak_asym.append( + CompressionAnalyzer.get_asymmetry_of_peak( + polymer_trace=control_pts[time_ix][0], + ) + ) + contour_length.append( + CompressionAnalyzer.get_contour_length_from_trace( + polymer_trace=control_pts[time_ix][0], + ) + ) times = timestep * np.arange(total_steps) plots["scatter"] += [ ScatterPlotData( @@ -444,9 +457,9 @@ def generate_actin_compression_plots( ytraces={ "<<<": np.zeros(times.shape), ">>>": 85.0 * np.ones(times.shape), - "filament" : np.array(perp_dist), + "filament": np.array(perp_dist), }, - render_mode="lines" + render_mode="lines", ), ScatterPlotData( title="Bending Energy", @@ -456,9 +469,9 @@ def generate_actin_compression_plots( ytraces={ "<<<": np.zeros(times.shape), ">>>": 10.0 * np.ones(times.shape), - "filament" : np.array(bending_energy), + "filament": np.array(bending_energy), }, - render_mode="lines" + render_mode="lines", ), ScatterPlotData( title="Non-coplanarity", @@ -468,9 +481,9 @@ def generate_actin_compression_plots( ytraces={ "<<<": np.zeros(times.shape), ">>>": 0.03 * np.ones(times.shape), - "filament" : np.array(non_coplanarity), + "filament": np.array(non_coplanarity), }, - render_mode="lines" + render_mode="lines", ), ScatterPlotData( title="Peak Asymmetry", @@ -480,9 +493,9 @@ def generate_actin_compression_plots( ytraces={ "<<<": np.zeros(times.shape), ">>>": 0.5 * np.ones(times.shape), - "filament" : np.array(peak_asym), + "filament": np.array(peak_asym), }, - render_mode="lines" + render_mode="lines", ), ScatterPlotData( title="Contour Length", @@ -492,17 +505,17 @@ def generate_actin_compression_plots( ytraces={ "<<<": 480 * np.ones(times.shape), ">>>": 505 * np.ones(times.shape), - "filament" : np.array(contour_length), + "filament": np.array(contour_length), }, - render_mode="lines" + render_mode="lines", ), ] return plots - + @staticmethod def add_plots_to_trajectory( - trajectory: TrajectoryData, - plots: Dict[str,List[ScatterPlotData or HistogramPlotData]] + trajectory: TrajectoryData, + plots: Dict[str, List[ScatterPlotData or HistogramPlotData]], ) -> TrajectoryData: reader = ScatterPlotReader() for plot in plots["scatter"]: @@ -609,25 +622,25 @@ def simularium_trajectory( spatial_units=UnitData("nm"), ) return ReaddyConverter(data)._data - + @staticmethod def analyze_and_visualize_trajectory( output_name: str, total_steps: float, - parameters: Dict[str, Any], - save_pickle: bool = False + parameters: Dict[str, Any], + save_pickle: bool = False, ) -> TrajectoryData: """ Do all visualization, annotation, and metric calculation on an actin trajectory and save it as a simularium file. """ # get analysis parameters - plot_actin_structure = parameters.get("plot_actin_structure", False) - plot_actin_compression = parameters.get("plot_actin_compression", False) - visualize_edges = parameters.get("visualize_edges", False) - visualize_normals = parameters.get("visualize_normals", False) + plot_actin_structure = parameters.get("plot_actin_structure", False) + plot_actin_compression = parameters.get("plot_actin_compression", False) + visualize_edges = parameters.get("visualize_edges", False) + visualize_normals = parameters.get("visualize_normals", False) visualize_control_pts = parameters.get("visualize_control_pts", False) - + # convert to simularium traj_data = ActinVisualization.simularium_trajectory( path_to_readdy_h5=output_name + ".h5", @@ -636,15 +649,23 @@ def analyze_and_visualize_trajectory( saved_frames=1e3, longitudinal_bonds=bool(parameters.get("longitudinal_bonds", True)), ) - + # load different views of ReaDDy data post_processor = None fiber_chain_ids = None axis_positions = None new_chain_ids = None - time_step = max(1, total_steps * 1E-3) * parameters.get("internal_timestep", 0.1) * 1E-6 # ms - if visualize_normals or visualize_control_pts or visualize_edges or plot_actin_structure or plot_actin_compression: - periodic_boundary = parameters.get("periodic_boundary", False) + time_step = ( + max(1, total_steps * 1e-3) * parameters.get("internal_timestep", 0.1) * 1e-6 + ) # ms + if ( + visualize_normals + or visualize_control_pts + or visualize_edges + or plot_actin_structure + or plot_actin_compression + ): + periodic_boundary = parameters.get("periodic_boundary", False) post_processor = ReaddyPostProcessor( trajectory=ReaddyLoader( h5_file_path=output_name + ".h5", @@ -676,12 +697,15 @@ def analyze_and_visualize_trajectory( ], polymer_number_range=5, ) - axis_positions, new_chain_ids = post_processor.linear_fiber_axis_positions( + ( + axis_positions, + new_chain_ids, + ) = post_processor.linear_fiber_axis_positions( fiber_chain_ids=fiber_chain_ids, ideal_positions=ActinStructure.mother_positions[2:5], ideal_vector_to_axis=ActinStructure.vector_to_axis(), ) - + # create plots plots = None if plot_actin_compression: @@ -711,7 +735,7 @@ def analyze_and_visualize_trajectory( new_chain_ids, axis_positions, ) - + # save simularium file BinaryWriter.save( trajectory_data=traj_data, From 26629d5eb79981da5ff6eba31a8137d374997e31 Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 21 Dec 2023 12:29:47 -0800 Subject: [PATCH 6/6] try to fix manifest --- pyproject.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0074069..58d9f82 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -144,10 +144,9 @@ ignore = [ ".cookiecutter.yaml", "*docs/*", "*.ipynb", - "*examples/*", + "*examples/**", "*tests/**", "environment.yml", - "examples/*", ] [tool.mypy]