From 2079ae4c7b992bff444749d3a99800df2a019686 Mon Sep 17 00:00:00 2001 From: Deaglan Bartlett Date: Thu, 24 Apr 2025 08:12:30 +0200 Subject: [PATCH] Add tests for TFR integral --- tests/tfr_integral.py | 82 ++++++++++++++++++++++++++++++++++ tests/tfr_integral_params.png | Bin 0 -> 21128 bytes 2 files changed, 82 insertions(+) create mode 100644 tests/tfr_integral.py create mode 100644 tests/tfr_integral_params.png diff --git a/tests/tfr_integral.py b/tests/tfr_integral.py new file mode 100644 index 0000000..e7a78f7 --- /dev/null +++ b/tests/tfr_integral.py @@ -0,0 +1,82 @@ +import numpy as np +import matplotlib.pyplot as plt + +from tfr_inference import estimate_data_parameters, get_fields, create_mock, myprint, generateMBData +import borg_velocity.utils as utils + +def main(): + + myprint('Beginning') + + # Get some parameters from the data + sigma_m, sigma_eta, hyper_eta_mu, hyper_eta_sigma, hyper_m_sigma = estimate_data_parameters() + print('sigma_m', sigma_m) + print('sigma_eta', sigma_eta) + print('hyper_eta_mu', hyper_eta_mu) + print('hyper_eta_sigma', hyper_eta_sigma) + + # Other parameters to use + L = 500.0 + N = 64 + xmin = -L/2 + R_lim = L / 2 + Rmax = 100 + Nt = 5_000 + alpha = 1.4 + mthresh = 11.25 + a_TFR = -23 + b_TFR = -8.2 + sigma_TFR = 0.3 + sigma_v = 150 + Nint_points = 201 + Nsig = 10 + frac_sigma_r = 0.07 # WANT A BETTER WAY OF DOING THIS - ESTIMATE THROUGH SIGMAS FROM TFR + interp_order = 1 + bias_epsilon = 1.e-7 + + # Make mock + np.random.seed(123) + cpar, dens, vel = get_fields(L, N, xmin) + RA, Dec, czCMB, m_true, eta_true, m_obs, eta_obs, xtrue, vbulk = create_mock( + Nt, L, xmin, cpar, dens, vel, Rmax, alpha, mthresh, + a_TFR, b_TFR, sigma_TFR, sigma_m, sigma_eta, + hyper_eta_mu, hyper_eta_sigma, sigma_v, + interp_order=interp_order, bias_epsilon=bias_epsilon) + MB_pos = generateMBData(RA, Dec, czCMB, L, N, R_lim, Nsig, Nint_points, sigma_v, frac_sigma_r) + + r = np.sqrt(np.sum(MB_pos**2, axis=0)) + zcosmo = utils.z_cos(r, cpar.omega_m) + mu_true = 5 * np.log10((1 + zcosmo) * r / cpar.h) + 25 + + gamma_inv2 = 1 + sigma_m ** 2 * (hyper_eta_sigma**2 + sigma_eta**2) / (b_TFR**2 * hyper_eta_sigma**2 * sigma_eta**2 + sigma_TFR**2 * (hyper_eta_sigma**2 + sigma_eta**2)) + print(gamma_inv2) + + delta0 = (mthresh - m_obs) / sigma_m + den = (b_TFR**2 * hyper_eta_sigma**2 * sigma_eta**2 + (sigma_TFR**2 + sigma_m**2) * (hyper_eta_sigma**2 + sigma_eta**2)) + delta1 = ((m_obs - a_TFR - b_TFR * eta_obs)[:,None] - mu_true) * sigma_m * hyper_eta_sigma**2 / den + delta2 = ((m_obs - a_TFR - b_TFR * hyper_eta_sigma)[:,None] - mu_true) * sigma_m * hyper_eta_sigma**2 / den + + delta = delta0[:,None] + delta1 + delta2 + + fig, axs = plt.subplots(1, 3, figsize=(15,4)) + hist_kwargs = dict(density=True, histtype='step', bins=30) + x = delta.flatten() + print('delta', x.min(), x.max()) + axs[0].hist(x, **hist_kwargs) + x = gamma_inv2 + print('gamma_inv2', x.min(), x.max()) + axs[1].hist(x, **hist_kwargs) + x = (delta * np.sqrt(gamma_inv2)).flatten() + print('delta/gamma', x.min(), x.max()) + axs[2].hist(x, **hist_kwargs) + axs[0].set_title(r'$\delta$') + axs[1].set_title(r'$\gamma^{-2}$') + axs[2].set_title(r'$\delta / \gamma$') + fig.tight_layout() + fig.savefig('tfr_integral_params.png') + + + return + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/tests/tfr_integral_params.png b/tests/tfr_integral_params.png new file mode 100644 index 0000000000000000000000000000000000000000..21080737b24c04125fbf56c97b603c726de40929 GIT binary patch literal 21128 zcmdtK2UJyA)-`%9ODO};Rc0huB?csz$kBv~fCz|W5K*G!Bp^Yp!jh}as9$&J`^rW%{g+i}I{#zC$7O796 z$WBNcJAC#^;9#@8lES5h#l=l#?{DvVyL{J*cjlG*EFSIOS2*l`L8+Tn`rMeTVdWOl zvo7zS9hEll3=^Q!b`3j9^XhC~Yqh~nxyRnZXKX{H=6z1@HuxJ0`DY)fC>)NJl)V;8N>-{7jKuTpc5+~r??DA;hul;+D11>df3 zTK47rzGW+Hz8)iW=o!=3_v`*k4nLH;hn;=m^>(i)-SUTh?~_7}vj=$xn=-3fw8$;c z-MlV@|2f;KXP6#-{P=N0lJ-r<&mU;!r<@kkhcbJZgGB6WTOLyneObD1ZDGq{a)+i| zb?P=jq6?q3X&Lz|S{{?n=x!WrTIM%YB_mejw_khx-~P7E?PZ8@bC&*#U%%h8d-oF? zJqo4Z+r=@l#ZZf$vS^uL(Mh@$6pC+kMb1q7OSit3-0X!>-~9J%aul&Htzil!I(IIGfKPWE~e2vf`6Go|#`?hnLtn1%-Pgo z(_7CSJNEODBS%z(q$r97oR=tm25UFKZId&l0&9nyaL4i60tosyAB z-?@lkQ==`7r7f%|^4Pjz)23a~A1uV2?~6ZUa*o!p?b?80=}Rsbt%^~IwfyjWMYF?H zv)#gU&ma>6g<@2ES4K2C&tbMnX0_dGD2F8V^)zxVdQRh0l{G^;V^d-YIH2gyYF7H=5+n zTMxG6sz%BL-}mzJV%q*K<&lwXhX>y{X57ct#f5Z1n<0)~^8*W0`HSi|>6pBOgEeoY zSq!}2Khaks+0xQtZex?4=Je?q4TGSip4$sL$|2QizZ|>1>gbtIh4d3$zM3kt4|t}Q z>j{>b(wui~Z@*9)B&6siGU-?9^^F9B@RcW+-wNl>=|p~{Q1k+3riS-Q2wwU0h*f`T zs97b;GI}G2Qbp9yLaeOET(mwmWm@dnyEmO{ed999J{j$?KH9CkO6A6hTD48|k9v1! zEZ_GsLOQ^5Va~SMwkIgHg`I=r?Cj*A;rw|_a85D@%q=Xc$nCpr+S`YF$_)t#;nuHY zwZDnL=oqr!cA_WPDdGC{>w;^B-p)Han2m{3BKYmbqgdaxwzfv(#>AW=&-;CnwpA$) zW5;sJBc4=~Rt*um@rrB9X~ygF7Zhmd7%T9=uC%dXO-11wXD0e<(oHmf`Q?}T@k>KA z1729N@}FmBPGw!KS1~nxk~*K+^RVi{L0#@k#p{K{S5gi&*5*4oIgGzMA&4_SwC&!~ zAKr~wl!SGmk?)v0l)tzjm|oaK|9%s_^osMY9k0Clb;V{MWLfn`baq~}8h9_;UmIWL z&tvFw|9%9X;?yMF_Y`MdHa0d1r|A++%R*Ydk0pEaebUp@bB`ikTIuO&1o`?#9=@@X z4Et(&`ovVGZpVVWjQr&k`SF2z`2b!cx%1~AIZkD__CCJJTshV1l*p-?cIGk-rNDD$ za*$oZuIrxS$&6(b)AN|?b%~mZw{PDTbQ}D^LhQUrliW&Y{h^Ewx~}f-2BEbS%BStG zJejIe>oj6N6y9tUT91gdf%XmgGtGXjXt5%4Yt*=0@C%uqVvEz|T|-00!)FlcR?#sF zD-F+5pN`}fJXajkoUNGNVi%iRh#%vB&SZJAm)W1>KJs#2#b5Zp3hQkMUIT-b;gT~5 zT-A?{(jZ{f3o%kCeU~=7lm!Ze78D#(w~%fWHZE(TrI_w~+x)hfLg}e&%FHg^)RaFv zBs4KlkN8ZXRP27`#iHoNEHc34qQ1{&IVJFW3qA`zNhV`f=;R81o3!TjHKr#nj7N$^JIqebj91Cvp~%U}aU0fhcrXc_ z)6&wi8!P3rnjBD$S57W!c}%_U^Q*~C7g%oEt5J%XT^q?4T%9mS6y z@0*(*BUy{rxG8F3epc13&(YD5cYc0u&T`;=l3lv8vT}5B@u~iPJl3)0G}~^qhEvow zwVF_zYmYo@E<3$odb)5iFw~-E)7I37UXj^PBgGF+FQ+gTmpe`#?CR?3MP8{u_~0>X zKi1XRsnlv?VUgTxviW#?TJWjUr)vWREL2!ED2vaATLmI*T< zaz)J{{n_az6@n*ROYZsIiLBb7AT>v!+-&da3hz(M*YIBd{r8X6v#om*o#ta7I!r1{ z$jQpy`R1E%;^y7l+zf*z7UU1=ygef)=RIMito*=ka-ja3cgu~o=1VTi*pt_arFQFj z_Kl1@d;4EjSUt7bOriJ{?^}ytzK5Ip?83rAadC0wqr*2?B@`4Cir>6>cYVd0zGjz~ z7w){gpzmS2?APB%qu{7KdCEFvfzd{dX}urcg*N0nQarfJ-B4ck4A}t2~@WT(5GvoTS2uG@Emt~P12N(q?6wkd# z#%->1p;Xj9uyzqT(Xfo7*H~Fq zWi|d$7K0v>k)c#yUtfD6+j9C^bh2(a*23emV6nu&Kt_^#_Fte-0$%~CO?(XXGN=rb z=&MT%#>%Nq)|KJZ%*!cT(AU>rO-Fa`)G1fYxk!n-J9};NXM;9!s$3jAc=`RlX8Yk> z@;5W%y`j5;ugdYYUQJFjX>I9Tx+FQ%M1gh@T>4;Ry3Em=bPAU*$E{wyx^FF**O{&6`=4y%&tKyZsF<=O^kUy^S-wHZloXo&@Y;;JvJLcamKKac-kqxUBZK zlmZ36$%?bO0+nUqM&xTLx+E@q`eb_P(xv(o-Ak9)_Uwr%v-)ZKcEgF*r0tV>8WJ}4ld{QdXe z_fo4Urro#loyOBclD1G#Qi^)=WZxwnojpL6-e>XS_3PglE>={W z!Dj|o?<#4RAAGZXdFsVCKfHYTQhE05*|WQU{PFOv(}A~vl8tjlj>l`{svr&-esV;# zh}YTFw8JxnnrSL35f*%I-#*i0Z$ges`}k20@MOxh5xDVTcWFRf9F_}3HDY~}K;_r` ztk*BamVp`Ym7t~=)xBSuvgG2`{+AT>+RfbMhpL#^xJ=G2mJQ~t@9eKS0&sy>o62h- z1rPCO&8w`8dv>tM*zC~h}y{i=hH2&W422}4fV zopy5^jbylCof!4oCbkW)q;PLoVY6#PXMf?h)SEf9=12133gA-eF&O&)50P!^CI6?) z-??`1BYAK@zqd=leDLk#Ey_$B?p0N=r+1G6&Nd2iEVWp07KE057RKckhP2 zVz!qGcASZR`t&Jbf-hgZhzbdL9?=?P(-Qss`H_C_2M_ckf29V2eZF|g`1=H{st9S` zL_|L;tKk-P`^kZKOg2qgM(Ii%qf=8RgP*A9&@&7CfRm^k)7Ot??ak^>NY-!7vWm^J z8BPFIF&OQ9vxkc-?Zg%8!G5099N4?IkfeTNP(-8;n3?YfP+H{ zt9$jTRqtkp^AnS_U)5&DgS?bJI2bA3v4p_Zf;RwD^xfL^umbBm(o1YX2}~<0-Q>_A zXYh8*&!0Z9Pd?bM?Q%*|GQMP_<5ejg2tnj+fBzUPW?OrEcBKBsG^5ybr_-lBfeE?> z2jg3E?b$gwOS|Xqu`3!NFZDHLCL^84fr<3vrqX*4+D~34ke@@?Z3V@doeHC+v^Kji zJI$FHy?=6;#bNx;k)uZwv&y37BQ57Xo2~fvTeT#wmltG_%Z)ZdDDZc8U%O-HPJ`Mw zWs>z==$K;OzgOC@etmVln3$Mm=y8wAipNJm97Rh51?1%;q~0On#yg}LH^%@}pJ8#D zO9D^bz{nWq=eIdX*tSwfm9N?I?&>vb;!6n5w2VW3JZQ{@2_QA{UfWB*)St(0s@bLn zWM#vK4HDjhR(-52EG+kYeGLKHT;1H{F=C#~A~9f@F^H~NQ!Vy`j3NLlMkDPn?Z7m7 zK!PRYkotK|TH?A&eAxvAGS4lhq;RjGp<&;#L)^PMQB#cMT(^y!a%yUUtiAE+El-to zyNd7G%?_qde7MPC@cidhwbu0X^koJ#lts(RtXls`U;AQ{nTkKcJ=cU zgzgOt3^6p!R}Qs0&Key&zMSH0c(7~jl)q%7oV+~G#n*q!1%h<-@Thrqg9gK_7VA_R zZgh8-Ol@CHY^wGvx*e>nQL5?2wc7`gbj8KR-!wL=uGz$OR!Yj9j*3+9GI^@Wp;X6IdS)FG$aSX zSxxA5zyww;Wk$eLz)6*#J`~0q-JqebMP96iUf=-L${#rchZ9a^)ms@}=Aa0i_w;Gy zOjcZ6+%xSJl!E#Pz&Dk=)P<{ZKOi8W)Ua@$L8^MZllq`{is^?V@#@*T!q4D!7y z+zf747d1iNh;aJ9VIdbhpx2aPhA80ET^1CJ)PV@jcJQD)Vz+k~encozn!V}bvIn2Z zN1m>EdZ}dn=FL$%B=0K)2kevpj{@YXgb=BloVZEj%3-jH`>a9yW=H5nKR?2pu8v1Q zvTnrGt&)kKed%-?+6L>o)k9)=DT6Kc!>MV2kNQ~JRtf6>{ha1HSa^*Tvr>#&a%>xN zY%M~$jJ859Tz8@)>mNdl_sU0(TnC=&yE^Q_*Ai8;k%=ij*M3UF%>INIb8eq#UA$^t z#is^gjFI!jO3q6eDWwA6PS5PSwgknF5VQlJm6<|)h!_?yCtzMKBkJ$*bIPUv1Ze#) zaoHaaCmNGl=@)U~#w-iK+B4|uc$snneo;;wf7X0#mlJ>F2dA=knSE{^D&?=L)Povr zQ_q~G9+t6^+;nJYn{fzey`Z9lEL(^1yNTvTzAL$3K;QK{RHV>0M~1=E^Sb?-KT>sk zO9m_TJK7&dQ@;}b{L6yS(jWX69==c9-2AC=?!-B(j}?a@Hq}AfFmwPL;&oVCSTRta zykYa^-AIij89}LiBL9MVD}HW73EL+Qn80j5u-{>(H6f}!WJ^(bx#m>0_qa;7qSSz60x8B=qWUGTWd3s z-s6;$xvdA3Bq|Yl*#HJ0QXO~tT*xx;{Zb5%_^=lBsEVI}Xm&`DA`L5!(nJb!y>%;T zGD(TU)zh=KEJ$eZhwapp>a>&N&kGjr^YQTsgP^Yt1q#3E^&Rdw9uV-fFGBGik7;099ioiRZsFo>5G=QG6pc?RYS{6~YV5bG! ziXto6rfO>0g_UC!6KkPgf}K4=^@ELz>$G=@RsT7W&uurL1b8DuZsbyn0_0-iH$8Ot z@L}nOz5XsaTPbQc^;b;czDZiy!@*H4&DR{}-_pc~-pe?~^ZQ-lBbLMMhq%YsfGQC0Gy`Hel-L z(W70%!_WXYXc@RqdY@HL2#27A&UYS+&hV0A7* zEjqo`QI*B<^W&XTeidQ!A!djg1rt>DB-?QH=!EQKie6Acl)G;1!EBdn>>!T&{7g!`CNw( zrmlE9k$xRN<pzElli6}t|QT1HXc5w{X0Ml$kDt-^N??&Eg!0zI%3?Z z-M>!Bcm?uxotzzN&Q?%YSML#26s&CVn?L2dYXPcn0GJGSu-0ZPE}oaxy`#(JQN@i% z0uB;$(hv|30HROUtiZCmQzUvNCt<-VRk9LDA=P?Ng)j}MmZc4Q>4j;Ks7kpq(qG(} z%<{7e3%W(k@oiFB)`NR&Y^GE)jUDKjr~W7cRpsfu`>8ZoOcToFJLvXJakwIuD;i&f zr~*zx1FefgJ4HUj&+ zZ_1wskE?atn-?V=z{~d2Pa*ZQ_qfD9<)9|Oy9V`2@#NtTEh;A(k^xl&DvgtrE+#24 za4*F+#>3!~MF4=m{`wu%bAy+z8&Ft2bMuU;7f>|*pQJKx;B z7Z4CTK5hvZTbosi7~D7A6I_SVR&}i6UPH&`JO_L5sJz2VnvjD_b~4kFj8affBp;N# z{Do{1h!Ey5`es-0+qY*=p8Sg23f?x*K# z{6*}oZ%kq!sjowZ`@`G5o(i8-r18u$4}JM1x5(92{l%V|fJkbB4fT)646dv|_etvOv44?RAY=nP?G~p>}JFt)RN> z%NI)mgb0aZN%*o}q#3^g>NMB;5eTa-b6NvDOeU{bG}hu;Izbk~Nu=-;T~AyuD$O@1L66Dc^?Zo(UY z1BD0!BW;KksSop?-!p*f`9Ob#%}z%EVjc)aRq8g);Sg;-m_<}ru70|fH{HE8cX44> zQEeAw;vhk*r(D2UmtR60GhMr7%NCM%0JTyeu^>&K#hoM4L?t^pI_4esUJAE4erKtO=Jxl{81nDnLpvl7!5!io0XL%JNtKJ zlHB=BCPvC5PDATw&Gls9E6j^aN+2Oc=Hv)V_$C^Zl$6MW^y0&bYIzwpNHt(#Wa9AH zjq6zrHYh+;N_1EbN5HjvK0Y5Y4x?jZ$@)`cWAd`HvO)bI>oxU@NN;=C*z`Z;p&%9q z)~r9)U8<649u5SC1H!>Z)ycME) zsT$MAuL?=^IiiqUO#4kCIpfQv?$0QRY7qDN0%yD6{0m7+eYrdw5OCw_bu)w-UzFdm zCM(kU;K;uqaa=a%fJPu_`SDC{x#qIp9~yHG@%K?od4Z5;8e114>NM=7F%Fg0(eess z&OG}3`HFgkS*2jSMDF z(xh!)bfgbtYS_9QChK*xhVy36o;w$TAx~>*!@~P1qDgKi9u52`if~q;a+G@ErK`%i z#W#0%r5e=ikpLzNm-NNsmmo35*f;?vBv_}PH#R25@>AR4{$RtM3lK*7XZjO;Ep6Lb zMC?+o-?*WmrWR*A*1aEdA1}4C|0A2q>*z#45N4G?awT51NqSmZ1$-89IwM5*$Y2#v z7RqI)Y`X1B<&jGe2T=`rTsGB$8m)>Eu?=L~3~xxU1HZ5f&izRObw0zh?wGVz5hC38 zmcSaWyW;-;U$XZfkjI}=r|{Yn$ijn6^Aq6Gp+M~~ytT%P`hYsk!Ln5zfv+`W;i=$;wh4Nd% zsZlh*-zm(ihfKeNBSWf<`o}4M10>L666W4O{Sp8Cr-huQ@A;1nqJG`rv+rPOAj$(D zA73ugH19B~*3dOcq8AkvC4>dK^7-KNRDiN?3-Lak_lN8Nd*&HQU#|GE`M;DO^t(q{rfs&3vjmq6-R31+P#725djz!40egPDc6P8i4A->X_lEM zsN4|M)k)toz{E>Qk$~jAXU`rvDCVgUEO!eEsvsBTwA&$Wp}L){jKYB7fddDKsY5wg zM<3Z34#KDl+S-UyXLJ`4V~D{M8bS@~0^AneQqm6&i1*hercVmkPsZ6z4aNIv*hxSk zz=b5-1WRXlRp{is79uXmIQED!Ky3~Z*s_(4)(!nD|Wc`l=ACtdJPlegJ003&$UxRy>PX z=p2-RDpd~e+R$^2f09M|@POT3(m+yyj9d9o#NME*2nGtpoKqZpe0AM12P7B-%)=0S zL0CCuWXdc-O>W)F51(W=V^M-m zz2i5l(9<;7gxE0k02Zv>p-=TxypMg+*QcqNT|X=R#eA_wF=82|;fE?_|6+ph{|Hcm z8@+*5P{~K_j&1*bS>nfq;dTOoP%o-ewZv_BliBWd9@OEsg0Iwhmv4>DJ-_^Xr=WQe zVc}||r;Uw?mebaz;o;%6U7trfXd7^=36JWRz~~um_#+7hSxL*lDnTeJ3{K;(#1>?! zGQ)11ch#y@uD||TL2@2QGf_kyYT=9(IXO8{$&4-T`LRN}QXkoq^!oK{mHG6`^&!9h z`s>Q(qSREct8Yn&#!B%$>rmrrD;Bbf(r{MXd7RkQ*|h|63Gf{s~3YJ4~gZW zT_2laFxYPYJ$*a{LEF(Vn2#*sF+>x~~)4w$Zv=}e7Xa8{`&)p;K4TDmcxZ6fH zgc#&)qY7fl;FMoY9Xp&QCGO{+zeUcjp4KR>xhVE~4FSyKsi;hRG*vnv1`IOT)R%Kw zJm%?FHwNVo-2yBAySoLbSN~UekN83;Bvk+X^!zA3q|K8z{2!<1idtFI69GT`)Aa0| z=6PMHA&>+TQ|d*jd}qrCciqhhYB3QX{aZd%>R0_=uTrR#qHaK?;Gbt_sjUFG4ynVW zStKPUz&MMb&RW70C+SxGi}r)~bHT^~#E=dFmvkJ!LaH|=N>c()A7d4nyQWrfZBEp@ zsa&}ibs*j0+sg`mwEY7kZ21w4@YBDSL6O31u=Gel(*b2t;?jbr+>wBK-p2Ox0xD~e z4HVn7ijEuXT!8RRyXC;K#JrgpWM?a&bXWrF`ZN!Ypzis1H1$fSoeP^ZugVfONQ|+t zcW5@Io75o@$2Vs6Mg;f-1sORE!p>cEn-Ys!RcA+U!VCADVA1A3adb8Kecq>dFbc@^_V&gkB;1dTjD!t@$EvSd z#lgUQg26*XiH@6S&x(};*RzwQ=K^=!jQDPJLhBnE$~ur*l||w61r#W&F=);qE@)7h zvps_F@IJ9c0YQo=BqK^7YU`I2?~u8D_X5tdxDo7oefyJ}KJj zKKQKr*6%TVB>m=dpGRT|LP+shs_UPGt%TZ)bh9SaP!WnionUl)Mn*>MEd!pp`T72V z;SV3QL=MiMa!5##j@F_nd4r0m5!JrQJ7RI4kwTdqBa#Nd2Vm>WKytZdtHTsDEva?8 zPN$H@Cw`8-d-wiw^Cmkxdr5={JG&%63bnYc9J_d*Qx!()9|5H7X4Uiq5ef#a2&ZvF zoIFT&*`IX?6gkgo5*aPV<16@zjftW;c8#nPopgV=8F$Fl1fmYz<)ECiH z|L4BTW8&ia_#QhnvW>$$5sPVE{Q7kTi6jmNboYIKXX^Y)=sXhnqK_duS3Lc@Hq7)E zqZV8}Q7#fbe~GqO-amYCGy-cG8BU5o@BH%mL$;@e`=yWnlwd*eHykoTtlnCakKDOp zdM4!zz!+W`8^K6DU9JmNeDhC~Y)@Bm&_A?7P|0yy6dEcxSzp)PT!lFcv1lV`b;0)| z6-%A~eDaqqj7^11ru%-o${}ScxH+`^j?HpC3E&%=Kd_$vVI^_4gS$hImK^l`KhOKM z=%F;-_CIJjr^kPg_v!87CYy!`kp5W|bIPlq}A{W9V6H7oQA@+DHCpp`{eCJ znb4gVJx`#-;c=dsF^6x~OjEa=DS6N-#k4x@UH6~lW zwkw!+rkmC1mL^RJhDC)LOM~u4$p`A+>IrlTSagH5Sd$QVw-s%I;fk%UAJKN}y2WKWVjqModu9c)Q!)KG+UwKo&`- zs3AOy%t1OGX^o`-lArBhGS6{E<$)-``M`lfH3-|eF-+_^^ zzej8_69s5vBghoN#Pz#mqSOs`I5(U-^$;Ba7meJ&zkrS&H#MnARO0MJK!vHQprHrt zq-SB5fF9bgb1L%M*h|r1(3i{}6@E3K+lD3C#V?ETQ+#?7!O{^#STmL8-$N&u!qkQTN z+=JY<`QGpeE_FH~KTG(Ez3$#4CLbYrPglV)Lh8~b1J>@#osJ92oMi8e;u6jxmG_?>kwi-e2E7YtWySnjWQk!bCnz;h^RN z)0`$Q{{v3xolTq)i>lTdi>ey76Q0g)vy)Ej4D+gK-BZK9ts7(J+yxeQEj!VsXTfJX zMZ9T~hg+JtJDcNQ-RT}eQkML)DHBDss4^?k3o_sY*EK-X?rq!JTpbM3>V4)<$t__w z`Pf{P?iHCYZWoX6OJk$x$=c*z9{vuvq%2T2@^qQaIYs&I;ab1!WpgLw9ZMF{R`pYi zq%FFQM~4~!p7_^boB`ybca}CsG!?3Q0|I%OZ~2OWT5cQhTIYXlwwiznF+S|HsHzUH zkxg#IZgzHgJQyojjU)3ro9`9k_lMYoi&&J};!r z;dv{GpILFs%@<0HH9U1{RVfqx+1mda>Z-n+OxjT>lP|xnO-h@vRA#u;U;7Pw?|$_g zSWdtBh#OPDv7{#p3o|@2QZ@KGmxPSybQtMch1uqic9HwL?BTrFfq~0L==gg^PUm@Z zdMm}bXgfL{3IyPkVc|0l1pW=rwCIU}JOd-vU-=DCj!C^-e_+MRl?HFTcfj?TJbBC2 zwGyp~<4qRjy|6ea28%lQ+`01@afRHbod$A-iBc&WTR zHEQFZLwEo8B?->mJUj~UO%Qn)R$Rn*o#ERFdr*cYEx4#jJ%SacADz&tc?&a!0irHi zdMN9{LHMOUcX4NBC;p|{Qr3k)n#H9$E7TZ_d^DZyc z5tDX3ymACQL;;02%W{#ak1uYmi|2PO>Dwx4WtB==9%(o4J2cpkDu4d`BN(3^0T`C) zZ1p^_vA55OPi*q0*xX^_Vv>+Z(zgj4DX`+6oH_*hNjej6-~O92$Hs}?suL)ksj!Rp z^E1q5=1$bAGVq%|O29{G3;0MH+hFrhVN-ck;Wzw@fp~uC>3v}}BMNp(O3Hou>Sx@O z7LDjhJnp{vF??y!(Emt#Ej-bxyOe0Ee4H$J&>;LHjP=dtC7hfaSsGUAmSn%lWU}E% z@V2;X(kK>b5w6yFZ*;)Ly$iBmR?kCwOt4+V0h!)x9RLY1ZkW%KV||nA$tABG&!4;W zN;SQ=lm1L?FOditv3bODP`4Mvf8|Nlk)!b|>6u-)}$u5t|WM zyXQEP1~=@7UCKimW~+L zhl^FoFnz4yIDLuQXea6^`0!6f50@7~s1E=i7=Q*4z5(K7g+_;d5dL_cwGtCvU65u_ zeoK4^cOhv5-^zQr0{-(*`=QKN5y7GE;#o7T(jJ%+CxRU&NmG|$ZqLb%YleEo!`Fy1 zK^$+S>^V`VX^44{HBz)yb-%%^;XqHC!>l#g#{q_iD8l3yrxq82k(|TP7@2C`IayPI z=lP|XlXQ4M+|&oj_djsSwKh?+_4ysa`$dk%;7+&5gMIr}iZnAW_AJh^5EBF0mEglk zCK+w9Yz9cp{;0@0PQ+@ap1G*U9=LyaXOH#Uqzt(-={%#wz_NShMvBGBiM;on04lQP zCmXGZ*9c|C8DJN}j7l0C8yg~@Y45N2$#KL@b2GZ?c#Zg4u7twLYAGaM%?)bi!n%Fe zu2-jOKGIR}e%4j&bK%0;f$bDrwD^aV)^-$k z2%rr6{rdG6uR^##4zJb-y2W44JlyU@Qx=<)Pb2!dZFwk60!3t1Z5O^-^S*o#t_fZ& zC!?w~H0;4i&&;_oG}Q3*)Xc$hVeKXqX?cx^#n!#iI6HpDc8klj(@lMq*kM5%d9$0Q z^BBWeJ*HUM70^s~xZR&^<~Kd>Re(eo`soR6?61p2OI0hzwyCjH1L{|I!FdRsPI)k$Ps-AQzk()Yl-A)%p6)|bit z7o@!uzSyRS<_ATiz@ISF9U*n7YS?gy)f4JUO5U)P_8qFLm&{XQdXcCk=n);NhSgEB z-hn2nuy`xjt5L_oX=ykW2V5t$QkL-y;^f!at;maM@ zIt;M_{D)^zS4749BxY|7q4HqqsmEW_J!T75g@ScDke|%lEiR9VtcBe|(s0!Ywcpk2 zBhc2)1pe+RX*ew_5R3>x_9&@sEs1I+_DeEd?LPdnX@F^UU|?)#b7zz1vwP;oqeeX? zy@3LC*7t;G(fMx&&`mbTAj%hUY=hn#84FeUC0$85*JcDK!!0h2oZav-2Uz1D^oZ26 z#S|S}OI*NYEnuUQof5U5`xzR6lT}2DVwKa4gC>WXV-fdj?8-5?n6(Deqg|w0gR7fB zb~P;aVS}tlR-L?uxHm}zD0B7|Z8h;hHxEu1B{ttHV~1@`b_mq9mU{894x#dAYik=# z#%G3qbgB@Y=sZ|lBI6ZDiTj4NXKM)}Tjz%lkI~>BEx-0do&kAn@S5s$W*}WY*Qg;3Lq^?G$I{X=wCH(a;Ktv(m^}wZUa{2qYC1+^w$<5{Zz-|{ z02P_;yxV|o;L`c~j>TJw%(2-Ge4%;^bJN~o*gm8eZZws=tC`59=LdV%YTw~%H8jIk zF~lHD8m6$R0jbSF)%E8$3rP)5cJGPw{~>+5SA15M3cc_rXXvG2XnT&Gw0V()(0+nU zp;Yk3B(qrbd@CWdfyirR0Dy_w`; zfrm{$uRY7k;}_N96lY}3;^yYu+IJr9?d>Wj+8Jxs$p}A&KrTaC{p`j+Y8k=3>&qIs zxsfboxLpWghS*|#+!otFFBZ7rW*qxyop(`NGaehXbWo4(7(_Zg!+?OFt z%4D}9*o!D+;{t?81uz#(OBoNB7lX*zD8vPcI?T8_FiHUX5WOG!Ul2Xh?uEK}Legc- zA7lyC%F`TK7=^5*&{ry-sFBMmPV76u7TC@X&h<8}6#}p{k=E8B!yqiP2A=Do(ICTpa$}hl_zCBrV1*rNu!;fgtwIh6Ep78#vFH znVA*mne{#@?aoN2{GXh`+BHTC-5;dO0ev4ht*9K3Wa%tHu~Sp&TU4>76iNk*#+UPF zY7$ehxd-N{vM+4P^ONcM3K|*@kzNM5wYR`mKn1mFCVA7O5!O4fe^@+g?(C3ZLel(< z8TLYAmbFY9sn$E$L<|(@MeCJ^%!%4b@rR;wW6Rs%H!ChHOHjCK7LnwR$E4{v9bP(h z80hn%ucOuuqfg+aC3K$ zn*huIt@lm-Y?&5tj2%p)#N+XBL&YKJ?r{^1!kGACv-RaZD3XIzGzuN#-=bSeB83G< zf)8M%kmKo|J(+3jUG*9o8d<=x{TR_|C5GAcirjC(OxDWC$ZQu8Y39emx&vJq{e5L* zXQbStv2{+1Q%+m47!SZnWw0JQPx->DuQE#z#&LA>xntBf&rS{V5r5To*a^UBleIf( z$2xXJY;_rP??T@#A7;2aiMmAyXdYPHB^i?qh91}=Lv@zFLkYB_CXysNP_^fQVaEelvBU1XJL8gT#`?<3!doftOXq-WtJdyf@1H#fgW6Hn-?BhPYIU`MC7 z->q7;A65nTU<~tdxC?TRX6=!%;v(HlspSa0RD^^LC-?wnP!q5r&oy1M!sAobn^xpn zFQ>V}@7An2NPY{D%79MEC${*79JcrpM(#L{F2@gpgGS>+&5eFMhWAT5+oL+BKX-Wb z&xswtHaCU)7&cCUK)b(Uacar{+E|NbtqyF�wh3Xutu$tvDz88GpcJ%%QkdM!j$5 z$2I`^{;1)p0