From 3d6545c75fb87022f647d9c4f1b319e67f063d8d Mon Sep 17 00:00:00 2001 From: Blair Lyons Date: Thu, 21 Dec 2023 12:33:04 -0800 Subject: [PATCH] Actin compression plots visualization (#26) * updating plots to use subcell-analysis for compression metrics * move visualization out of main actin script for reusability, debugged actin compression visualization * log clock time * a few fixes and tweaks * lint * try to fix manifest --- examples/actin/docker/src/actin.py | 106 +-- examples/actin/template.xlsx | Bin 210581 -> 209772 bytes pyproject.toml | 3 +- .../actin/actin_analyzer.py | 52 +- .../visualization/actin_visualization.py | 794 ++++++------------ 5 files changed, 306 insertions(+), 649 deletions(-) diff --git a/examples/actin/docker/src/actin.py b/examples/actin/docker/src/actin.py index 97b7619..1cfcad3 100644 --- a/examples/actin/docker/src/actin.py +++ b/examples/actin/docker/src/actin.py @@ -3,23 +3,17 @@ import os import argparse +import time import numpy as np 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 @@ -114,98 +108,10 @@ def report_hardware_usage(): ) -def analyze_results(parameters, save_pickle=False): - # get analysis parameters - 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 plot_actin_compression or visualize_edges: - 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, - ) - 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(), - ) - - # create plots - if plot_actin_compression: - print("plot actin compression") - traj_data.plots = ActinVisualization.generate_actin_compression_plots( - post_processor, - fiber_chain_ids, - temperature_c=parameters["temperature_C"], - ) - - # 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) + start_time = time.time() actin_simulation = ActinSimulation( parameters=parameters, record=True, @@ -217,8 +123,14 @@ 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() - 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 7cfcace7e7f33695a69ef37b2f54ff15f3f084bc..94c77a934ee8360d69145de9a547887feddf7283 100644 GIT binary patch delta 23685 zcmdqJc|27A+dtkpkx+z6Whc7~DLX@`WJ@Y7GVM}!*_YEMX(3KI$ucSNmTE|utdk_k zR#Ym>SYomc24l?jJwx^G{r-IJ`~G~s_kBMezu)&C%Q@HUTAtVQx~|vZ%yE?Qr+!l#zGx zPq`I^Z8wq~U25)riXWj-sH5gM$BNS%=^44r_Lu8xwq#PBo!=zI)zowSQw5aBEf+5( z`Kh1iCv6{rjM%J!qv{`9+}^5+F1XOhdV5R7W~=0iM{(v!OIPGWPp&^uS#)twY+$qV zj={@2TmqNJ96ERM*QJdMu&{+Rd~r3M*0q{YtB>2rgrGT?tkv?sC@d`D*gucudkh?Rqhbdn*Ra zxc--PIt49>a_6hAFWUMZ@;Dy&{LuX*(I?yV_tZ+Q65&cMJ^wUq_m0J#0`D(M`6NQ& zPY$MaKFS(O%w{;AOjRF^NIMa@d%aa%J@#UzTXVYLM}0Wie^0d`=>Gk*5=E1L;TyKA z&mZEi-sLV|lvR|?Aw3EB+_LT5rSem=g~u+x>;MGMn!HO=R1S=?vG02ZUG^}a&8)<3 zpj=Ro5C3v%!@Z_uf;*Bimh0@^ZoDg#E9JORG3F5H9S0ASFF@B-)xAzHvp!lKKE*LN zws5~+&PTX9(Z!Ao20WiDG+kf5HLAFxkfJ|0eqL?oWjXhIgsMuzG!@H8{rL0#>kUOz z!qzX%?mn9K%}~0!45q3z{kgMOVqKA?9#@9T?znG@p zjvhgRKJP9v*^JCOtPcJjM7O^YFkl)!K=$7N3Lrjx88@@vvO*2lc=cB}PwxQ^sD%ZK zfAjS`#t}*Xa@kSD>0qji@rf%5ns;8*O=bAq4Y7^gvPK{m4yvh!kp(GUJ& zFjB@}>@{-sT+tVaax(1|-^ZfezVVlc^+$3NGxk`9`rkKgGtxdHz4^!~OOwFNYeSZu zOF*5~@hXQ*f(+-}gIi!=a_WVJ^y54Tx42PzR#SSF)tE_Yzx5)f}tJm8c;ykxjZ@X$JJ}*c;V|ZJ}lpIapcPPjVf$!7jIqIbKzY>Yc@txV>a9G~32y5_6*6>nYtGfmlS9 z-FK5iE5-{9QR?Dqa*AZ#Nkuy6d~qU#CTJ9Xv080SYE5ZSU|3QRZ8rU8pc!p4ck&)i zq$0|w%pF!(g;_a~6=>x9Fm))-&LCmt1*8AtSol=&8swlJsGuyPAIn~lb8<53mFw-) ziX3exV^`XVk8Vz;yJ(f2eicEwW$uwaJR@?e%Nd=yDRppTL)o+Ak5#oBB3bV{BhUG| zFD1C*^Zb$R>p^+NkB0(3i~2q2uj+VEY&zg}jnwMqhe(Qt8txg3QKWK9z{<-4|Lu z?RM2___|@1NQ%|ycEcMYH?2nA89o(Bx1#D8UKN3>hVu;{ilkOWpK;L@B~X*=TkbVB zd~O@1P?Vg=+s&sBaB|#6q3**ok(AbSYB);+_Lo|#bEKsH+@8z6vjv0Qbua%Y{~pGT zpb|`)JhUrM0FxPy6GF!WW+LgWcOM%>u6wIoc>vB-Q;-WP}#>G z8SXK|u}`{EU8^X0xCj`(MokV`J??0BK(f}#$((KK)>Jl`H2LUCH?5kmk+F77uRXk7 zSy4W_t}Ha?_G7!gOcl8po|=kxV*`_q^15l|xQK~@4X~y3S_LO$fTzi|-pVNyD#<9G zoODVEg=YKvwnMv_eL?I~u&;^1mG1MRk6;}oWwlc!4{wD$7%;vhdvf+W z$tT#1wiAw$T!fH4B4qKh!zX7Vr-pK1-zPleq#{bZ(n1CYjFIS(frUlqhEgO+LlvmS zu^@feUu&(Li-Yf|!Bwo;Th5Zk=$0uOEyd1((>)IYrAZ}8X{~hYY0@#)QpDejjK_2tY)@5R&TW4z$-ebF!1|w8oQn+^s zhYrcrq`(E5ikf@qZ8hY`FL}?_n{x#_l`5v$A-lDPX8qMl-Lyo)u>*J}OEqi` z3I!qUq3Hb#DcobM6>O=Xxn$*-LB(c>VsIb!PEhsCN6cGQ%()6R+bYr+ z8VaMAVL)AZMJ@Z&ph1NyDjk0_(kSpDpXP=P3AjLBMLZ_DH+77H{kEqT(W}|Udc3#^ z?v__tN!d0^Tt^z8qSNF89)7D!Iimde>UbKQudKAa%})D{UxPkK)4!3rAcf;ZcJd!i=Or+)oLh?*p0IIy#v*1|S-IoMo52iOANQP2U! z+F(FMV-dO_BI*@fu7m;pysIg+YhpUOZK&m*SU!PpzlRKNgc@O#vN>~97V;HyeH@CxO1+87_g?9Lf{ERn>T zr9f=*7?BLIM|@|Y>(T6coTNNFG&>bb;`V*Sl4=@--~@FIAZ}V0ZUeAtMOMtyB>m=S#VdV9X0+~&iyV>;kGsI*dM5C) zjeG)a9~fh6_N&6%F&08^4i|g#M}=NneFG`j24aOpjuy1*UpFs%y=fcb=ZIbqc(TW* z%~8(%h^aE{u+C~x!vj%5K^}s&1ucZnY9+U=Y>jQbn(~*R`^i^RUV?JT;Bg0jv z9u};v?#Y4A>YKoK4EThL(xDBmzKlgJGQ;u zYLqK~GwbqBc+C#G!L4SV^4GHBHoz;Jw|(8pjMGZ4bEBiXiV-@gEPRG0RsNN=3(9

TLLm@{z*qsVXk%n&It?2TKoq0 z>ObK|wT;zIyhcb(qVm5*|F2rv2K;k3w}{>e)JibWe(pq1mvq}cwg5i0QEdrQJbKY7 zc?*tXjy;1|*n(I-G$H@!`hTSKpTzEDT0$-BI(JTB`_{`_b=4K_XI;Gk3+=Fduytdu z!n3US5AdSqt+%)8c`7{0y7m;0lBYh%Ao5AsR3DYLa4YPL!SeC0{rm9F`?y}6lOSPI z^LSy+6E)-e;(AJoE%&VpT&}j)^FPKppU5HPz07cd(Ln4yxSgN*UyCszjU*BfJht-B!iDFF9 z5=6d-fJ@C1Lky~auc>}MKxK0%dYHs%<(Vk;?O&)&zTHVWx*;bL;>-IdH{ch>WV9_aJTa+up|RPO<>^dmapb35><>6cWxDoB`*QZ z-Qc3jCdBSUa0>{~;Wc=hXX*>Lf`AcxEzJgi3P&>614Nx$pg6GXQXc_ITlbtB6;jaU zIIA-NA4<|U;L!2tWnL^hwcOcQ)!_V_{7LtnCvO^S7UrK6I(PYw;-wxYb}K$#b)0H4 zcxyE1Nbv3@6@~t_20-D9Mgai1K88jkTxCX#^8+`7OGv*Of`KbM`oRRi5z$KB48o+S zIWZuZ_#zxs2av94@H4!9)OravR(KR2f(Vb^NzH?88cc|~p=dC$Q?4LPEjI(Xe_sH6 zgzKN%0B*$VFTo&yp8#L3&V&7Tt^zw8;LPTA@ImYIeA-9Qo~ zc*?hfJjKL{a8UD?8=IDi)s=FM_WNXjja}eC@YkYjVs%^6{UOpRLk+OBiTTR$M``E zLK?CqCZ>Ss#%z*v^1p#%_m2U?HGy&8Y(Q?{-{=iuMd1RP4S=671E3XF8s7}m6(+i( zjmjKdgVNyVfPw2~Hvc|a;Pz8XV2^KRMCaq63kaY5W{h;606ozG1q1wuZFSEs_`T&S zU=-HAsSa;%=7THDgw=8XPA*pxEep{Ow(u~63BT*!EWlvSAp^hP<^#;a%UN20D^ZLC zHUmH(qO%a#yx5-+5w!7-J&}d}u+aiYa#5DNm)$QL2N1|*dd<%QAjPFS>6;JLiL1A(yUfKy4z^%cpJy}CmaJ;#6rM8Kd+yS)h9Ci?PZS}!x}LvWU8s8# z3p!D{&oOD_zH=tur#{{pvNeD^?Mu=pgPTjy-9^-`xf%(9M^&RTCXY2&p~ZnyHFKEN zA7ExtY%y@O+p>a(#iLzMSOXLE!IRqyW{3)Dpg9O|kk3NEWk=?V^}c?RG3!)v z+@mG$6!~lup~;rhXZ=+)c+{POacz`JVrUKY2*oS)9E)3W z!>G768>A6&DSV(c6pk2L#}mlkc1a_L8d1mYwQN-MUlVx+mXff6v)_K@r?V zi>E?7FFw*UCkOjz8aCYEdjJ6o7I3)>cwe4&v;|0gIeWj`HTRE*wjNu7*fud2aW%z6 zN??C_&B!WAvRHBuqo!$=fbFbckE4GKap6vJYUDV~ehvm=_mVgyCy2=) zK+M>&2C6T_qTxB=W&Myern$fX%eflVL~MO`7tC6t$GOg(T(2p<{dtdZj8KY826vVLD5%sG|HLh z1eEj{t%$8zdUvPv?=J0LoV0V4v|ppf5l{E`zl5u?IBEOj+d-3< zfHIdgZv}*pnT&5zqlji$8Fh?FU#{w16?r?Yw_{B-v)n+L0vaR z0!n&|m-E!`_>1f77boqQbod&g(HE+Usm7i5_DCY!GnnpZXS^}G)pD(=($Lgr8Gm7A z2g8O5<3F1m6%jxluK-qRe>vp8M``Oh##qNW)Y7@)VxmV~aDQH9T0+`+qmCKrt320q zes@h9bxo({EG1=2y#Yky_L+BHA#E3pt(e`$4(NBQXLM4I$|qRfJt6OO-E8f7X;@5c z)vh&mYnuHJzI^*U3O`AQn5_`I4O_7!k^Ee}!>jOSh52TQdrjBdQgk%g?xj0eq|^@T zD9a?*ua!k&5jc(rrJ*KZO2{;3Ht4f-j((*i9b8L{QS~PON z_Auw|_nBACu2n9FL*$ydBhemYe3U89FYHpfdBc4jq{eTT+bu$%lVM@((d!FgOO;hd zM{I>nMD~$M`Q9paLRHA?k^0@AN2m!yp_5p%S6a0R206moF~6-&2>)%hF1niecEhYx zg}u2RJbPs-b)??Fw8~M^O{nU`8lTzG3DxzMJui#$dI(9{5pFjvMc{GSRfgA;*9|{+ zispKJE(%v^$`^nXiqehDsz%3`6bhRUii#33dB{t*;BNa#k7#sxA)!guapL+c6{Vt;vKb7cV z1$&nl>~RO72uf2vFWLz1@d+P#sZ^w1SV>RPt-8_Uu%Ryb3^rEJ|8Vs21ijVY79eLe zf#t{)!xKfgZRDQ;=mNWAKTJ`%u)bV#B~MB71pMbf?y|@XrOB?SH63gSY7A?~5!8*2Z+0 zV?RM*g`HDm)9B*7a+x~v;Oa)Nm~wXc2##J#X4o-l{FImSrSju#jCj%uh@l(qDc5Kn zWBH3Ax+p^q5uWsr1FA{SuzNY`X^g3+d>U%l5nARa`{HHFa~=FJFWN)@5Us|Z=9lir zUXl7PK|@N>hP0(k`7MLcPk>=I3b1y4)3CqFJk7OuxyidaLzEu(jd~^O5zT%TUeUGV zR?}=LPf(TF)NSu}3;(6pevy#{lDQUGzW$=t^^M6-^Icbu`x^C-XD`V$mgzW|j)^KZ zp&X*^g(Mqa#V{2!Tb7$$rPo>3DauPhYF2<=3WiS{X6um zW(&)&U&=l9&?_vtXUFGY740QT(ddFDO6EWLqK@(<;-(lkS&8!~VWzv@piHFdB2-a* zU&m{v+rDh~%wg^Yy@;(pgVFvO%$G6W;%@TnMY+bvsv#Y(S;~Glp|-|XP3Vu^4`x_Y z7n_bTgfN0%_~V3x#iMRsF#R+6>hsA*|7YDoqB&sE)Dx;2dDWRMhXQ-3MbzW2$Eu^E zwekY6zv3#b*GVg(fH@;GY);8 ztU_sphHLM(N*7KA4ENRUg-?}$t$h7RBR}4T|GrL93(UTKvQ9N zc+XT5{kds%{$Q-$cnDn4jqR?NO`7jXI@gpPhIj2|n>?NksX*`N5MWkBCr9=bS|~g- z5a(5nUIN;T@&|i=Ez;yQvfmbwx3AjGrmfRPt3eXG1ru92RXa4S(p%q z99D$fzN*zzAYe_Vz5rVRN%DXBA$!ym4DN2$f_Wai^Fjb+O} z1jl~&pBxCNWC#XvOp!Y}A7b_1O9(yrITJLF*q#r+naRnpVw**8rUv)lztbpJZklOP zN~09>#T+nCi#z{P`)^UH?P|Gqry*^}_BNR$)#WEb<}6g9?PdF|?w^f9 z4V{M{^6JOvY2DPQwM#wRg-(`QiT@DJFX0LA>G&yZ#1k&eE9Kw}XN3!!SXOGY-Ksf| zU|)tXv%Cr4U2iw;6|%Jl7(wTxxajnxpGv$ygz?WBMn_a{Sg;c}JM?)vnVBC(N@onJ z`Jq_HoURz2NMcTY(H87JFxKtpzbRDG^yezH&?LfxG2HA6&##f00qCl*0Y)I_yRz&% z?RKTvy*`#bzycS;#$koOOb4_hg=VC`B&o|bcG5Y!dzmAsb7|P{IDx0a^h>>!Z)F-^ zImte3|Cv+mYs%LsUv4r}d%EwkyRJ;L!)`|hhrjY;qsi&3TUj*hrxrai6Ttqx%jnLf zbZthg2XWXpj#4Mn+%a1jGUV4EquSJWkD@B-i9$|kHMP7O`r8gf8wIhcuyo3xy6W0& zrNLkj#k%i;cvPpwJOC>_F=PxQBnKf5r%?RL*Fj6qk9u{wl3~rML3y)xQpGl$W8=w6Cle&lr`h{eJ!TLv=Ku z5;wm%0h2((+Yub@Vouj$yPzS`cr3Eqxk zV&=~s1rH-eg9&XRx`v!mm4alLwx?hu-7-u=4T21pW5XhMMgO^XQQ;-Q7wuWk^+;3c zMrmQ{vZM8gu=e;t50TiO(0si=TMOEc;wU4iEUCSc>v@D5Zze4?HG{j(3T2EjQiITP zPie$GKWELb_`TdW&X+qLN$w|*?o-xL`ce%z2X-5tKKYuXJY2U0Q|pH_hwj>~4;4La zC{h0D?BDpSwWs;0`&R@~SRMb?z?T(~N5mbtiIFQ`TeZt5?egw)`Ei%mA^oapU#xy^ zY%~k9YRjH#a5t?_OY-fU3EAljRVSnWgBOh5A7Hp5_WU|@a^p9CPWoBQnkWQ~cKatV z7pee8|7*!_sDF8=HHTt}0>TIL?02ksb3KdkAtmlpYS ze$nNz=R?%c4r3S^yZ&dnYaX$!3iCAL#P}IXlFKX&DdjVoPSO=7jCJI6X*CgnJby!p z*z+G({Tkcr3bMM(p043;GEBSd%fjtq4~&E@$WEk$bT9rA&RM&d_=SCZ? z57dM%&x8cIy=(s%O`LVp@A*m+{|H#^Bu&1E?!d++YX00ztzTm&-dKu!>LJ9^`qrmj z@@193>E1p;%x7(E?5k&#YM>k#XuOBJ{ylDC*tvg9q~GJVII>xx1bK#ey(%+bqL(Z{ zikxX)65B)bDA{m6(=aVgR{z%%is~A>uk}GzRTs^E+u3+Wy$zVty&v> z{h(#^%gmo!aAkCxR~@u}th8liWbVt@)`af*J_EYMQGu`}6Wy$PmFB(m5=TXmvzx%x zYuDJ#GNFzv=IDk2_h@dE4SQAND-t*Qcr0(2k|ENdpz8>HaLpQ}BeoEi(B?sAw0T29 zNs+s{rr)}lv=8cTW$mg&KJh3w_w1mE@}aBb{!}{2b2>-f7IbqCJ1O~nLW;APvE(06 zP}$S0_ISTDIeqtsx>kNKyd0X9sAoT7((a-dbp0CTl4x`3c@T4&E&?CBJ`*AiEC0Q> zEA43}*AA)Fg>d3q@aNVxMK&gyKC5~pd;XPkS3};PiAE#28wzzuGy0zNN{2=F3&h=% zk5U}*%d3%hUVdHm&CeZ&GS6a3?l{}?WSRBZn`1}NZnR8(PR@V%Oy;6LH7Xvlz4zS5 zE;3=A?$(K{7vM}>PVn=ON3I(Rl|5Si=cAw!x*1#E>JUfORbFI!_XX1mJN`k07`{Y& zOSeg%#NWe_=V8+(Nj`#Gl4M;er}<(YZE1d$&A*kmTtSwzn}Kmf&G%JnIt*sZ{%&4N8VHUsp*!&pPCfs zG)ZkFGnX)2P+b**@XtuKZvTmg+#DY7BNJ@&wuz5kXwS7Pzt6_vlINN$e48?Xo8 z^HOBQ{SFv?@sJ2HuvEE-|w^42Rzr@seMz1}d4P#(AV}()Bg=!4kEv^<;ZpX~R?`A1K``$i! z^u5y+GiO$$psV}$xSZ*%tEc7rvd-;H(%9i}a_hQWrKzm58b2TU7yYx0TnSJaCa^Sm7X<5!G}*P!gEQKsH@FXGo3 zICk$T>ILcwDLs7r3-WFwut7@kiF>?V=f-jrv*xV8u=J4i zjp=}hz6xev$!YS^8)x3wK_wp3PK4B&C<2Y(!Rp$o6r77FX9N!IdPd%H&V3V zbrWLYiD@f13U%{%7h&q$d0S_Sk}0>Y$8#rwuHY`Ugh-GjgekJAbX~xXBzvhMJ=cwl@opv@G=w zc@j!qTGLZ5GI_d0r!&M~qc?XU;%^Kr*FN67X`jSgmRB?(rQX0*Jmze(d(n~Fin>6Q zPxR2cS&}%#k5^hUeswdYf5tw_FP`IV2dkZsv+DuO>i?%`{;wtJWSrVaALnfpy)I~= zLP~8QW_sXak;&dN>cu1D?gHV)2Ubx_4R4Nr_^+#h{aZDV`+lp2!D>&Z9hm5CNVg84 zF8B9d(L(!4otmZ}?Hcp1G3Wo8wj`z7#|;di7$$GnA8wO_zM%FHo9FJbz&7R!Jt`+) zQ5xw*pQqTJFXm2}c~ZenYQox)n-KpufDAFr3=kGZUf2LTL9QFo(hcwemxqQ;_|m64 zrb>^HqwGr&=zRv~3bM}*Z~`UKH|?GSR!ElxAc#0!1dNFC-+}Ldm};R}gcmn~OsjcH zEfpoPAVl~ipbwlvrec6(F<&l?#HDg4AWjX$#yN4;5D)Qz+W_QF5pWqP8wSJ>YAzrw z)RE|J&hcY%m*^q+Zvf?+zM!xu{A|G}8Hqg|zPG99b=?U>Y!PUUC7I8Sq74#A`<~T>uxLhS-n{_<#Tzp*;n50!2jMGiXo?CnxAgcP@H?Ky-`;K!9k` z2Y3Q}#^tPDWC#RZh@Yr{H^}FR;g%sfPf$ISc|9Q1Ij?u_cQ4Sbs*c2MD~}3ckA4o3 zYZ>jsJ)wR7gedeF`ViJ{HkVEZb0?h8#)+EQn4YfMwuf)Yd1!O^ z_)HZV9>KlAmEeJ3f&e9%#Gbx=6F!kRPS~IE9U-SvLkEhnazyRdXn>1If3q%Aw{HQR zT4!Zk{?MbQ!h}5D90-|bSJ@$3-EX@$B4RrM3zSKI-hdb=@lRkz(%t}@WLp1Z`F|BX zH<*`yL?j4}zN7^(l<_?Uw39HR1CviOWSFCoHFiyG#w`)d0mU%Sa3fyUEs#&;fEn@K zeZUq?_z~dp;&gWoyT5|TM5rJrxl}R%+fMHoAN&|O6`qY)oI%qfqwfA#nfjXEM{Wh5 zLpa&yfn#PS(;bX>B!fxQm2oezo@ky7Ky)_)7N9wv=*$E}LAf&4xYKW4gVd72NuT!O zE%1)$u*V<=Cjd<)rRut8-`Hdl24c1BTR^1N9pM3o<|L#OMed(4)rhEL4}1kwi~qt< z?ydhT+S`KM86x_U)O8@)4zO8^;t>9qYU1d&71k7?9js|ODkO^s3?8kbJ<0uRtUyp!oLtm`q=XJ=BD!0F zHGC1XHFCtv6M!9{&iJ=%Z$xXE1v!pNyAtH>2;X<|Iauu7RqM6DHvTg`igI$%ZBk~? zYHFt9v~AAgQZa63p`T^+A_DQnVo(` zq(&lo`z88**t{I^vM}fkO0?8r|1LML=10#~sH%fW?QJ2g)c7d6?ua|bCSn%zwJ>ay zapdcCNs4>pBqDke&|b+TK)4Lnv2|zK=zY4p?iNL6|I?B81(1hIQARTU+B4azBuauO z+V7$n__+xvO8jeOBjhXS^h6u^rsc}eC+-ZHOG96d;nGPQCKnyLMBQ{CVLthMX!7R< z)@p$ta}i4#B<_JJQw0eU%(#VCm-**kG@rT z3|J%Z)FWYJhVG1vTLD%eEecUM2b7@mkjR^d5(_zW1d$5@%!!L{0MQ`w*ag^3^v^*j zD*D2pOb)P<=xoKiUf~XG23m<{z@u7n^MTR?q%9hdUWF<_%zQNWd8{`bJ8 zRO>T*JI73$3N_5`aAORDKO2=C4|EPIH9p3jc>6Ax;W~z#eFiB0+IV~cl*m3PKRBc9 z=!+KJ$BI-*P1BPF?DOQpAN!Ck;kr0vyEM2OsmKJxgnl;o=>16e7a$8HT_YYULVq42b`n^N_{m zKXm(&)lwH8JX}hXL3ap$h3DtoK>zAc(eTO+hh$y`qy%|FT=PNkTBP7GunCO9A(fYb zKsDaP=dym4S~PxS;7A+VXg!$ZE&tR3A?b(zqJ1$?we8=H4&?s8brE%rE(oa;Z~~Y_ z@*8xYxuU90jvMYvK0rbTX3f(G2KFn=y@#Cs1ei%v(OV^dsRTK6Z(-zaY+{Z_&>?E_ z*M0Rr(i^J(H|TwDq1=DPZ!Ee+`(NNUI{W@x{Qj5M{|Eg3FV_FbZv^GHhb{w3Yd4BN z=8)|EK6hE#g2ED(5|AWxkIHoHTVe>E9lBV2uYTJ0{^OEE)aOt9Vw}c^-4B2l^P6Wx z>j@x6w59;ZRS}InfR$Mdx9M&9@V?-i9;+%=d=oyiN+(~*;_@G~%Af$5h&KhB*4VD~ z%lbA-t6vjw#i)}!VE=%5@?oL^?Zl(6I3;4`al8=TBC25Rk{Ji_%`Xekxx1J(8zkMNBxLn)p0lO6WH9^nQp%oOwb zN&jr8{g;s=wWVY1iK%@Hn#dTUC;dshH*?gq#*WZ?zp%Fo9sTO=w?jR9?#Cf2b3Y)2 zRiHYOE(ThQaUbDW<8(52m;lWTIvL;+h-5_&eQwBI4%$E{x%;Z4=rju2A@#l8kC%q% zwXsWsJ41XE(&-(EvK>)llU)3vY*jt8^now-O&+ts*^+wYX#*ea2_E5kEg^d4=>r}1 zRF7GmmYNiEZ*}*0W$H_M9i|lw^&4N48 z2Y}P`N?L`9qghH$M&qPcHbpPWq~!{4GjKAyRj)i{pp^Fr=aWs>i@G(?YmYBlj;Tos z>1Ou1o&roEQad{MjEhD}Y0{YT(aHl~%TC9AyFV~)UpBcM<50F?ycube2Bnc@vY?n~ zB<|w5JeLVF3T2>!SJzm5N87ZT7%9^0lz^sJX6-JiH z$*6bv_sH4s59EySL-gmLVoW74X4Q5QJX8@+MHJHvB~T2xCI<=uen^5esDhj^2l$Ci zir}3;nc46kGP7iknPt>?7pMN*96htt0fs%faq@6BqE~)IUW$Z%9R^HD}5&9 zb8+|ArPp4+JFz7FaoxQHqX9d!a0i#Rd#+#Q%x-T<=zpIxH1YjN%gl}hcjlYv?>DD; z*C@%7vJhkgUSmp(7|yU&if7uH?e7gK{1!VSg@ z5#mQ;G)FMf4Z=<9Fs2;u+57I#Iq9Xp1V=Dsj9N$u$8?;HY~@h<+n$#D6-5~AuM5? ztXHBqT)M+@Zg@DlnKIeWSX3 zcIdm==!r+uO?$7M;C`V-Ia_J6rH)ui1n57$xfKf?ct-zpRx4`h8Sg^KSh|zCKlgNz zrhQH#I<7*BhMvw6&NuZUIaa~O;X%edwskAxBTQxJ-Sn}gOW$q|27TjzrA`3}$2yohcqJ}oU?vn1E5e<7nH_pfar6(i| z{l6T|e`+LMlYWitDW$NNF!H`v(`3@H5Q4vZtu3A&geNwQ)lV=!B^V<5J%WemCz_{A z)WQ)vXEtQwUR=oEsPEoo|KUQ7joBUBHh@GWRXyzfX1qQ+oy8(P;?J_s1mFP)$i^qi4c3QFPO13%J#c%J!2RM%bIM{-a9OE zS1Y)#Wwqg6jnRbr5$ioZZnw(zfAAgWvchbeq-akZ2^^Q699!`)MI`A-h$(Cg?i~5{T_-zpCNV|F$M~WOG zsYRd=HgK{UcXiIz)_k=HSOGx_ch#EV=?GgO$$f5^QWAJ#?~& zIepXi;Q87L^rh2#*9}HPp6?=N0l;Y9r8B%E57xy!BjgiztN ztT6A}Drbxj-mcTw=4z+y5Si|^X$UyRuVvGeFLK_iYVih9S0S9q?Hjf@pVPIg417%; zwkZ!)W{1D}w7SF1#>#SS$7lZO-kf`Jt5=8Zlw9=AvV@qqDdo4-Hc)Fp1ITi~7PsH%y?nxt?A zquC<}j;j&L?Hwt_vNv`kEvGnrlXr8Fw-C&s;5)hK(N&BbEQ^C+RAs|@S{gWSy00k_ znM##)iDamFaJR&n{C7$jG10`O$e4X!Hs4zNnG+wE5yV)vCC4L z8fpTuIVSA!bPlyx?n>`C<5h>tegi6bzk4it2$lONH2V67y0slI9;cS(1qGJCK{U?2 za#~qhSh(x#d-Cauu#KU`;lBjsd4i&@L;a_7q70_f`^WoU<=#$B-~`he@+#{-25v-1 zx}d7!$M2ptC3~8po&!K@yTNPwuZFP zr@sB?H=;rx>K9gEqi83`0+9Iipq@AZ{kXBAjoT^yaC*HjaJ}WJ&y8nx0&^}?p&JT5C z&(#NEY(l6JZE{P{en`ZCT(TyOP%-_j&Et|DdTNuLQjjJy!%?&6W_bAOwFYM0%#OMg z7pIkN-ThLM=V{fgwNqz8B}#LXUg<`|5ME$x&khkLClez%nM* z^w`6hDsl*6BPb~6uK(V_ne&Y^S#{Fs)cagDR^=!!08QMf31|IRDd?FFr>0r$Tdx z&!MP{I_#jpM}2Ow5l6%4S?+@IJ>DmSKw- z0sWeNofX`HvD2F}JJRZQm4DUF7%V3a`}ib~XfoVWq*u9}sV(~m9xY^&gCD0OjU0j~ zSc1wH)Tx!+{n{Lk(vk#RW8|YwsgDsQS2ndcWnxLxPM^W3cTNh_&Iqg=M9FJ!jgF(F zntx=9wyPL&ecXqtTlmwEN()fYq&1>EN0=B{;zxCoM317{S5ACDKP~w&PYE|@2A#|g zpmv6mqd4MzOo~)A^P%kND3p|h#NVT&Z$?PPo6cJ zy!6(HmjQLET6E5ZzPc2`M}~6{>Vq-1`~%|<*K?|hG6Fqa{(0GWIZ4uft25b@}pP@r-@X=ntk(Z~ACnwlmLZffIh`RH8(k zmbpx1-*l5fWH^Z};3JU+sS*&J5x7o~qPlt3I+Wx2$s}^$=L`nfmyN!1+~M76DxQ{5 zmOK?emP8(-a_9jkl~_<+`G9grYsC7=(|K$-Qvlk;>8qz>2>#KO zv%H$qL_u+Bjc^aFp==c8q@nMxQ?f^++o7^>0=8+k1%amdD;asAk3^dG&d09YF_lr+ zK5kTZ9^1mHdU&eE0(pW>Y&z$%>gwYYei=yV=|gd zUFo0Xz^GDYKnpZ8uNUDOgPQW~sy zSo`q+wf40P-ph1`)`fC_wYU3PQ&0t+Q;1pIh^syo6skB+wzVkP1|QR{dq%IRigTXQ zBl`!aw^-k~+2n3s>{G4TV`b<^+)!ExIvk*eRi^t6CLm|GfEt3eQ$-9DZkzP{uL8AG z!P6$g!@PUyOI>SYhen2-h+jFW*7Ns~ccvg8pJ>JQIncEP_LpF1?y?pm7Mnoy$TTLOF-;9=psZ@p7R ztvUNgiqJ!HOhGL~(h@wxKjcl2D?nagK^@>QGLC**_-a~c&fFr@Tm<_P^41*GL5LP$ zys&F+4l0hi>rI%84GpHK)#yD;SmxR)CCUBYkK)ZllINB;_|@*-6UUg z7OTz);^4bl$Ln2Fn^rpmO>of_L-0Qy0ecG^9_V2bj6Azagfm zexcAn7k0GRXD4hi(qROuEnzBlwmYE*pA|JbX(Nd!LxnCG5SyuzXAexmh?o@$uVxRY z;lO06jcK_T@hlG;4M8zqpF8LwY!MV>nH8ud*f4*RT3HA$Vokwf?`YEE;al{1eZrRv zLh4a{4HKPxh`AZ4hHS6~4+&1&xE1Z~JW_aJ&J9FvGswS?TCjIcKff5_olHlXtWjl; zFU8jB*X1B1L|q$jJAlkuqbKUD7-u|bm>4)#iWpmi7?5@KzxADub+thpS%ShPa`|}} z^yLJDRnl?N*7!QI%??ysTsK{xN?*0%*GXZb$E;S?9dn^HM0lnDJSj~0iNck>Co3lV z&5mzJWVV6Ui!8!vvcl-O5g4qn7kadmmq~e8!uZ@-WmW?Ed7Wi7kLS)Q+q77FzngF8 zkpSKyWmuu5=%HQL3Jj_Z8tx>I}zt+87sBj!Z-=O zzFtA2Be6Q*_7JC+f}TNI6Jmhv7`)ejo5ZU5-NSl9Oz2owKYF&;mR~1~t@(Mv7?I!z zy8WLqu@zMV4|e!6-PxPTkhlDx_35?x_5J3sgTS`S`7rS_P2Um3>f{ABtQ!AMgLE#8juv@)se0L;+~rE46IVq9 zzp8is*EHwBIbq3{GVC{{IJRw8>=9+$u3gKNqdHl^mPPv$=^T|uVbHY}w zJ01cCma{7d{-%>D6wtJR&oc@^ocA16Qb)t>IKmC1EG@CYD&waEfUS@ll; zFG_lTs+8fnsH^?SS*j|027(VbCa^FX1uxKge(nGC@?g*XrOM^IR$1RJZC0LcSk7d! zolTtiz9h@T$o!4lUm7s0GuC@&A2twZdEfeveerHx`vpr5l_}_2ss^m^_Re>@WwGcQ zv*+5!_3<|g?|vy1G5-Db&5D@$DdM#arB184Qa!tr0(^=W-r}+~y`@x3XIOCpJQ@bDzxl;9aGt_wL;=ExF1_|K#5C~u^ zGC5X6YWg=5W&vr$5hduV^rkDCGHc5q28tLML=a{GS*z8iXP7dp$$)#6{}>ofC!ia# zS8MtzQ)U;Y2b$AAm@?Zit|7W)#$GScE*3U#WJghV7|EM@GI zB$cfQm1XSH5QD*({dY#|=ktDlzW4Y0yC3)cxPSLQ$75!$*Y$c{&+B`pt{ zC*D!Cw4OkbzTBLC=xM&H9v})=6*JXx(ByzyV=f6zPOU4JDFG7)lWL>d8dtD=Q&`9P zPi+~Fk~&nNP;Pjw!imSK8L)Z`-@1ub6ZvRr`Ot7jR9=JK{$1f3%7pz@?+!n{vUKrd z_b1BLCsXZSZ`bl!EVjjImrI)d<_8Ux@x$vIKW`rNtUg^JzH9qV(aFHK%AINZSSCJ; z_MtudsQITwj#i5sdCMq->`QjSw(b6KON4}x3r&imEZ}QKYJYQxdI+VhOxDX9_~In8 zuIq-yc#7Da__PH9nTE!)D+n*Il=~%Cnd`0RYG(JSotcrbilSI8SJ%(HB^JfBc^2HL zzx4oo#^6)c{t`9MyToVt336}b_a?ek&ez+$E-JMnH%+o-L&0t7gIb%mFS|gT_N&;l zI*>$(mZ)t2p?OqJwp1K*?1C;oREOclvJ>W^@d+Z_d<>>Q7$P*_4CR*wi%9oAkQT zl{r%=ct3fD{G_Td-`uD01N~j@7Gh|BpoURH8i}8D;h2w4-+r}_BkvUE_xe6d;lzwn z<-TryTL7sO?DQK}CX zdgg2XCUJ|%iSb~?dD^KcH%-R!HIm%V9*bGy%amU&&%cB2tIME79KF%7EC+kYjm%7WkFKr-qby1Q} zg%*?%@O4gJ{1!*urSbFIOP+P}p@=D2XxD;8u(F}(f`W{B4<#f|zQ$~Q$q^AD73@_w z2@>3++!6z$2n@XjS*zs9Zg<>5HlHBl#+TlYj(8AXkK_O{ZO z=`c4)>UAIsrbrb2LjFV?E`m>^3-}Z)mxm@A5F=)I*a%iTolU}$BLqXyWd0Vx^mvP` zz(0gjOC<2y1oH_TUK@q49X2yQmAxi~z@EZ4VsYO1vVHv0NpeKQ49!6>#c)8LAHLjZ zHcs=ikSp`)pJHvB3Fvpgp_E|(8q|Pdm87v#c3EYd@v71TSS~w-Ea;;MyxHvn7O5t} zrH(o=TAVTr@h5BrVWTy@Gu(k!i8X}=va+`ijPYuFni(mjkp2=KwV}cK&XQt_wDg(zSLXSazW1g9GUZj>t4+7d6jb%*ntI4^0Q;7&yuVPVyznU{7Br`N<eFa^m!K;xBXez=>oXC&1EpNWLLf!Cf-mWD++`p+Y6}6n>dPJt zy9Q@ING`1LwhW&@l^P(q#7cc`=`KA`Ma`z*Cs=DCb!9E}wk6~C=H~i0qLaB2`L)dn zu!Uwg{pd{CxZ@6s%dd?gs7Fm>c@E7K-IbM*`EX(n%2pKt%vhNq57JUh~z4Dpi-MHV5XE?b+c_OHIvi;EShbu8pGy zxikkETv9?EUJfl#*VJkYa5kjxZEjY_&lJYfffOR{P-quyz<(GgR0J_1G-P7e1|(ul z<(+NU8TOmZM;=n5*pYeNw7MAYvy|e+E@8>n;16aXN;7_bK5AfZ^DsxhFD`BZ{jg{F zY7MI$mz131<5_QTsZl8fcW>5m`K;x=W?tcVTgPj(EZ_M#bqiWSKI}Vz_KRQ7!q-@V zm|PlfQxTe|o+CzAM6|PO;yz$|w1#Sr@dY|BskpC>2M@maZmp-%%2S+Z!zrE_x!OpM zW-8M0HERVEO)2t{_HiNO+h)T09S=D)hV#qK?!mrhh(Zk1=SaeVL|b!_7PCDS$I<1ovQ_k0euD=uQx^4@zYanEFb)Q825kgwKq#Z23-4)1jlVIKq6Mev;__Q~AZxVd8yHZ(^%SG@XMiNs*aOaprf#mJfG+RQ@UV}kLN5y46jW$v z0+6EyB(?W#DdI6M&eXGju-pNFtd~QQD?dqSW7-d4}NpWK|sdOb#*(+E9Exd_cfe!$5$>N@>foq5{6}!cTiC zr3kQ2GF?`E2XQ;8F4W$3R#s}k+AP(glwdj2oZP+ z5g1~Nz9R^#Ll7hmNX#OTQu1v-!2%MsG=DFVvg{~46!GGOj>n9MJ z|Gs(3-$U2=GX}|jh~ee`3vT3pNMxX=@WVuHcxKBDg?LzV8>k_p<(bags#xn8IyeBF z*M?oKK&-UNo1?{^LuF>J!AsQ-h@&OF)}j8^EwKXT+~)4-%yB2femdzGnsR%ZnvLksTDoD@5a&bJ#dc=h>mzN*v? zycTDUS1_War3~)JD_lRid$X<0k`+3O(r@g&etP$o(>AJEI`h-7ue*MtY4c$lL#hqNl3T44$wO#`qOUq5bb-Do9D^pL2=~(lhqgg^21??>&>JA~I$IzK zaxk(5Orafywqk)frCM;e1GrTT31=9hnAieskinX5qWK4NAobOD02T6BZ7qu+99CB9 zsNpn78|W5A&JwKngJ8r(&#i&y$bM+`9~_DcIs6}NU<>?4W-WW)f0ehFlMjlHK^NWc z1|iq=y0UG>qn8fcHmcx;S495Iz<&k#Pmh0#65AYcL-WfCl21vWoKr0Y>F(U<%mEnzkpeaapCQ%-O zmLQxv1d2h;o3{XQ1>+r7AUyO8l>WY~hBV&=+ycVegTVPvg}FItXBxN_TwuBCfWW`7 zn@=f?5f5tO*uV_|pple7P~gOitEEh(lYo!lMiA^$g|v2C0+@nSCrc1?ntcTS)e4>x zM?+Vtv!Tkvi^ZIA&NlEaWzh|!1@lG!0ssN*?FL$j6NHJT$x84r2pR2OOFH)oOab7rAn+l~jR)t0ycYO_Gl&6d zNMsTcDJU?9ecV7xaY4xs_b6}$8)yGb6(c=xIPE0dkT>=`>46T~M*UnQ->5(5r47&D zIwv63`Sa36VNj)|`U0gqCte^DIl6myk^mlP z4cZwy`3TBNoqQymaq~+O0x&eMCSH294N5!w-sGt!str!t+#Y{-(u_8o?hUzUm)o-~ zXf^318>|lR?j$GbFSq(Y32ZIaiQn~;iXo`vge7nQRs}_~#V6Nf9V~-39bPQjG?fk` zl&>xDaoYw!P8X%k|L<&wDd1-yQW=hxhS)x40Ik3`T}WLJ6ya|P3j7NIcS~>{LDis$ z(B*7UbUud&wMSt990>ts1p;))`jn+Ow|#)lE7(-E1%zVTB%w?G=AfVmiaqtXU}ISV z;IOc^l-!`$A*d)l=%eH8@V`lpfYUfJiiXKdDDTpm{7rShjjXKz2ig*81)L=jM`jal z47psioI@z$mJ^IfJpSpPteA6G->}7B9D1w8Kn*^-9V{h6RHb( znSl`(hC(N4fTGArAVN}93aXAX2S!NacaTV+9i>vD>@hFs+C9wg=uDv0E0%KvzoNr} zqRVb)L#{d2bB0g`1iDUv&+bhi&>W{>lLe1Hh;^dz)+YUnZ zS!O^NbnR3&lC6h;@S%9l&zB7VajH#q?>=bUVlkMb1NrPm!#;<>8t5G2lJH|uA^_=M zlZKY0V8mI!m0({gDa!U)4?VeJ3KB|xs38#m2{Mp7Yzd$c-E9SBkFthQ?$V z$YMKGm0`3&z;j2Q6$Z4ZR^WeB!ts#$p{16%tb#St*%|7O=B|v=tw0WnlA&^finLgo zi_r=9>TX8Ck-xTqV5bw*_Ynh-3o3098BY))=1VI8F0iOVx{1nuQ1XWX(hdvqURVmr z>L|BFYS{2OqK$b3)K8sztAOZsr`UGUPpu3u#`yBssy43R3KLQC?jerzF7`+AVIjt7Et!p#3(^nC80ri z-S#C3P<*wK6i-079VF<+gkr12hQ73<;P+3)!6)ybto- zw-_je7WV4F@qnnMuwJFWSxZyU?_USmtmp4@Yzrv**T4{D&#{syX4ekT@w*GIXCc}* zhGIso#EO266Hp)n`Oqu@!K@y(Pc;4plfB6Z^3n;v6JSJ(XsWSE(Kg@CHlx+16QPP7 zv#mJ=QqwWt!Cp^erc-Uukc&>a^Mbr5#iQNstW0F?ky0;%+d)G-Vf}1hkcCFII~~}@ zv}UB~r|5(yOa?ZQY$!UsJCpuRed(DtWNltD(*4~?#wS5n5DKN4k{0D6odFiaBYBzg z&HU%wLOPT*6hbotk`UPA&(Z1Jj$cbV+thAD*Dg<) zJTjZ9&K)HwlOav8jVzL)@TAFgwDxp5l&n*}o$b5nRsLYf`Pzx>T1s~Irs`g}G72;- zlmy<*&WjO9if;w@7wqCAWr@0V2e1U_&$|HRK^&|l3RnSEBoE{d2s_RPtRx8~-13s8 zdskJEfNanV6kT(ha0^DB$2$K89tS^?S@TSsyfD}|MWP&Gwj&6 z7)Q|=v4zj{-hAQUZ@a}3HyOfrm#$0R`+DWzIMLZHr|`wQujz#28|oM_+YM7~VpbPW zUGF@WQyHD-t#FbObuhkeIf113yghA2gv`C8(+{&GIvy?6xV>459CKfDSJ-rJ)9d_) zugXXPb;0@1XHs4y=tR1vq<+P^K`p1SjeBCBoX*hrY|-{fcf4UVcPfJ~c8e@dSRw6qo^3c8PEZEsAMXTkS`!w+~+8H$I&&x9!_`Qf%9{ z4F6Z2Zv-j9gglH)h(k^JnNLhmV#oARlynG+s(R-6aT%{`Ty0yj)x+F_T6>>d=f5$w z*C}Vc^|09NXl-hP?FIf4c^6;KYk1eGY8+b?CvHR2+}Ue&RKs`2&AHx&DlUbd zcAqYxb~M+1SdIE{B(7Aal|@Q^v}f{Fx{>Pmn-47KV9D*r_=aClA|fIJ!R+N(XS<*B z=Y)(N$E-*G?Y;KO*3Xw)fh(3>=Y$U?)SSJ_aqUobR&!|VY$QmyfyJ%N8#n*DG4v`{ zf}s8QK^ORzaPB>t%;$L%IMmTb0hf^y!6FJ~$XH$*nHPq9gihvDQv~DLdG@nzS_Su7p%vFsTf0neyK z@uv(kar`kju)ojrfxslr$r)F(B_&y<`YAW^#vjW#YN|iTzvfO><+t}CKplvEL z@TdW}lprxC@!U^lke%f^8X*`RMfRQ`MDVBtF1?1#;|SQhKR^qeJCpQGyJ(b%2qnd% z3-#scWU;eTgV)wi(*~l0GNT(iSB;$=K5=gNsdBV+borKEAD?ZiYoa?+eRNe$2TFJ4NJQNLnE zxgDN)GB)T+T8KLjEq`VmNwMBecC5>NWdMHL6HY1C&OEUt_%ie8sB9q*00I4>Z3d?M|S6So7_eAPCKOj z*XtiUtL%=h(+JX4v5kgSE?(d*#DJaofq~`vVU>L2Rz}nLQpS1eN&CjLu_w-{Fa7WE zbXS(5f447*tGq7z#FCJVyRceE63Xm=Kfc6&y|h@1v6<;NRk`T%kgw9xwZPUE6_PcD z_I$AR=@Vs61hdbSvE3MLQPQr~xqoL3nbWt&3-ib_D!8?FL&$}18j}4YH*}S-gFAcg z$13YRS%YSO9y)y}f|h+ew`@|JTB7cFd)@)Qb0a>c*#S2_<`o)jH#C{-Cg6>8f1I=bBv`;d0q<9Laj}vrW#)qK8ETbA2~s)_XyM{ z$=psE-g&TXjlX}yuFzL?YBx5@Ku2^JuanWr87>-U+Xw2`78K=R2-Sx##S*ZMqlGyY z76t7m3QKO+zAOnNGegQ(as9ro4Xx9R7liELXYZSY0Dut zerwg)^U%?`*6oOM%mt#eLEP}Qm&B_Ukn5l@;tfc|C9$Qp6=%!Y^LEypo4HPZXCyxS zz839Nd(GwUXHVwmro#r(%@fbRQK$09F4nqF-?34bfUYlFyenERGqB^N1;cF|`DU(- z%pXDFx}J}{s~XI9b66UF@mbK)Femp`#vo<*g$Z|fh#zD$V?Zb#CgJih*qYnjesQZe z;U*$nh*I*C3Far{-T<(M3AorK^_8R=GmV0+bQ@l(^;e3>?Swr$jZnShQGz@1T1t*f zNabi)rDC^V%<7G}#RO!z=%356sjmRFpr;h<3%B#i;(wv!ZYOBH7^a{of%4pkZ5zv? zhi@r;i=`Vh%NEnm)9H$mqLTv*>eS=a8&f_X(>|vH^=K|$qhhPl?d!?mi8|e`+!e{s zxtQ%XxQG}@EOQNtxHPiJr`U4weR>|2c;F(yu2Wn4T|7pEJ0FM;amH#z+|n25pyv=s zka%V4j~6ei*1KPRJOV!=Udj-siX4!n-2JQXLZJ5C1^LlGa<^g2Wzn!6;~pcinOAS4 z$cer91QYR@o(8g5yWL=Bn@cVRmn&Xm=Dt{62GyGVQY8FDT*6PitwZz{iB_@TUJ^+$ zjvfoIomx{MHq%9rsVEkx+LyY7SC58O-&VKRZB1~kmG~`oeI)f_IOEir`H~a4V0-_t z+Xaxs)6((|GfhPN-igR(I@dP4R`##cqlNld@gsP>+FmAa^n9}z$7HrP&R z2%iniG19o^p!-V^-`H~r=9=G$#0ZIOrY^=sB+-n--t3t#=p@Cp)z-;%w^Q05(5Lx9 zGlqRgJ{lTx#s6xK4b-0ZzHx(^VZl#YRwJ}l28DFbzFNYK)%D38K1x z0$&}fF_oY5RYy7Q&gU&Z6=Z-&fFEBFdnP&1iGvt0t{iwV_M}hQWzBzHv-029jARe~ zp#86cSzr!-0@I%b=C0aiHe57~ygryEYr^V{*w4Yg%X*cJmekl_ z^u55a;%r`WSxVrOwG_U2_Zjr3B($+ak!|j27p=JqCg%PQ<0mBLpJ9wpL@TDFI$2Tr z8}A{@k3`N`O4t2tU8|97iV)r$s+r$tn!|f?EQvKu2Y%wQ7vg!*?Z3Ao>}%8VKYG3x@URQZZcal>^ag6zx}>9bVnLJ?oN{Duoaucp}A&~2az@{Q$UJ(HWpFHnAR zepOJsGFI}CyUVzFw_cq6+1Xs!yg_pD7qs9Prd8iD*}R$Wu-R+ulrzy>@{l*Xp0aMV zLCFJgd^|}~^p(xrKkyiQ$D`SQNk_vd0@=T;H}u6!n?(iNfY-@d@=usn{vFevu`>>v zXEg=IBbr*Eqw$A2zSD}}PrlGGgFt3nL3rXgqW(FAm4Al7MG%Kn!`*HC)|!6?v%%!1 z?u&MWp^1T=&j0b%`VhL{X)1}$h4CE78h0}nq&qRPc4rR9ITdI%)XS+xuKvs5YN^( z2s_(?KeQNb{ALJ!A7^M8t=Qq{djTC2>3yvZhz?C+zb5uLo!9MreH0@IIOj@x;`OH_Z(10_;lQ$h5XD&RR%Fopg+q~%4 zyfaF~IJP`lxsXjJfrp!gLaNzgMMC!Pmtx>ap zWOk`&>eyd_gT;lnF=lP8tPfA0Ht+HB?D2otnmMoiCDr!GkJ(_LV%yOwYjUJ$K(b1o z?@Dw-mbl;pNL^U&_AdE0fPu4aD zKI`a7AQIZ2J6dOEc9_15Hhl_<6L-b`?7Vsi%nwjD??B+I_57#Aar*59CFVTJYT*GQ zBljsT`emE8$MEk?tGDHHBHqtY@HnF=fy}7qcCKM`oN2a;8!bFTip>o`>|@9D=KgCL zshSrD5cXMjJg2a7cwPx1Iccz3r)Qg&>Jf5<%th56BgxgjI{B(Cmoa!hGXHUGksX=g z%)ZFD|MAWa&gmi(TG@YOU;?Q5A;}}W|rTiI-KiP|ze?ZKuXU)z3-JIEr zm}AB;lyivrz;E*xv*xvGzq}*L7kckT%sHbNN)RQ<=u4yE6yXP-_73EcTe%!W#Z22z z;Xh9it5s}2er(qXq7Y9F)vETIGn44WZ3*rb2k_$`I;X!Dr!n>d9!0el?;OQ~VmksmRNN z#jy>JV>%nezOz&)jN@D3?=0&eSzym+~)uRtWi_l(k3Hc1}A|JZC{S4Td6x z7kB~LM&ay*&JBNbdGz-fNxa*oBiTkX=kyQbI5iewhv|7|2W9Dn(WINXT?icH`Wud! zjFvsEcxt+qZjG-^nr=}JJFlw$eB^6?Gd^aI1)1A}L*7`*5)M(uo*_s4p~~F|#^K7{ zzcC`St0t z*&)|szfO4QXt_2<#|^gLl+W7mi`e0BogpkPw|!Ic-8}lcp!T}`a;+Z`P{{QUQ2ME_ zw9H8vd|R829SA7pk6q^DfY$!EncjCkd#-h~Ik*?u2D$#;#vob_p6qzQ*3&Z2ORAV4 zJu1g|XOh}9Z%!ioo=thG`7H~H1iAi>1nJavB_0mE(aYxR`MRG2b=i_At)!sOjje_4 zqew=QrjU$WNLv4g%!G#<8+hai#r%1}w0uF}*MS=ip(9GRri~t2M+oH>hYre3{~Nhl z&V}R%=12Y)N`H~e)68QL#uLKN2imI~v7ZhjQ6SgfQRqg?O)>cSPp9D8SjXXJa(i61 z?rr(-Wh2{|MJ-h;Hz|~TTe#)tc-KTyeWiB8dwRGs*Z0;VDVzkoHFbo-H19uhOO{$< zacpp5el&9T))M570-5f7t%*?F*9a#qsTW2irGdKhUbw3vM>r9WOTH=QP!DQ}7808X zJaUf0v7hc+xl4pvj@>i6b5EQ2lOTk;**7zq$T8uCTWtrkVmq(M_U4O_gimE@xxFfH+66i zUMq)r&^JeppZMyE{I5=+ZE|yl((k6%z^pjREwY_VVTObk^QiMnYM!d?mrQ$6DEG^B z+3bFWx9H)HGsu2Y9!rAm{WR-c_{&Ud)=V*u9xg64yTV93N;)9*-rQ>{*LJaZ;pIPs zsFAvLS6Gop;>gxXWqmizV1<~i^Cm&pX6EnPRDn58}6cVKM36^ z@A%~|?YBGiS$E3yvY$@o$XQli)6V=E;)j0^@!Ft-$2wt%R?yE7Kb`y?qI0yIrSH^W z{;fej=V#4kyFOLR7lytZ4JnP57x#=;vkBe3*))FlZmV}T;#r#e((itZKfK#|$VT*( z=8kleWy1scp41lQCm06>=F5vzdH*4C@*Tp1M@p;AO&+?0*}}Wl0($CV$IG1x=a%Eo z@A!JT*QH20{HT3-_knQ-g{i1;-OP_Mu%QvKY8f2)1g<@se!x3d(|lKlbX&SZU6vAe zf_M3WOYX{noA~h+aE}q74DU1o7HY&#`F&7_>ieDX%;|wbd8xlw)C-=`t2GsdEKi8Fg#VP<1B!0oht(V{5IO zU0b%%r;=pd(TkyNI)fjjyX^2t9N%1i$DunLM_h}JsDRr{oMxz)yQzU94dZvXbZ%IUd`s8O@K64YUGzPTc zG2}bh<&tmk#U3%)4zYv7(zgro;XTWkl6%lI`w;VT*M!;?Q@wj&!*$3P4)=@TY-wgG zE+>Yp_Ee#~;{93n9U3}I+@Gn_9+ArL*?I&$d`sEf@3e#MHi_M~@3xs_sozS!vKLyq z+iqx^`6=~l>6h18A1_xi3K6iVgs-93v#~D2?f1xnJpu=81exo*n;mentm26RGyke3 zh3sQrJ~Q$$hNy~#KE^m?!E22G10nq#{mxCb`M^AjBYn2JRLJ77q338*;u$X7vqRNw zvG9>Kz-nzJohiBBrdRQYwlxOww<(6}W*eouxDHi&XS#%qJQAjAwVLm4G;j!Av!hF^ ze3kNxNASg7#NuE5q5kGicNDSs|GJ}9L7a$gtDUntEd0=lQMV=1+sj{kjG>J7Y~6*fxV#rW zX9VaktF|deB#9&vp3L_UNpn%W5_&F6=?;FpJxd9PiO|kcg1Zd>eKq?bE%o7nqPL~# zu$$=Osr;QCKinay{jzJeM&VooVD3_+Tw)W;?C^&+MD9W(_e{R{h#1ahAFeb6)+_%X z=j9Z5CA7HMhs+r$2_~2B{wVDkWLb(mR77#Ec^R(ybv)om5v6+X%Pi3qT5jq0_Qu=p zwtQ?eKT9hxJ$YTcRnykH%V!ypcENewp6@!t&fncpyGC&&$Rrcdxh&E~#&@LSeZL@; zA*}h88h=a2A0?4u^P}JBp<*CBQtu{S69zQ5&mvaks}z{v9`x#|UeX?P#o6!Kg_iqUR_~ zZEkn|FzTQ<8q}8OgOk9UIWHa!GQ8z8wDCXO{y6y8EK^Z!a62OMiXZ)pedEKKC2H#tp)+{QsW{L&b1XU{dU_S%25y z>4QC$cSNSNrh-mHZk`UOWYtBoq@N_J9T})_Sm6|yNY=UOb|CmdH{4SXEQftsfDQ1e zgTO-Afd*KS=9>ehQVWI$D>&33Y6T}=23=szI|GLcfUN+Tbod+K0WMN?N#?efhX?US zKLw$hM+dN{90U-YYzsJm4qfo@W?&2KY7R&~e*|7aZu=VmuEA4cpu7VAJcUv7j9w~B z=EhC>Rp8e*4=2E>#{dJ^MjTuS=YIn>!9LpnN%*Z5V4~3P%_9k>uzcj3OR|7L=F@cu z*GS_Yz!w1a6#$E*IHTLpC_)jBa0zyu4;sLdb3qv~$3C}p@TzUVR(PWgFqfp256FWc znF5C$1NM+U3ve<3iBo?M$iU`~$lW0X#~-;&Ndv@T z^JlvGbY#2GqiI%+WnU{S1#gF$Y=uB zO+T`oIbfLN(pnHs;bGvs9l#}$&uQQW07o_emMb%1a#u%?fZJJ2w@9z5^>7f}8TovS zqfDm8ouUMxB3PtK?^F+KQkNlhX z7o`IB8{4n58iYEd|fV@50!E3aXz#r`hLN)cc z6=o9Z;D%ek7BD2ARCy2~Z$V9cIJ1$G#2MXUT2?nb83*s&3~VJ;3IH)+P6l!9c+&wm zZ!54=40#Dl%1Z{e0meus8OSPuzaO>e?0QGMW$CgZ+z7La4M-e@ghaX!vuoMS?Jq4!`OO+g!^Ug5y5{Iw|ORnBNw4XmCIYz7B8 zu%QL8NrC@FFiqyy2v}$VlYrcB#OtS6lYDOiwuo+C0apPOyqpPWfNSgE?qWb+-lzgC z$Wwg5<4nZcqlZ7Y1~0h=yONN8`oap3B+;HB1qIMylV=DC;XgPb!g)D>JiNsb*d**# zNKZ=c5H_{lpxGgk+BV@%p>>GZ+KvE=<}|9ckJHm&-;GGhsT}y+7obeZTUelE46`#j zAUI)!27OL8FXYkDFl7|50KC6v2M!wqWZ@Vo(268?Ti7b`2a$WQ^1f>ZX#53XBk7cJM!Pl3~HAk4#L#5czk}R`(@EqiSB*OC^6&X&TY5sWWSOtKN&TAh`BRt;t{lx zf2-lW{eXLaTtY;{Cdm=na4PnjK2={f%r`z>X8J zG6lRw0-sLduOk5rA;&u%fGaC9VJrpd0yDq5Kt$pQ?N~cPEuuF}Q%nuYOMtD~PCNlu zh5-h$)aN3j6gXPg)(cg(bTkX`Q#4fra3r=33m}QNy$lU^8rb5ov*T&4A;q zmG|Fbo3#RlNv6%hf!$n7`45hO_FpEmI*14wEtRE>t!alSCqmh9k}t5Cbfg+-b92}m z8Erhx~uq>JK#JzY_^R4FfhS{CmCGHJXQ#(B8lyaFGm}EExXdT!MP0Be)NjwgBex z%mF6BS-N{VJY;(&Jmo&zspO79N=>yE2?coh8nrH(K!dm!)1NTt&Cw>F{^RX{T z4sm?*@;gMAhDq3W2QvLA?*K`Va2T%q1f+@52}?;I_5iJbPRhR}qY)7r(yyzKxuHs{ z?vo`6QjTOCSPSz)00Vi(@98W&y13v+Ka)H>kgU!@&mw_`ng9p*3mK4}dyJAq52h8n z!>GT2jl%Pl6v@mFkOSus|F9H}mH%w{o%UIudpz|cIm`hK;XwQ8IgU2gq;&(%fQ4PY zY#nd*D02vUcpzYH=`8{h750&m(}kdWLsc@&|t{4Q;3yY~q^rnxkst8wy3 zM%>!)yEMk_-hqP?24Q6GQT<3N9Zq}}dZ@-8IX(r`0-4^;gAeUA(L3T!^ghfCKGpm1 zGOP{1*#@YYn~@1iTXE0aOF1SZ36;jf$Shh zs99BXk}9h1m7N3Wh!SWHzzAXQ31`vdk|KHptD}P6m0n-#5k#*$OD{}fb|xO9mV2c~ zlTc>#(!P&L&e=$|z5!)o=H9#~U?s5-jQ_eKQ=^hIj2QQ?cv1 zNQZGx#mBfjKjghM*TnANQ{qT4-X>&Ba2*~!Df<_KR&G6+ZpZ_QG4 zc`V!zD=Q>a2ww-@CrxfsEOJ{L(Qj0J?g9p>AJ%mzja<*&O3Lng=WK>}6&66{e>m~y-+}E0?(5io) zSG(#a&t?RzEl;6UZAZDY*ELi1t84T6oo6CgQ-9NbDIswK+Id;XDf~Klw0X6iJbFjj zIyhSmGy+_OfLFuLi@{Aoab*C6S}O9kL)+8kL}{pdVJK+?DLWMmJ#e{s8LvBpqME7J zxaWVS`V+rZUl<;gjBwIP-jN6!TrXQ^eixDVtma+r{uDoYFR%!XTmmkbcb7g}BV_?+ z!&bwoY6!0jwgBQJB+&?ktCxe8>%Ww9r+X_n1OA+fP-p|XptF2!dwKiWt|0H&BxcEV zwXTS2D z@uuNXN=CgMyv^mVT0J6kc#|@yDb%^FY(2B~WZ9&tiz((G>vGMA%cxUF5IM|?66D=9 z=vG)7bfPq<9cdZ$c?=h$rmA+?WM-%D(74&uapd<7%oAxW*dhv0A$YHFiPiEBBZo1Y z`v|}|Pr#k35%3$MvD_N6fK@{w^mPj&IN8{60a`Fk2qMf3GqGHF8o4K!(`iJQ>11MA zwRAG84HMM5GML{`Y*U&V=+PcVHMwNpHt?b5-sE$&@{FvJ4|*eUYV)}t%-8|WYT+JB zNpaKk7PKIN;KG%~*Dr5RJn)5LcGoDyQ7xFeoULc>!wU@y?Y8BJS~Fhw*NVmwo`*k_ z$RUV$$vT}Pbe~iU_6q49Z@fU9d@+TN9XY2rBCAp_+NdvU#&==gQY+$LDK+H$ z5q=guA+fZOyomcWi82${+xecLIK_6_qYWIqw+ zhIupCrV3T@@X?bWTtdG#M0u<&+7t>&aVABkfu$Oe-C@@f>dpT` zcW$nF(dp>9QNyo^aJE|wf_`xA!Fac0X!$AAE&oOU|3KikT)py@K~O)92dR zm$=;`ohKvLAHJwNqhLu|8kc(P>D=Hgiu;}z9%JvDnFCx48z6jGrl78HCYeIqv(&3q zQs%zfm87lLIEe-eD(4zXK4yg5o~#y8ta|f)V{hxC3qHXUn%-@iUhfypPrEqr)yjQ} zP~lD>2yQ>uc^K}JS==Qxb!?Ak>Y|UQy@oxT<@Rr~yLL@%UciS8y(67>=eM0Wv%##d zHTUqN!jtv4cIIj?uzF6&%G;@S#mCdI>YaPH#e3=E`yDU0z6j0v+OL{QZ+9Q>5!12P5AMpsJbRQBbb7(HmWU_TlA|jEhZV!% zsML(ydZ~tTYVfef!~;o3w{U~E2g6IkEMNP1vc$7(({ndH8rOJmT|8ORdB@kCleGR9 z&D$?)KU|eBaY*s=fYzMbq$t9dV?HU9HrqU?H*qW3ZM!WSZ`cL$V$GCDH8+*~Hko*c zCdO^F8H+8*Xx?mQn{6tcd4Fx&)yNVid>WEeZrTfHBBQgVlUTrQSti)Fc7*yMMSuzFkaygPFh z&x9D*;g2{0r;A1cfJDt>dVckrzZ6m4+m=sysJLRU1V-%G?noK7&t84G=gFoyd7p?q z&bLK<737UiY&f{m`ksOQ}M7V#9ZUgJU#sK;MLcL9Y+~9Pmnw-@V=Ho2AXjR`$!)6*FtM0@7LiGPZmclgyuP$424!GkawlsBnp8y;{7tn+nd6r zQ4MhpL>Al}qlQ~4aAu>BW9)qJboO3Lmu}QunB6Mko75vveXFsAaZlsgiLLcC`_}qW2|X&SuCKJoad>E~KD(09;%@{~ z4L}|B&JO3NiOosVogIhQH|p58bcBtO61`pKCliyiednq!gK&#AmKo({TM?Fyi`H8b zY?DTx5no%p42!VveLTo353gXgR(G>teKfdk;cJWL2;}+Gg8;!IoN;flS=!q0RT8RCV>$j-P7+zK9vh<0hS3>YK|R<$X2dhi~1{=&nY& zl*>=@9v>R$4r4uhMCFt#49NBD6aFlX#<0&6~2bXxUj6ZN*@-u!+ECQ zGUfG2t<`!hlV*D>J_f{PK5TS%c}i1OIsy#O3DH&U>}{%1i2x zl{0U5){l%^Q^v{DcK3$LO`cBSH}H7vVpK{!D_jY>uC1|Uc>na2i}yHrpqD_85x{G4 z*-WdcfB(kC5EVj0~cT(XGQJf9qPdIpJwLGOR) zWLEP6`O3u*TxML{78>hYOI^(5tMt_-5N~l@1`r&6^yJzW(xgL5kT-=b`yhP8%YjIc zT11p;z)|hs;AK&>Io^om*nr5~+JD<&s~Y$68zhe5343Q`y}=aeN^;jOinoPqStw^* z(bs~PPppWE2xd>R-_=Gj4sdp7V*MpyS3^(}L!C;_?U<0P73H7((wsnSVQmUbwtv(R z*BIC}f_>T*{)KOCcTO$?TA3J&OS?!8jf`(!WX%);(cBZRmZ`1E-H-0#Ghw7AgUe$Y=*F`t3bgv5=M zKSqdau`oN3?I04DvUlpbeC;X=RkLtOVg<#*tenD^^^(6w&_To44WRz=OdTus)?ogr zR!rF8PsuC_k>@WN-RbFKE_W}k?B-M$#r=D5vaM6e`#aeVY{DfKb)bJPbbXbEPM}_? zZ5NLI_I=;eX@^T{TuBzZ7Wty!HO*l&C)<@$taguZje(JcH7r5R6^SKXjJwP*7e;sE zTk^e-(?h&=UQKooJEAqEf?I%|+kmUK3k+`fo)KLnBSSgNdjsYv1!gmAN73%d+5fAS z`wwd3$^!sS-kPKdKP3E6tWpGtkdTl76@h^KP%0GVmsVuJqGy3>K}F?<2!!Z)tq2;- zB~Of0O06g;R7G+U0_NZjX{vx|1EmB*`7t06BVq(D$!%zLE=;jWSOt*XY=UTkZ| zQ)L_O5tMEM4t7;y)F#l>`1ZUQk)rh0%%@6zdwQcr_j^;}0WUyX?RPmO=VgTieQ2>; z2Paq!cPDfg3Dx)!dDDldMzCI%kc~d!hzG~A+l>`-S9L33egL2lUG6D7^$NXbYcfn@ z0UDMycYJ&tCinr0X-r$+bjIt}xK#bt@^-@^8eTBPO1S2~5pw(ig+P|a)vI!%bS6;g z4d}SWDeN1MK^PeToCz_9CFGe$%X=59^|KOoE~>U)mE)nEuSUZjmo&d78Y2AXO=pt~ ze$L+%uV=HW`hG75fARrbT+Lv1?{hpXV*w_bKJR1!H`wF{+|U`e+7B?1`-=*@8dGPo zBq-9gLp*`bbG8XWl#fV}gF8`80>_dP5d^Mss%bOP!ud)wD>A5#_)2hpC zty0C#nPS-8x&x(tfND&tnY_`@0K&1&5X}{$eb`bZ9Mhb+pEuhHBLV=Gs2VxmB`Bulp;-2jtU^AqBkVE9ihK6Te-)fEQEI3zS<%p37P`O;LFQ(F#|wQ9 zD)#(_Q*9*PEL*Q)o?H;8{P8l}V`1jb5CiKJXu3?6*gN`=Rao%{Ty@+2%+T%F>0#SN zB9(CuDpK_fWKWs#TnyUOCSsxF*q7@B2H6Uxg@X+y@pWBjJ0imt+93!7c zBn-GYY7B3e8vn+7Gow4~QgxO3jhvOa71{|se*=ru;=DEY+l_zorRf9;)1!}Z3B*aNvUb2%dWa8L@cy1%1B7Hp0L z>u?$a1BEdFUP16o6e!{`-EdI51q{Hm5n#Y#V(&fiA!%+-U!D|edmsQa8iEp1_Lk|J zm&x{3i=UMaJ%UZ`>ldslSQG)|lo|b{Zto0RFF|l8g@0iMsl=6BAbG!%JmSywiS&Lm z7mLlT!K-LCF?W=#a15_{yP2qD8*n01UN+5FtiK!F=_*P7p+AEfZGF99-FGNV#CSlM z^erxal*YZu*&QZIO9wU>-o-r@BFk(LN2coAQLIDG*&AzB8lCjr1(Gtl{QdL$yF;%G zzk5zdQh{GZ$ICz&nIg~+wJ2_lMOHf(WQ=MH$E_SYmSfHvUDuncsovLF5rZF(c=2yW zNd>$;f!9Ri{{($xs*LS@EAxggmsV81t#Yh#OL4UkM*#@x3;KZBKmUB1olw=>_xvxJ zJ8UlV)Th!!pG^RZwRD3e;v&by$5?3{zYuV}?r+Lp6JW9A!4^*Cr8s?pVk=quQ_(M@ zm6IuQ)3ZGFk1UbG5>c#JQnYh+NOYf$IAPeKUo#`l@{AiHb|v1)F!LBsuuldkc$)`!X!^JxpvvA{9zn+Y6EWB!9A>6Y72B zVtCS1c`nMgFI5PG^H0hoB0>4vOc`G?gT@t% 0 and periodic_boundary: positions[d] = ReaddyUtil.get_non_periodic_boundary_position( positions[d - 1], positions[d], box_size @@ -1287,7 +1283,7 @@ 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 +1301,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..1a4b0ba 100644 --- a/simularium_readdy_models/visualization/actin_visualization.py +++ b/simularium_readdy_models/visualization/actin_visualization.py @@ -9,18 +9,18 @@ DisplayData, MetaData, ScatterPlotData, - TrajectoryConverter, + HistogramPlotData, TrajectoryData, UnitData, + BinaryWriter, ) -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 simulariumio.plot_readers import ScatterPlotReader, HistogramPlotReader +from subcell_analysis import SpatialAnnotator, CompressionAnalyzer +from subcell_analysis.readdy import ReaddyPostProcessor, ReaddyLoader, FrameData -from ..actin import ACTIN_REACTIONS, ActinAnalyzer -from ..common import ReaddyUtil +from ..actin import ActinAnalyzer, ActinStructure class ActinVisualization: @@ -262,392 +262,19 @@ 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: 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. """ (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) @@ -655,7 +282,7 @@ def get_bond_stretch_plot(monomer_data, box_size, periodic_boundary, stride, tim # stddev_long = np.std(stretch_long, axis=1) return ScatterPlotData( title="Bond stretch", - xaxis_title="T (μs)", + xaxis_title="T (ms)", yaxis_title="Bond stretch (nm)", xtrace=times[::stride], ytraces={ @@ -669,8 +296,12 @@ 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: 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. @@ -680,14 +311,14 @@ 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) mean_long_long = np.mean(stretch_long_long, axis=1) return ScatterPlotData( title="Angle stretch", - xaxis_title="T (μs)", + xaxis_title="T (ms)", yaxis_title="Angle stretch (degrees)", xtrace=times[::stride], ytraces={ @@ -702,8 +333,12 @@ def get_angle_stretch_plot( @staticmethod def get_dihedral_stretch_plot( - monomer_data, 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. @@ -712,13 +347,13 @@ 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) return ScatterPlotData( title="Dihedral stretch", - xaxis_title="T (μs)", + xaxis_title="T (ms)", yaxis_title="Angle stretch (degrees)", xtrace=times[::stride], ytraces={ @@ -731,15 +366,12 @@ def get_dihedral_stretch_plot( ) @staticmethod - def generate_filament_structure_plots( - monomer_data, - times, - box_size, - normals, - axis_positions, - periodic_boundary=True, - plots=None, - ): + def generate_actin_structure_plots( + trajectory: List[FrameData], + 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. @@ -750,118 +382,148 @@ def generate_filament_structure_plots( "scatter": [], "histogram": [], } - axis_twist, _, _ = ActinAnalyzer.analyze_twist_axis( - normals, axis_positions, STRIDE - ) + times = np.array([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]], + 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 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) + 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.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"] += [ - 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="T (ms)", + yaxis_title="distance (nm)", + xtrace=times, + ytraces={ + "<<<": np.zeros(times.shape), + ">>>": 85.0 * np.ones(times.shape), + "filament": np.array(perp_dist), + }, + render_mode="lines", ), - ActinVisualization.get_twist_per_monomer_plot( - twist_angles, filament_positions1 + ScatterPlotData( + title="Bending Energy", + xaxis_title="T (ms)", + yaxis_title="energy", + xtrace=times, + ytraces={ + "<<<": np.zeros(times.shape), + ">>>": 10.0 * np.ones(times.shape), + "filament": np.array(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="T (ms)", + yaxis_title="3rd component variance from PCA", + xtrace=times, + ytraces={ + "<<<": np.zeros(times.shape), + ">>>": 0.03 * np.ones(times.shape), + "filament": np.array(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="T (ms)", + yaxis_title="normalized peak distance", + xtrace=times, + ytraces={ + "<<<": np.zeros(times.shape), + ">>>": 0.5 * np.ones(times.shape), + "filament": np.array(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="T (ms)", + yaxis_title="filament contour length (nm)", + xtrace=times, ytraces={ - "Lateral": sum_lat, - "Longitudinal": sum_long, + "<<<": 480 * np.ones(times.shape), + ">>>": 505 * np.ones(times.shape), + "filament": np.array(contour_length), }, render_mode="lines", - ) - ) - formatted_plots = [] + ), + ] + return plots + + @staticmethod + def add_plots_to_trajectory( + trajectory: TrajectoryData, + plots: Dict[str, List[ScatterPlotData or HistogramPlotData]], + ) -> TrajectoryData: reader = ScatterPlotReader() - for plot in plots: - formatted_plots.append(reader.read(plot)) - return formatted_plots + 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( @@ -936,14 +598,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, @@ -953,43 +615,131 @@ 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"), + 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 save_actin( - trajectory_datas: List[TrajectoryData], - output_path: str, - plots: List[Dict[str, Any]] = None, - ): - """ - save a simularium file with actin trajector(ies). + def analyze_and_visualize_trajectory( + output_name: str, + total_steps: float, + parameters: Dict[str, Any], + save_pickle: bool = False, + ) -> TrajectoryData: """ - 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, - ) - ] + 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, + 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) + post_processor = ReaddyPostProcessor( + trajectory=ReaddyLoader( + h5_file_path=output_name + ".h5", + min_time_ix=0, + max_time_ix=-1, + time_inc=1, + timestep=time_step, + 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_compression: + print("plot actin compression metrics") + plots = ActinVisualization.generate_actin_compression_plots( + axis_positions, + time_step, + plots=plots, ) - 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) + if plot_actin_structure: + print("plot actin structure metrics") + plots = ActinVisualization.generate_actin_structure_plots( + post_processor.trajectory, + parameters["box_size"], + periodic_boundary=True, + 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