From 1e96ea46d43ac456e75cbec2d4c83df35e612ead Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Thu, 22 Aug 2024 12:56:24 -0400 Subject: [PATCH 1/8] Add routine for modern dark scaling --- chandra_aca/dark_model.py | 62 +++++++++++++++++++-- chandra_aca/tests/data/dark_scaled_img.pkl | Bin 0 -> 80177 bytes chandra_aca/tests/test_dark_model.py | 31 +++++++++++ 3 files changed, 89 insertions(+), 4 deletions(-) create mode 100644 chandra_aca/tests/data/dark_scaled_img.pkl diff --git a/chandra_aca/dark_model.py b/chandra_aca/dark_model.py index 0f30833..0d0453a 100644 --- a/chandra_aca/dark_model.py +++ b/chandra_aca/dark_model.py @@ -1,16 +1,18 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst """ -Routines related to the canonical Chandra ACA dark current model. +Routines related to the Chandra ACA dark current models. -The model is based on smoothed twice-broken power-law fits of -dark current histograms from Jan-2007 though Aug-2017. This analysis -was done entirely with dark current maps scaled to -14 C. +The canonical model for ACA dark current is a based on smoothed +twice-broken power-law fits of dark current histograms from Jan-2007 +though Aug-2017. This analysis was done entirely with dark current maps scaled to -14 C. See: /proj/sot/ska/analysis/dark_current_model/dark_model.ipynb and other files in that directory. Alternatively: http://nbviewer.ipython.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model.ipynb + +The get_img_scaled method in this file uses a more recent approach with per-pixel scaling. """ import warnings @@ -340,3 +342,55 @@ def synthetic_dark_image(date, t_ccd_ref=None): dark *= dark_temp_scale(-14, t_ccd_ref) return dark + + +def get_img_scaled(img: np.ndarray, t_ccd: float, t_ref: float): + """ + Get img taken at ``t_ccd`` scaled to reference temperature ``t_ref`` + + This uses the more modern approach to dark current scaling which does + per-pixel scaling instead of using dark_temp_scale(). + + Parameters + ---------- + img : ndarray + Dark current image + t_ccd : int or float + CCD temperature of the image + t_ref : int or float + Get the image scaled to this temperature + + Returns + ------- + ndarray + Dark current image scaled to ``t_ref`` + """ + + # Confirm t_ccd and t_ref are just floats + t_ccd = float(t_ccd) + t_ref = float(t_ref) + + # this comes from the simple fit to DC averages, with fixed T_CCD=265.15 + def get_dc_exponent(dc): + t = 265.15 + dc, t = np.broadcast_arrays(dc, t) + shape = dc.shape + t = np.atleast_1d(t) + dc = np.atleast_1d(dc).copy() + dc[np.isnan(dc)] = 20 + dc[(dc < 20)] = 20 + dc[(dc > 1e4)] = 1e4 + log_dc = np.log(dc) + a, b, c, d, e = [ + -4.88802057e00, + -1.66791619e-04, + -2.22596103e-01, + -2.45720364e-03, + 1.90718453e-01, + ] + y = log_dc - e * t + return (a + b * t + c * y + d * y**2).reshape(shape) + + + return img * np.exp(get_dc_exponent(img) * (t_ref - t_ccd)) + diff --git a/chandra_aca/tests/data/dark_scaled_img.pkl b/chandra_aca/tests/data/dark_scaled_img.pkl new file mode 100644 index 0000000000000000000000000000000000000000..ee8e881f1becd7180f222b3dfab310f865450726 GIT binary patch literal 80177 zcmX6^cRZDE*iTak*(;Iir$kB;g|EHGvG?A4@4a_2B8eihDMeO9BFRcBgi1rCX@~cD z-}n4+&L8LDKA-1)?(6!lD`Izc#*Y7gg=VEtcw1Y#*xHyoy4q)@FuM7=dISnsyL;IR zxca*II9hsnSq5gMP?>w#TD!Y>`*`_U`($NirR@2?3$s%8x!L^xC7D?&J7TkB$gHyB zv$V2!RmfCqRBYn19QLH_vGED?u>IdZ+G%$oD>32!zu%QYb}36;U0r#bc>Uj-$_}3_ zt3>ID3uHUK^qSPF+3Dkm_WKbHOKGsKs82iWmxsCd<@Yh?lCX6y;1un95_*b5rY}v$ z!?KDgl!DF~;fEfkN)D<+0OW5KYh$oVmbvaxXNeEq*5kXXl(4K2)b9Q~22@nM-Q!+v z;8|W-8CA5z_aBu?`YmdB(naSr#cz&+tDM@?8Zu}Z*s7Z2Pr#f_$^|APUyRO*U#+^4 z2!U+!49apXWYFlpYf_L$pvw5cwQ_k-oU3^E&0GOl-wTfG+eYEImfv7|Km-b&9g%Ky zj6!zjAF?4?3qq7LQwA3E!N2ZP{l-k%}H_xS&ZI|c0*hI+^P3Du+v=?E^DsL?C{NpFKpYLt3#@M#eeOe1s z%q@G}()|(SIFq`T=m=8C3n9vqc=V5{za41OfJZw+o9QVFq;u|(6uobYhswb^$6uI3 zbS0+fAhQX=)gnG^j2I$x+hKDeHV7RLic^$dWMh#`^`R)SKlKfUp2~U#h}in$S)u5R zt8Y{M&Sz;MxR57nC{qnHaY0(=9$290Y4|0dIeT2qTe&|{sSb~`cLLMHB(UDLcSn7U zC8}QOM_)OtkA*$eRs!RC_%i7}zbX-d*69BAn6s|vZF8l$dC&z8oeOQn?#eLFzo3xt z)fk1JN&l@w5zD(8PfYeR82+vhWc5?RMKx9$7a7H*TMmMtpdL zc_ODj3^jDRBz$79%-^JZ<8U}qCl^*~-$$Y&D{WJ@$rJl&ytukXNXVsQ*3&MNhBr&Y zinOvAL>a3R_IhZbitNIkl9G$aIQZ*~@30;AF{&LuIueP|>4Q}x`!(ruMgjBrLz29t# zce-{_TY5TRmzDH*`aT$v+x_iBDh`-lDte>#$sGA3H1#Zn&bULG-*HsQ2Iemxk;Kx~ z5M#^nlbS*Up4%pQmun+&|L@V{5645G<38Q3z zi!lzY=bSsnA&;FGSRzY51))vTFw8?<0g0o&)%D`x2!AmbdOJ-WbdtN6&YW_D)ZCzw zx^)097NzsBM_R*QwO=k|#2ap5;ahivl;Pbav%S8liOCF7YR&VgcNbLC`>We$AaY^Y37Kr&wo%3@k32zoF zX(U@nSbx)RVJVb?7sflxn0j?_>vxEnZUywct6{qUKt zJww!98;jI=>$%dl2zN_NiJ($L=MK8j$pIbsDcK-z_aDrgeS{-%G-TJA5K9s3A$t zKZTMj8(%tkxY-?bvG}T&%I8p<1=_|*4q+$1;;BF85Qxaz|~kaSP>Vb>$M;6 z^~2i>$sCUgZ806YIeckR8~SeMa?;+B7zn)o#@bI6tS?kDmc^|xFP^vOd4dA^Jg*%4 zZtse+N6qiD%Pp{Q@Xl%3^Pb2nx0^ZgQx%e*^Il}!4}|M%TItKN2-vs|KH~ptj~=aB zOK%lp*hs0G{;p9*UA(Y^4xcR8Wd_piDdnIwF&W;SX@Ll)43?ThuArf_7k^l3hk_~Q zP#G_Oyx!W_U!<=GQ{B_vx?1vRwmDWR6y^#^mftm}l?@=rn8$tegB%#|t5LG^dVu=* zqXoGjBI4uZCL^Towb4^?078)VFMk|9+wl=gaRN>YP%vFW&9Ng^^m&1k3@9I`nelnXLMg)R zBGqmUNTl7JZ)%pogrS0k8=n`<&wcgoUrB~;U8vI;0(V>xl#rFk)Pr_+Ax=#b{eO^q za6!`>UCDiiudP@?Jny1B4XY|vBl~xkREmI-wKBsaK^Dg0RH7w4%6O>oQ~sNn2hK)H zScL9zgQ0-3HkElG&eFa*WiJwjF?*97{zP?5HB^NO5zpJrN$iJBc$Wbg~X49Vzf`VcDhap7Jq5pzw);O#nRksryx1ZeKWWwZ6^)U_A{S{Lp-43 z#n8a9q=VrA^G7e9SYzv}Ug@_6N^vNs+4(cAeq;yUY4&60qdAVW-~?1 zws7qi8`g*5mu^O>+vYI1T~|J4tbv;4$!DVJ`be!iwWoB%3TlsPvXTY^up?FPg)#9v zBAS`ME2V0qnn&rgS-(8If5fSsw~&J}Z^2yB%?#W=HHQ8yeUy0JY4PSJVcxykQ|^co zew&?`U;i15!w(`iz6qFv%_o88zZ(wtptAmo|8p{^={1dAf9W9qy`_NTRs<%7wgW>$ z^g*qm^5wsurns}6rL2}@j`NIhBtu$vta#gh^0*xjZdLgN+%~|fe&Z$9W=l*b(5MPn zI>EARFSxH%)5wGfpJYxr<5&+hr_h7vZ` z$RCOZsQJD(sxVE}y^79~q0af7}(+#YF#m{`_moum{XyHbuh%9g%$6 zGLO+*6RKR9r<{e9kfvxp>uu}|y%N2ACRHCS(l_i+lCwkMQ1#&uBY%tpXV)z*SzvU} zpUVvk>i9<?nO6kE;Tz#dX;WTRXz@Wf zzh+LhCkehs1lH{|tl?WZ=x@EHg1q&gV^4_vs-Wnr570G)USfHik*pWiOG+il&q~44 zB8dH$vlE_(b%Y;2pbdIC>28MxB^=R_Afq`{6|a`8$$58*|hl5 zHAnn;XSnGfr;O&YshIn=%8*=ad!KUK3ldWuVPTkT)2ZbOp{xEW z^3~~R9xXdwHR_08KYZN^9_wK-UZq3!2vM&VZe)1HD?yHXJG{+Z8l$UAoURK#SZ%NP z!WrfT+Xb@sD<+z#_-N!{`!52zTTMkTiF4C=q$>I|$rOKAhMz0f68)5Bglyd)4t*nh z&hzh7aQAB{$A<}7#6CQCH?htULQ^d{>DTOFU~>9|_>c??awtwP>Jq$*{ZG~D9YDBO zifa~04(*&~cLR#daq@=FaJ8`_THO5}WdBk?*T+LQB83F^Y8&`L_hNdscmH@i6EGLXx-M%#5;7iEpZRrgBe5pug>+k;6F;ya5b$(~uD zVDcgtYp^W@R-P_DoAiQE&7_EwunmrzhZ;!jb;X$a_VX4gKj`E){=2ImiPw3{hd!OQ zgVB)PcK9V1yx0}%$s6N_kw=ei-HlYjee_|;H zG-84Vm6}Kus$e7*{dN`_AaH!OkZ92rFK`teGJ5dN6vfk`qm%V=*gK(2R`0HeYA14e zaaViDR`>qQGmFN`4w1#CL}zSwo^6aC3jpt;07a;pJ}QD9vp7jQqW()a=fnM~a1~md zl;w4T&d$8DEDHxTI@Z+(uFHT*DofF{?kbv`KfUL!i$qy}=UEb6F@O|JNM&J#t~PPbgYtqA6Y2pkn9 z&nyGf0TAA^so$U&0q(M;#v)f8_#O>z*VyTgXddq0+HbU>^gxe=@1HEnr|Tt>?E?{d zsyyCL(*vH`trZnij`%tM^vu=V02sZa6}DrDL!n^7FGo`kd}MJy@leqW_nzCRAGqfR z&6>kE{)T1~crHH5>vIg)!$Pv^QoNyXui*N`J|$e1%55xMRf5>;)k{vx_}wwQyP|J9pU^2vDL!^A7-3-j@+&@S zb)4xGIIVEa02Q9|J%am5(09*_7t4}D&3W_fI+0KuW!+b^w@d^7s6Hm@4_V@MVnb~J znGRf^Es;I{lLm{dX09MB5!o=h#CKcm?X`PD{!`8&}4aUKHyDBKFFV&oRSD)Lg9o=Xll`y!j)=CL=LenXw-C zSq%bP?~w;B86L1(>)`uN@ICE+l*eNdwZSbS?l^kN1wYrtty7(JF>G+@>c#00+~HB@ z=uvP&boj&8@U{>bV=s~0fU%?!n~%dZeSLj!RQbBA5^?T@c0UyCv9^aMzy86u#PdMc&^-}&&lV;b z$H06)1HonWc|TS|Az$QPrgTgja$K-1LJ%3o~iCJgl%SN-GH<&hT}JHG+3shuYrvv!BiW6WJP}cj?zTL6r)!;jPT1W zS}lR74}yv~<+N8BM=f&#ztbz>!6`OD)iDJ;5O$V3b2uFjT?EVL$Nb@5QIVm&nuG@) zj$*V5M))pUo_~$d4+LAbU$wYu;a6Z~ocwhalyz&${S4K_3wBNkmml(AQ;@e$pz#Dd z8BHGVv;%t5O@_Z7^MQcrhLA$8GBmyqHJ!ZWg_mjxV+R_{u$V5-9dtejHZ9BzWvW_m zZ4V*k+)+cL=PqkIk!YmyF~qhIy3Zem&>W6AOLVsFNf93Kz++VtQO?I{NU(XQd+1Ia z>L?s2_w~iV?X4jf%MA~RR2$GqX7B+Fkgp!iTp66B|+t;hVA+C(+@-|-W`e?)o1r2b7YeV_%Evqkw ze6d^L!Zd$Y7?!oSvs`ouJ@qUTtxJYB0=ybbO9>sX_3b{2(E1z5P*B~}k1++MYv%Vm zZQ7U+$|_SJbTgUwya8(ZHCNoE`TXq%w*q)d4bFU!*24y;;VIERZ#1h-Fkd`k0j}(u zo5GIq&@f7lom{noiFiw%Y<>v1)*|9mm8@{mRYf@MngUo$7b_I+8sNQ~WvJS;H&RU3 zd<0tpv451hmg83(&T#gh-{q+TQ)B;sj4PISUqrv}(q|1!J^qyvamp1NF1_24m74ev z(-WERCJiqoV|}mPL9n!)5w?+xK)ozIiTtJqlJ~P5X|Obf;5hT*-t|lzjPe_43bDYJ zfAA}=ug>sT-xcE7r~=-H(-!7;#KAB=v*+JuBSiFUcs#I^M*=H}Ia*E{C%wBj42m5P z8FMYgWtSaJYjJBo{h|xpbe&OlvctFQa(|MId_ZkL@$-j?7sPLBUG*Grz<%W#{gE*Z z>>57KvcRVdZYItA?c!*-OJrPJV9-Pyt()j;VQ)wmH)QWfkU>+_AJtEOX1L-q^kQW@ z0>8>b$>J{>A@;B{ZNYVQw2D`*X1q7Rl)*Y#_^1(__(ft~6?@^0uzYE>z5$k2!<|oQ z`(on84jE?aU`QOk66sQ6ft0eWOT7--z=hWznx(?=faPgmhDaua->d1rp3_0bjU)Mo zavb6MI8uS%BnIDZHSgq3NWxXtVBdX*jQu;%|Qpk=XL- zQO8je`1)zQ)kq8gUv2DIVKLF4AC0csrs+U6o`bV&I~(L{t~BL5o+$k{{N%}uKL)CM zLN6R4?t6H`!z^wka8Vj?yNdfDl-Bgv!&}N|`c7TaTdxnrsVq+2SZgrN)~IrC2Z3cn z`PkNfnuwB`H)Z)u!q8UrFrSz*LcfhEh0`aY@5h0jNA@8|j%HUrOy-E)_ls|F_$gwj z{h8dI<}m1aD0>GFXd_}3 zW}F`n%I=d#YI2qU<7F+d#uo*-^`wDKGd{eF&j999EzCV-5%_mbwPWe!6`Xp}$>GUr zj^l>@d)dUSaQI!s36?Y^bgv{mS|Rj9?;EC~ z+z>9tZ-du5#IUbmkGN942Y6Vwb5GyVL_5Os?FRf2@_gvo!mt6fgrih0)obJN`6oZ# ziU+`y+AB!;xH3F;mMONMwgP=?FSXR9CQd$Cu$CG11oew!A{PbxVHy0KJ8aVf-8LK7 zT=OEpSj@xGQlN_x<4Be3xvpRq?6=(Qs}CQG=sBklLKo91VSDpU8)5Cd6M>o+QL4Bz zr1Du5?o`U}WmeNfmyAvOUwJizIZqs~KPrd#1x+8bC0m4*53!ue3PstQ^wXB?7C61{ z^@Wd?mPn$qu3UzM1Or)pJ+Yc8lD8UXCt>&PwHUnpR{U1geGoTxcm-dmq$`ONq^Bl64WD+daHy_ z^ChXW-~^d9!prBLgglMIH@+E5p*CHd9f(-rA2kJ=NN+>urZ-MTQZm{XWg|tQG(pKd z7_$W`S1Rvk9sM;70=4W8E~7elwO-+~{g#Br+xzKC)?^5N%CNu@ zsR%AoLxbThU2w5I)!XM2kGdHa_nujQRD~|gEnW17_lx4;JC6bp_t*CP=SXdgTx~cd zNUYeSIa5s)wFmWH zYkXD#zhCYnvkoa7t1Na~Zk9*uRXGmc}F6S}kj^NHl1c7si6$DQeT{Ydj`C$l6?F(}aui*58mCB{aX> z(7Vm+h1BOg3MubhkmEQM7$D5y2HX@!DUDEG=hwcgO9j4mA_{W_E_kD0 z5T*543;x9=pQ#HHVVx6k?3bzy+TRWh_disJ*TitHB+PMRfs1b>(Fsn?f=}pP0xVWt z9lkoTczp1mm#{@J&Uh3rGkLpUppdMU%ftgc_rss(B@#M%W-+a4o;Cipzdc8569KjU zo*k&O$L<|!%mvEIP;>9SlD{qw`s8e2q z1Ys7xGs|@b&FDVzFbQDp+#_&Fk!vlSjdKzUB0ie!OpDtN& zhTK4r(O5zho^&uQoKQAFbK-#$^?^xaxQ_b+NsV4qKw8uws{Curc^NBXV!m z`Xln%#>0?(j`*FQE}k!gr*y}Mr%a@*(jFPq7M`UsuZzF$`8zI{ut zj$IbBf!?AHezG|B!hdmWEexMJ$~@v(<1o-jpVV$Z_}}O2cZ@SDW5uDF(y>JggFjTM zn9IGvD4bcATBw8&>z5mc&Bkk8Y7gp<3H1pgc6JiR$6V_J&-*4(5l9%!KYPTLvlb>3uy5);)#v7lQ!n7gY zQ26DMuq+Brjg5W%XM|+WJw^IWx-g2_Q8{Vnj;25!MJ+ZG_Kv8w_Ed!ty29@>E4IHucLI9?A39LRaL^o$sTH_t#YDsq7Yc3=b3Ho2dDYCi?l%+2rsB!{3-5@ zHczkOAu3mhu!iLowYp&;nNt189c9G-(pIKxOUBI4d7kW)OvvB(IqUao|d$vJqaWFxGip(8BA};XU+( zuOOEm+rla12{->OS<<;{Sn*}uRZC$7(fqLEp_55SN4?`jhwly=jj7HojjD z(t5!nxb|;Pg8`n0{Z%+X@F3U1UtP%3$wRdEpB=nP?6@d%5a<1WuZ4{b#EehV>s^ zw8!nOF>S#;osr~&4_~Wu4I8}S9(+5iZYK$RcOn_nuULb(VB_+RfoM32yc#zYjYQTy zPsbzzzubG~vwHlW3nKb2eU{=;gwyFv*%)DQ)boa!1~y9({E1S4w@M9`Oef|JRRSh- z`x?ldafenRvUgJhyUWal{^%KEzh=C+HJ>ze7}MANE)u2I znj-8Cm=sK0IY6U?>E(s%KV#)! z6Kj*FeM=Jl42Qe;t~lUT_Ytj54SDz^%8~OBID4GA=DMRpJmNN6Es`F1fa~GZ-=hcB z;n$Nz3eF|)L43W-6K747o2{PRwc~&Fd0}j=#zo|pE|3+vz7=_txBHX z3?258;hQu$EpldbfH%?c*szK=;Qks?Ami}(U@uc#-UUbo$qMm zcc&}Clj()sR~+p2)_CG_ za4+@ppb08=otJMq6N{#uf6Dh|1)!hiV-fkOC~UkLi+?+!jbC){b_fb8V84BYP{S~g zaq6GBa=r`R4;YSaUh~9~Q0hvczn00!qyD$EZR` zI3BA+c|J!Cm#U~vrd=WYC*F9i@=Q(KeN<@uX~q@nES5*vB(V;m9n*P1@t3?PaP5FSWLCcLdDlfFy{o&xWR}?Hy`IeB zNx?AjYICjnY7QZ;hKC;aWuTg_RO@CX4(|y`v7-kqz{IH|?M3v2}!s45D*TdzOhz$AuPV%P}+#2&Te^>_LLehbj z(*#bV7OdUlvqJ|aA32x4gc172N&PQ9Jn{IT8HO_r;8tlfZ3P zs$>11L_?$U1hY)G7y5*SF{eeQTrh;) zjmnt|?TSdbY3bqPD*;m-J&qxT2u$i)%kJ{_K$gT)^0oJbAEfoTxxO(Hoo7A!-&vU= z?zxiLK0Q6~47M&WOKRi2BZtxIO$G-rHZ7EFSB3jt8~dIitnVhvlGC93-d@_YR)* zz~5>=gX5nJ;X2a(tAoo5f0l0*z1Gsjx8KK6_0bYW_v7P?rBn#L;IRMGCu&GYFSDPJ zwMM<5U1q$X9VEM@)I~$3z#~YP{e2|_!*_H)z2w!xX$lj`TfRijz=G|}tv)SCzUr;D z)6qxyN1bc5ir(PJsvhYix#1zr?9SZeYk22i)tGwK2+ux9Rt-@pptxWw{?svhy#Dg# zRi#NVd<;LHZ8&Ct9cRU-I*Z*=+BNF?v_&5cr|0aXW{JEcZ|k(Uk_Qg$t#VB$@x%=_ z6r}9&MN!zFw~wtfv7G!;M@BgwDdmSLR>-aJEof}Nt3U{zQFHj4%4%cIfRmFePZ5)p z%=I6~i5%Zt#eYjXbulEC6|2yuigk^Keb=)r@TPk|LkXiSUJPG(`24pm9Byd3?|7hq z)rK(ZWG7iz>uVLn6MlWUgVI3Lm=(_c*;sd9(m;V_Zi?n-X_V|+`ZxAe0pI8as`f4s z`YPo?3yr(hFh5&!^F#p&S0bEZFUh2$Vk`Vcs<95v;`bKi@P8$Igz;x_CrUaM;WPbh zDMMEqClfxOIsP&R$1IN6buW{k6xF|FNa*sPQ_}9m2uFcbTTaJesD|F?|9<+l>0_^g zKN&aq4_ias?*&9{Lm5lurh zbF-Cf61hKjhW71v7HRyKMa8M_orOEA6QvEF1{inzK$@E+^qASGaiIWte0%n-y|~g4 z3o%8E-;ATd_~t+uKcUYGrL9_gnh5~;qnyCAQf?5YUh?hzlZp2FDiw7edsI@66y^T4 zgHGDhy?xBmc+v7<>BrGjjJ|2L)^{>Q!B|lFOTri1@nl5I#>)!zB0&R;9>B?$3|qAu z8mMC_>S7s+0+|hSx$rR^jD%cBII?AhX!8yqy+bM(^Eh65WY8Z=O~>n#=+xl3qy0zR zg9s$+b&&bGB;t!0^JscnG!_}Xm-P#@@F)P`z5)nB?M{yu)@MmvLk$}l!J7GhG4nAbTxg}7B@{>H}bxlAx?5bPxrkR za!2KC&aQjl{lTHR-;;J&G%+WW>>zmOBANQoKMkY>KNWEQr-NU9Te>0V+;A!6rOmTD zwqQ-nleIK*z{>fa8$IoQ`2IB@>haD%yq&B!6S=K~F707yo?aEKT279fQPRYqTq@

rK8Z~_R^@d-V+pj8lf*vx`F)Vphs6A!2@%jG?^zGf-a}Q zgVrSi)2TZ5!^YfTvCrw_%Xn8bZi+uksq(-m4r#0x5c>5Uiq2bGG2p1sI74|e1_p=s z_KK7HV>o8ZZ|1ODK^N1pXn@DE_6eD5dux>_D%>Y#fXUBN^xA6(tsbVB;=o!@mWHfHN^I^l@07q>k z2h6Sce;O!Lfc5jdbE#G&hC`uHCv6>;jt1eNESoTj2G@ z_Z>Z7J@NHU)1WGi7xBIhYuTOh%I>_eFw5<_&Rki#g&B<|Jh;if=iIgbfx8mzZ$KFJv zX0X-jtzal*C-a6sY-`}0TNelGf+O-8XqrtE^{`1sHP0X>h0k8Tv+SpI&^)hKDKf2v z&eJ0mDXA7Hcb^C-%Jo3oYOecv7fbNh|GVH{5skF(_a1L7$HL$dsWY=x3j_JjEeB`h zFxhQF%bgR1;s-8aoP>^-;g-nPKIa3U%S^+hJ`%LsSm)yiT-wODP7`bFgG)4L5}z?z zZ(}zSf?^H*JAH8`I#X5`xHM1p>qxS$GX=zC}Gl3l6&4) z60|!%{azsQ$deT3zbx)EMri4UBCTc@T-DH+I!aQ)+`zs0#w(h*CYYHJu;z>-LpPqy ziP_-PQ?8{hO#}EJSDWFS(L{FRhZsK{5;XIS^qXB%aO)+jr1x%rT=5^2NOCsD&7%J9 zp=2NIwU!%VO;Cf2XQ!{GnIX1CZ!TWDt$@EKFBgZ2xdK~?jjDe<+Ng32lI<*T!Y#(a zE0QgiSn#$~j8arW@t+I#-_?a8yzK`~MTZxZ(+>dUufbeQ#>+u?ULkpPUvV< z6MM6BiTsENN#O3pHPF)Cd$y<34LKs&BPotf_&oT*=HnxEc%Fa0ovLO5y+=goVNDJE zLsr$!gg>29m=NicAq{nUrMXIha6~xTKFUcqK+==7tG?fLp=0V^7-+AGShBA)B-#Mj zBy)adpVGvSM>a=V&k;Drp*rUAZZkXx@aUa;76^)}&)!)g{z%_A(s7IM*~V=Ib3WeI z#K^b67>~cDxOM*e_RkJ$q)J7bE2(Ip?Lca!m`(OyFzFuAwKB$E+2Brpf zLF!;_*{!A0r$6{g6!dLoUYwkFMGoz=AhBM3OqvS4ImYIJilg*0 zG{X_N`+d{Ja+f1g-0t2lCvt@;{$wh4d%f{>SEGBYdK5nIhNG;tDHNm^7+f-Zu--1^Bk8-iQ)|kkh-KFSB_jE&~dSQ0@X%!USRcp9+QWsj{Cj+;4+ToRT z*R@+j9-X>miD7`q&ouXaaXR$H83#nKYcLSJ*1+#k_rB+0s1i%-cdLyB@1F-3%r`wD zJrTrr;ZzVV^>lucqz*==g~PrtdTi#3 zij@EM73_Z7jLUWg$7a#=BoC z0agz))siL5aa8=C(Mp9R!2=`n=^L!z8S|oB>WV(d=O1%=X#xp61~fE}rO+ptGhjrV z@6(^T_nej`@ZCD6;M%SzJgsa|p72n^X#DsW`x=To1IL`RFWR76XXU>9 zHv;dxcjyimCH%8+fhPtZ)bRJaAOD*tMrgP3>0=v{Ms-o|y+lm~l*QeCwYlw&GH-!C zp$1QII?TLsN>xNmp>WPPhbAcWMvgIsM`5ay|J~+Q5BM6zv)R`AU_eU!*u*`3xJxDb zv{Sm^yrTY{l{ZArHIx4Hsn9eW8hK1dy<&y6#8)4u(_)}kIrr&QwH-b@HnW;4i9|OQ zJL6KP4$`Ksp5^{d^wke3sxPH9vFob+NA>`A}CKE%4mDzwcF^9GH5iw{}7ezh6@_I=jgtI&ro`LCz8D zu7B&__&DQda9oXoMGSln`K<8~_-TD)HTuXU!aq4sNKVJ^jRWtJ`!h4O@FAqQ>t2H% z{;Q6!OOaIrF6}lw#1x0n{ok0meT{HdAm_xT`*JvEe|&8-*9m6?`TedjBtTqOy_(g* z6XNz$rXd09Xp_1<5#270m5Pbxl4==nsyK-V6j~rquzu?1R~cwlFdiD^_rdsE&ns3f z(P&U7H5S)cA@5;hczA38fesFljT#%`=LOv&nw$_=RG6F|752cw_(~|pdria$7v#U; ztAy@Ku0FF}6sAk2ZcKMxfy=3T#k42v!4SNBJeJTmC69D&{``>!A)%|!oQ1s*KlRB* zl+fd5_BeIXP^sd(w^ZQzrYvarmv#M?{1Dr;W^~TS7k~C#AzvMFz_qmq#qnBVuJUNJ zG3PID6fhY6y8Bfhn)%BPjC*x)=Y8&M7^fdL+d`9am3Zx zWnYrtMWh!M#%$EZ;tJ)R5fYaNVzk>YF;KdK)7l{OE59auuXB}<>-*y9QjFYuSqR~$ znSYUHPKU~3A7$H6Bs|-HI55*{gQ_Gm{`xu)|LDQ9f(OocU2J6*GOtg3&-1*Xi&$@* zH_)zI71{a{jVT|)F|qgNv#eu8|4)1!H@wRM>l}@y^Gp$lZhvT6{zn@|!w28pCH!ig z(FZ!&;i`DNdd7z#QJl~>T%Q|g$l*)h^}Qk8#<1`4H2-fs7AEb^)ZJB{aCo}IUhBRk z-u-JW%giNos6DGUT=wN&L2e;E7^FgYt-42S;$b=cP$ z9T@iw-@W-;3$x5$rC1(#!)3eEr2Sqzl)d7YUEb?pzr{fbhBJz=+uoxs0%{3~#z|@)Ol`S{Q9X0V&BKyfqr>5Z{=PXhAAJ?YYsXJ5npk62p~#(vnD=CQ zkV{QL;7gU#`+V2)ePI7QM4_(58dd|^f7youVI>uF@$Eq$oF7hfF*i|$wRFtXwzdNJ zl+&4>4%#5)D3{@7dVPFo+q*f=83vygzZ`uPQ*a4Cr7x>;LVI+2q*SFbDy!Kxdet;w z8MfZ^fX4{xt>j!DJ`p&b7Vq~6}N!2=UfU5R79 zNV)PmaML&i|0Z3S`y&kDO~0nyx9$VmGcQe!HYDQ+-78CTHyyBO(lyt2I$)45O_-KJ z69O>_`@4w!n4P`%Ss>7-)Z79wp|FHc!Zj?o=h_Pb2D1 zoN!Y_G0zL`?=-5|W7rcO7oLLo148cQRlzV#e|Sxb&?#QBW>UGZ2gAXP;=}PlB`nT6 z^(OB0!Uqzq)wFy(oV$W&xrx5vL`vNsMC7N~hc)>kLw#{Zv|s9PU?hJ2yC%afp$3Lr zc99cvDG1KlzxPMAA>KXkWjZsf3Tw_9#l-`P7%Q6lkTj@+*LNk$b81y_HF*B1RfH4* z(hrfb1UNv4*-2_j#}2+Fl@kvLzf{z!ma;&egkOeBqLBoi*Y-R3szv%LxP@lZ!_>VI z$q}%;IIV(gvFhzh)bfO%dGDnUksq6oFUlP!<~*($c6Iqbk3x-{uJrAZ zG!&&c(>uv(Vb3$0BVUY^@#l}E%d_Puymi@0ufUboeKM#IFP%t?oQ@Lz&3;>#^4~Kho8}y`-k`O~R)>`tQ$Ata>1K z?**!sPbLU!tddEI3`0V};nxnvdYCr!r@J`ehT@!zsbpPUbSq^%-E9~J3vG(BUMY7x zy~4tMd?F10X!5`H#)RRjFIp=1Ti|QW%GhUS3;3Lo8nFMSf^I>{7ZQbpZ|KXsc&Jnh zNje|nHZ44Hyv>zqBux>mNloA0vB=|)b@DFZS!vW?pPHkQ@`x5 zg`MBo6mHQ)z`pRv-uZWW7zvQRx@%M$vC2o5&JIYyRwMU`R8}Aa&!ovQzl*{brkWPN z{|sR1Zb$BIEd|blIh710=@6e-y>#d?F=w(q5Z^PL4#k6QDjG+L{A^(7<4JFVr~N$n zgy(cLnwzp`AKfuShkJ^Y%>W55OoJ0QZ^?qWh=zVE*bW@=)?#ghugGBU!!dWo8&%(p zG{RgRP~0EtRkC1-Cp`ZptCX1$I`a)fjSX3l30u8=^i&x;@9b*~4)O=>uj#^rPl)wi zQ;DwKP(ojqz25qHU*xI{fB4!#{vV3YGZ4$S4dYS~5{XcX zWR!}KWR%0+d+)vX-h1y6vNI|xWMs9GvI-^IG?2z$r4pt0dcXSgc%J*duJb&PheAC!s6zbzvfgnw`K_eHqHKviy(KeoUD z;Re=^p45@^e$3b*d8JY1S{&GZuko&!lJF zsI~CxWr6T!YjwEpVUMT2sE4Qa(p{-ZmehCfgH z(w0dcV)gY7Zk5^)s8;eRwI*wT@peISD$xhsepL5ZQs;tiCw<|3PXIE8mg>@oo|Me& zbkug)5HIWYp7IRNMo|Hs&+~J}kiGtS<2R|3hO5HUhHn`o!vAE5CE+gYd~%8MzzWf` z2OiuH4b#DmXZFXRd``j-N!RcED_+>`z&s~;)E8CUKWLf{+Tr=5s>2(p0T6g7`TS$E zKh_E@rb!Woj3Bp=eWiNvP3k|${V@i7S~5zx>)zl@=ehoq=wHD_0xac&Ibf+bvob&C ziIRfqS1M(yh*aaemYksijT^ofxr+S374~rNf0qKW^Q+8%i+??kogwneIF9&%O*M15 z8Uw-S!m4?L+y@Q8ai_SnEbyVIy)Wut8loaQ>Ns{7!s~p;3bB>I%fIr`!!{-Sb$B(D zz9SGZ<#9hZ#oA%oLuGPs*Z|CzN*ARKH9@;LT)Vc1)O&*zZAahhAYI$6{b9N{1WYRa z_Fgt4e?O#KD<2AH*CQ45ZYmh_EPnNC-34P(&s|lU-7%ayWU6z3_~UnEUu4Zx!^%bd zEk929VcvCp{fmR5SGTK31rcJS+CJ0xH)t@e8>EOhO zc+K@b8`v~uuU;o}8jQ9Ep(FlKPKEA-;c?=Jn-gs8wU31V%bMaMF?Yyq zl-pi!u|QssM`b_BRnHz{)bRCjL7|w=dn?i}IZ4vgY-5eZTT!E{cDJ>09@m){MKWRZ ztmntQsBEO=YJVF07KXs{opHVv1`w>saqeDFgO-nAM?Z@eBozj)2a@?}#apMxJ5RVE zc#5vx5kxQR(zMx1Z;RxYWuYsBfw)p(SzY`j9`Qa~?YIOz@mJEWZbynDMw&ckxVhv| zy8Zm>?z@3dXH|D|WF|Q(vwm*-I$J!VoA!v?rH78$J02EKy>O;I*w*x^2dI)#kDIG1 zql4m8{i2&WHq}^c>yp$UeXaVt;%9205-@f;Ekn-Z6?y?r!ri!|T>XbS+aH_HKBA}G zk&PeBos-k?<`~%Y;RxH!(@23NW8QOC~Yn5mIQrKvGYZ^te7#(-k)II z_azi6Yg_E5cj*#7N{blp$pGYrp4ay+Pek-PQ6<~k>PYij=Xju@i~ib8i5E+BNMACc zJtD0Dv93S*E;&jN+32^pbSezpHtf2|DP|~nE`IpEq7_oN(SF+?`rlUS-*fpQy0~`y zx7-yr2ckzwKjUHXg_Xb){#{aOc+Ace7Q|?V6OXzN+?))B)M4iHv(9?pG&}7Xn@7&k zeW!C4Lha!FP`6JtTOOYAQOAZhd7&e{%3<-NFE*YodC?oV;&GJmbbF2yo+dL4IuPH( zykLx`*Fq#%9}Wa_^7x@PWxOPULIt&()s%xR0#KKD?z{RfEzqgt%y~q}*k zCO+R135(>T3Ani~A~}D|S8dp}f<^-|rvYJ({*eY5@BEU?=}rk=8+4jpgDs;W}LF%vmJWxsd4H{Hg)oSte6wTr@)UPf1UUsYGP3`-$BoeOs2yVPtr<$UC(eSh9peD5i7kIS(z z@w8v7?sP%Pl^Y+}{^{U1$Dxe#COQyO_BR^tR)LZ zq28;_RQ8`Q+?DPa1o04FgxFab%ScV=L^Y=-3h01caYnXUIv)Ci;#)&Aqi|Mv>n-mV zHS}jo%70)pN9pSIgAyY?ggUZCCIfVQQQ4V- zrbyV?_?c}~7J6e$mM%1U7Ke_r|9=f4K=ei6bDG1DR*8z@TvrvyNS9N z?WfR%fBBt;%`^^#!xQp-rzM%U@46n|E=h6?zcyblc^#^pLG9TczL{26>8O_0;z*5IY-u+=(j=FBn>C3o3Mw zo9%gjl1~RMDKBYSrb4lE@q+011uamVx_QukE)%vX!!=?d2Dr1e)%L5o9HP7(k57Q) z17ZWN|03^SePlMUEZf$$C4D?d`@K^-60naFH{NBsNJs9eJE%5!(jtz-vW9#;3= z&o2-9>ZA)h2~WsQ%qjN!a3T5?D{o(0_JPj}3bs;F!c);Lv&(rShur;eec2m=K;9>x z%97PEv>ZLAewyg8Zi`!QFNC7NR`;D{St1&GQV-kP#NuL-H|I5kAmhSi{bnm2SaZC$ zSYeTckk>+Qudz15YZISLCkK$<@$@AXH$71OH?aR_uOSL=6b4`aD~tOkWfrl(^fZph&Q zG3Z1jDdI63u5G#<3$v>ERSgeSM6Kk{l#n`PA;tQ3?QuQ)-EAS~Meb=Xuf^|Qft56NWMRs_dM7i zUT0EXBJUv??wB_mCwZGJgX@Y)P7t1^pxrG!70>&5=I!%HC-IM}F{JNbc7Z z35MEszRMcrDbNVHjcaQ*sw2^9YS=&m_q?{wpw`@_@fBL@G_Z z``$RCf$HK?rWo;6?YoiHJL&@tX&NIZqDRK%YY${?4oAO@y+iDEL&)8|U87$ei4-|w zmWoe4WbXenJZ9`10{_H0&L0XvFh{{zu>dCo-Tg259mzFqnrLL(ay9}YwUQRn`KO^O zGWRNfmjYzd-mcRgi^N;Kk%E>?JJ?^*dy)|3hgm&k>5s)07+t-|9;KxNEd!0#nI|N7 z^}a!V7pZqXTk!-<)+pdH?J3`-IcMl*f|k+M8JoE)1NF*dG1Z@cJ$#&;Kc}lxa<-FP zXdks?7s=I7NN_3TG+N-SJo;VgqS@|c$j`(D=>76DN z@PYA-%CUE{xOq+b9BjRS3g@Y!v=#5dMV@!j9W5vNUUfpT zuJC9l>8owNO>BxWz`@cgj&DzNL0fFnRBWvQUFMaycS*f_tJy27s4W4Sd+xYw8h61- zS9d)!Oabp^`|5KQ8n7t)uU^K%AM%}_`KNYA;?JFmiP!i2@YGi#duyBuMvNbQy!|r* zJ*J)1M*{VTe|OiN&E-CrJ#+p7MM@+xrg!sJ-S@_Ix{z~bAKc(yta~W_mjRRaf$bmFhkWVdOFFa)jLwh52VI@r=k$;E$C3u6uW z{q+x2@pf2}g-zT8{98+!C!%!`nwG`2L_WXV&I?qN-!t&S+EDH8KYzT}e(Tu~Wdt!5 zZDG!5>iE2=I&N>i8`xCqt-FIYQA_#IY=;@~dvo8Nt*Y=xUO}9?hOQlGqPcfB?j$)y zC9UP_q(9SIJ(3@s7Yue9)`Q%K43RD1CluJ|55-9v;s1!=cU9J&tMgVER@hq@!*t1c z`k3*1{Hi+q4UPoMN-1FD*x%{ z?;DCMZLKuj|3aY2bhxAZj|U{qx`~;P-`}Tf=~F1-6>jLN)N~aZrGNXEDL-pi zM{F8U?a~IXvC*Ea(I`+dhzh9i>SO%P*&Do!HgNRr$-Y5&X4lp-HC7**;aw)n>X~p= z$ai!sRS-S>BFl}b^|@R;=2hI$08c#3=~2HXPx`U)UyVtn?zp_&tT2Zz8ENqg?$Z;& zxOe-%qv~WnBmCYXCZ)*^!VJ0FrwO;?(<^Ie9J@GATpsN{&>A%lM> z4v{{$wB0vrSslN7-Y-QnYC>Cz>22jJ1vu(>KC%05hi_+s_dOtfLze7nIs>l&9Inye z{Dsrj6*qXi zFSy&X(*`?Vw>=&_=LGS_u^?&dY`p6qGNV6C>g|PBMID=?asLP3uCWD2JPUI9__|*e zySB?uM6>C_b)Z9v@q!ujEH<^rUkbwTXQm3#tRx(CpYi$cj4rll=$#pN)PiRNMU;O~ z4n+7B4c=@cbFZDu?@kUmLGYoc-;pUFxL%c!YHCr4|YFFX^z2mrR*nt3~@MV`qNygTm#w%?|yjrCm6S5 z1q+&}HSz4imQ6=r+kL$OnH|X z8tdzZgcUSU{oiGdKm6vnKl+I8Akl-wY1aOXT+hYe-iV9!1tizn$Vq7$u7@(^P1KfR zHsJGgV&fouktTVqup$RC$yivS3CoB@3$+KGVA=okgXVkOZpz~0E@ykS1vix2s_1ur zOF`bQ7d*cz?7<)WUDf4vJjQ#QDf9^6vGoI$+w?v&l0VV!xL2=%F_-!&qr(;unYyHQ z;DrX*7@n5#I{LtRo0g|sRsf#w{JwaxDIQ^01`otJ>Z9AC;roQKGS;X^6+XBbL2v!i z)pN!{ShlddbN_u1ewP*PI2>+;dKIelrhkO1phgvQx=IV@7cCWr4iHXHO|^@{E#l7{ zvD`dNA&Ztw;d3*ad|)4v$F}fS6>kdOHoa+;!Bu&U8;RMvaOzw>S?!>RAG3Tr^p1$* z=Nq+8YwLQT{rszKzGD%@u>=jBDfK(by^4xy^1q85HmC=IXZO zz(kYz`K5{^9J|@6;(JXW3<65`NA-0;e`Mo1RirJh?2=Zw7h{Exgg>r!=X5Z#Z+xcl zvmORysdvej24Yy{;4l3=UAT3;wjL8T#f!p=O49pWiH>f1PQ_FnZ~S%_(YU4J^n&?f z_D@HAzE%3J%*vhQW}lo|wG6;=tkwwiDSxaF`}gqisbU~D^S&UZE9#e}vupR|;8Szs z(}zFf!4)po;&1N=wU;){tgbpRP){(SF!n=F23wzkD)E0zc}nEnCVue?wPJ^cJRtOb zm+@}tSWxLuZh0HvaNVIPjK50>bQvY_e7|h*K+6A4HsM8FyD6a5`4 z3v@+tzWVck1cUOngm}>92t03?wumK})qIyQR$xiMkuAD&X33!c^Ljspz8*$6(*KGsDM0bNOMT(ZLNxN)xtu>_f>g)47UeUJNKRVy zOcT_>K8=6QjK94w(I|JtgygEuI9QsbYwO}+bk}Ov5R2xzFkC3hncza&3My?bSZjplf)1H{elP>aqsTfMk~5k< ze5Ij|NdcW+r}>e08Wa5Gm##~Z`OA%({ZZR}$a#K-_W#@($)(cj&?HSL%3k!XGxkDL z87=4AI3>8=elf=56OR^0;Rnn7no!X^@JUA_65YpQHkZ!(L*}KT<%&WKo}GR;;Fjit z!>-l-#|sSc#;ok~Gv#nNt+wrrA>139njAkD!ZXX_*eq~4)B+n#3)gL&6=7DS;XnDp z8B^k~C|7!nuw3Jt%};ohlU=Ks!zXQtPG~JY``C;0gIjcQ!xa<1qKg{aVj;y57TV4j zfm|kjCb_+2Ud7SkG#=xF(uMMva-} zsp{CMe=1T=aw7AZM>4W1jUd5nLF;{n)R#4HmVOf-;OVqqZBaUV%6$5D79Nfzl>UU11u`@t~x z-|uT>{s=!XX*}2%i!SAjz^e}lzwii$b#}cRx}^m5@}1N$&+j~5*C3PqX9BQwixCwz3??OeLo2#ekm3L74pxD}FVP5(Lya&-OCVZ=Y3 z?Dc*oH6#Y(?Y4(`Ey*0|(N*(UKSlI0T8_Y11zXd8q*b0UfPK>avIpmkpw_eLlKeg^ z!cQEj`ERQ}9Gu^~?ubtYJ5R980b$Z-`rkPi5<`4QgO&UBRcygvt+?BC$sD2NM1M_w zx8;J5b+IShVRx&At(@pAc&HDNB%)`2J=Q+Mt^-YTdSh(?J4iCQJ&zj#m{;#rDt%N% z+OZ=Bv|0-2Nn`lcMfj-wZw70XA31=eu3a2``~_fko{wLV{4C1{a| z@)i#-jk!eZkp2Dn*{T)x<*iEjXXjzZg;BvbRr7ezZvKjYn;#@64iQl zn=_AU#4*X7sS`ozC3tgN)K&Ej_heXArDRix_Na?7i4wyWu|hfLs?9s zW=B9O)c%S43!6D$PugQAiS`&M4h7_VagN8XsMu|qthNZRajo3Tu8$Kpk3HNyt$^LT zz4(I^$^9`>oxJmx7Mv?P&S>v2z~e)$0&5wL(2i*SW4)q+=Coj^4SM1uHcxB#`Hl2> zQlY8x_lOTnNK$p%S}4X7pRi=~hhia%Aw!Es4LOr{O-dQ;i4ORDxc7|{`q~tgq)D#( zf@ShAQ$87R-1xzfq-}(h^ZVNRTvQhza-8@z!a!IHNCY_}nfYnD>||knqJW-6Nme<{j{~&+fVTVST)Q{!Mf6 zj52D^PhL9g>j7FRMh!<-eLUUX8z$-+h)=8S_mpgOFsGTa?m&1IdDMGGEgT|n$VUB7 zL{TaxSSnbWZihq1x`Dol(uMHgmCAXD-ox-^)Txlvg_T#@I$8;LP+_#qQJ>6t1DL*D zCf*Vp=}S#@T2h482Y*FL7g?;+X)FhKg&^vA_|(v6;!D3)X@BjB0?wAvY|wceZvz=A)mtnTYTY~L^XW8MVoNJf33Y>)rUbY z)lhw2D)g)FU!~bjI7>VclJPIy5aOBD=gSpB)<2b47lDPt++OEsO zedar*txzlk*I${$lRT9~M$V^tW=&LGNdGPWjpT(V6{xwENlq~9!GAYMzp#1RT_*+7 z$5##TwF$+$K#{}nKC^`{wpDN8U*4*XX94t<-=z&O#9eSo#zhV?u3ibkI)Uh~Ub{p^ z>Z8ax$zx-$tgwaMZ|l`0!Udx%K3&@ygACgN@f3O;a2;#s+i9W%WwzkRhC>cW&ls4Y zjSt4nhwH^DACmC4dUU(T9Wu`oA5k&bY7S91+a>qxmc5pjnO+ zVj>(;_AH+-6}lv!*mWZd150^xWjEEB4w&^LLTC(9AAob@Q+Qw5zV3nvf@aMwMcT=RS$>ZlCd{^zp-A z=K3QJwYkX233xXa=LTcizlSaksbid~xI^%WF2vfzH=X$BjAugr&nfr>*2q(d1ZtvTj6(&&llxvbVE05KE_7zp6 z&h6{Dx=ZP@9oFwYcDdb5ctcs#ca1bnz>}IBsoa`~d@nI$xuy`r-8*!4`I8qy2PK$> zg8lHimD20Tc`F=cez_z%OZ08NsW+Ki#BVyh!`S$$98!#Xw#{5k#kLE9a?w$GSlN?Z zPhFpd!;u*QB_=Vju;e-VFUJ;4EKNsh#0??nYPIX_EUBAE{PRR=Np4&EzSq;^ItVgU zzx4UI8?IaHsqItnf}N*}W2=od26E`?Ht)B`jfruVtsR6D_3^pNnyNjLJDpAQ6+M(D?X1R#i!eSuma02-+_OEfL6(9^GZB1NT!h~<9@uNBGr zm-B)5ORW*6YF)fEbqHrt`Ri4K^PVtowdLy$RKdnq@gq77srWwlnaXiB24_z0ZOGQL z!N0D%Z=OvO9c`X7W8u6qLhR0tN2_|{-m4W`I}>xfZ1ob2?9{<(gXF$*dOjFB{^&jz z%~`anZ_2!|FB1Ho6^-nD!H|5btu+}=bYqP|YG%SGnfccFaor^oeD@Qhr#iLJCZwXs zNzPNtRdyd*OBE#UZ&3>(_xf0v-tlH{MHoE!-61;@j$;4!tDGl%&=Ita;m3>{c)4zW zn_)4<>Fg3^`?x4{irD9V=O*tj?ekny8&~Ar4$w@_3?x1zia~F}0o>i<7+)A0js-{V z(`L2G$ZDAKD~piFUAIGKPdBu1EH`+6l~*k65*)H=3AdxHt0HKpYygh=J8CICB>b3d zl`nz`ck*pnqXtE!7d9UBUjA!L^dV8E#G40yAfV=m{}lDb?Wj^2*a6SGVB@`&sq z*;HP9?h)a*57Hj;_iDC4+`cbgKhdC!do-vAWT&TWjY`9wsEUt-vPcJBA&yT&cX`Fsg)@XjX zJ4GGAR`Zv1GVL*u@G*O-+!FOqR57ws2?r@yBkZ5q!Ot(vS|XLq4->HS_e41CR<6}k z5pG1EKJ5|vM!plc^!hhl^^v>M>JvmQFZ+&;h?)(rn~sG1wx_j zgwvJd#IO8%ruv+a2AF7tr|yzm{qG{ym$&$&q4iioR2_5awm5Ta zQ@B3;o78`++rO|Y7}z06YTS94UIUiD6}1w5^dOS9*roAC1yf6{jY}UL5u&UiTU4M6 zoBfLkcXid#qWRQ4ZaNaPF}KyrxxLXUNkQrC=71hiX2o#QZ(KGNc&IlBl;<9|cFqXH zeGOhazu_?O`mpTuB6@sc$0lZ1l_d1;$a-jN5rG=x2{pBIHMm6ko85ur;^~^yovs;# z0a^O{KoS59K@s15g>!rh)B8htuoMaxObqE!Uiwij8Wc} zk|91kzDuIfB)^ib@9?r#5BhgPs&nVg;k|6(oZ?qoL^!k0dLLCq*wc~tKYB6v?}4uL zEvayP73COM3ysF?+h^OEl9CXrA{|Pbn~36}8*cG8P2sk6zoYIb5r7kDzgckmL+J-c zk+-r9!rpt*wy6?d$+z4;ZvLToZPfaQ;j$&9XGNznuNss65Y{d$MJ9v?$Ps_k&;gW` zuNDdx2(MC<@~?!pIvTI!-(aW-M(F<06{-p!ROfXH(0ZDqASKk8V!)E{iPL)+-;zH2 z&S?=UZ*#oz>;V z80ygexaq7CBLAM7Rl2SOPXTlPSCyVn;BGFPQHz7-#GcnblGGvgNbA?z!!{^iIh{a#MJ0 zHDZwaRAk;wz#3^?dBQet_3`J&(3gI5qW@PE%5#vq=x)~Q$#YgI_%69r(7h6X#U4|( z%TH`@?L?k$zcb3ybv*k*lJY^AB8oOm##T$X;gII) zuUSCh(+qmygk=$BP z-iI~4Xp~Xkr{`HGIW!&-hNncwx8e7&US*fXfl9}3T?a@G!FHqdyt4w*#@_uCB7TL; z&}jj$gX&n(t6KcmZjY+{!s&;JUwZeDhrwrQK+uYRIZZPV^ivHK>zDO$yCJn{$GkmW zo!7m3o-Yj4E&+#MeGLY6$AhBLCN0XQ+;+7*ATmqEP&K5GHPW&FFmdzvQO5^dw`-PxmhYuU1qbw>8=?3*p_ ziC8hQ;@}`TMj4I?9ti^*{N+|IncZ&yp2~RpN1CeGy)Mq58fS#uZUvQWJ}+qA?_L-u zdQ?Tv$E&%w)bKAjwX)%WB2FUDioOr2^@;yZDT>7<+y?t9T_fyu6u_lKv-g2hAR5l>oZUy~fLC=} z26Si0UW|ideR5Ss;5s-zBfD+`<-PkyU5|RBdUNz`jR((+!hy#t2 zz-9W|aZsA$3D7b&N3E#R4wqKK?aiIpSn42idCFfxMnr#2k`$0XlWv8zVf)i}2q)Q| zCHJ#=xDEmsxC;FkWPOx1F&+UU3y@jBWQQNJN@5QMGQSu zO5^MfMNEbcBhThAgbWC{jm!t&y^rf=U3NF@{@qr5chLv*tqnW#YkYCaGS6a+N)H+S zxfK5qeqEG>*+b23R;1q!Z4U?vgtkOmN#Jh{9Bt9R`;WpAoe_csvCoun?QU9I%Hv4F zn-ODC86frHRlol{NZqwDDRE4AQxuGkJ1z{8x!0H6-2P8bBeB3-PxG}a7%~Rx0xKGX z2gbXTDxNQyJdb}2i-=E5Vk{arFA$#O_x`xY~0qIDL;7k{-nrhEpkHaZKeK zf4lFu{fIdJy5MK6=N)?}N) zPpL3{q<(4Ht(Lxc88D$4Fn6z;^Yrkc|xj+;@(6{AZ^26J~Qr%BVWKUhjir1S+GC!VV z7TmGb06GGuOYQ1}H*+ObJ&Mi*2es=NZhPxs!JvFwl&d20I=}by#R(&MFv-E7%K^0W zapAMX55}fGIacWC59cG*Wn#W+*fTX0vz6qKsKm}R9-_~M!F#_Kx3s)SoiDJ)N9yQT zvTsxQm!e^kkX2QG(GxVAt!8BiUt{JJ=lymg4;)DhI@v>M0=oB;QWUS^aGA{0PYW3# z!-<+H=bAF%iHiP=|3rBD?N>DPh+gt9vFoHJ$xj#1Zx{`c+?|ZReNZ#AJ>ex5ZfHN& zfy|O0-H%g>xPAH938!{HFzB&1sJ(TEh$cHfAB7VxCw%-cXJU`kflpZiiU#Pqy655P zlWH)nq;IUb=K^2Vr~@MB&qDXqoH6g357+}UU8N#4aN$AMIiqDABxZIj&XNA4bbQ~& z_X$!@&ZeHGDv`sctL3|TMBQP)_PH)liRef+wzZl%WWV0k6Hl@j6LE$k;~(3fNZf83 zGp

d*^sr5Z7sg&i{NB1M@6V;46IfG1)Jqu<2P~Qn?Z~O2lJ}2Z^5O5#LoPst)SU zT1N)%ket<$(2Xg=-TXmmTct=Xi^)BHsii~S;6C}E`a4?US3h)T_LMLAy#-3>`pJxy zQe~A_XPXg})O%9fj+1-!wo=EGg915!RU?&_WO2vxBFf&J#>s=3{~x z*sd0rF4k*~aF6};G8+;hNIdKuG!$?ZicXsieOm&NAUk5mblSg^6G}U4#b0JT5=6CuqUtafg;IK=eTcH z@Q`uFT-!=*hFz#GDYVC3iL(U}q>easQK~nG)T>efl#?2pi7#z9utIfjAYT67 zr$kZXfm%zuPSG@dTy!xsQoQYig0%yE4$l4F30BgYoJrfy;Nt3n#<=N#HdP)u#fGO;Ip2&c&X&!_&j;7{%UKOb`}+#By>O~8Bday(l#(OGi2ljqISL2=Vt^>LvEEXH&9 zeZSF;GM(*>nFU0j+MTZXK-zbcB=|q1kA%- z%Sz&jzj5$T;g1`p7-zXm9oj>DF6~}@Kj&2OAnVv6w@+%QT^>87NBZXy+wx$BMJ+T4 zdsUYG=MDWXdq!i@$1y3f>Kk|DV>;O6&|i|f(mu1u!#|?}{zz7w9kOd8 z=Yf}a_a{!G*NEoN9|%at?wCMbA={M_GMXwDv|H@FigLc?%cTuoVndID_em?cm znCzn%wRVm)^aQ(CY}wFvb(A?i>T+2pT4nAuYZ~4orh5{#1_X7&wZ(a-q77O zXHLQ=>sQ@!b<`by1Uu<|NokVzX#C|6;ld3tCoo^K@PYs4%O(nmrevR`<=uCLj{?nY zF9i@z?i<#VfgzswlS$rmc^h2yW*^8QKAhqwHyP9G1z^>j`Cz>}6xwcm0(X&>RLjZ?89`@*b}iz_1H>s=|mEInnmN=25q6ytG4?bwI3pSLl+%* zbnsQKW@r8-BYdXS${e{%=0*Xl`!|#M#SSq;g|dDFeAs=PJ%Q}`Df4?q$A423`vRAn z7j}^PC;hi~jI2It+98=Cc>-T4#Lv-FngN`>D8+5 zcLRq`)37D>MG zd$Utj&_7+xq54+?8=tFp>@FpJwa(3-kNb$J!XKMf3rgcL|#`k-{El6Jojc&G6H83kH!g)=|R~wIavE0nKOycLg|;CporRj{3x=K2t>fiK#7;u;sj;4x6pAQY+zHPyJe_ZNMklsj?pI72XYZO|^= z2qpQeQOX4gl@x4p9^}Y>X^-Z}1I(H%`Z%Yui?hNa8!vM|-;F%ZPk^2+~i# z(_m3SZr4f%PI#i{+SrTM`&tkR8c?yOw8G|2mS-|9S%A$GQC^A!jEs6xh)5{Ha%y{p zm5eX6h8o}R_X#7s%$kd5iSLr;Yk`j{@!t-27)_Zu`H^`h3WJVX!TiS;8lgp1Tsh^X zlzo`&3AxO`Z&9KI4&6KINusKRf0y~^Gnr58@!Y);OXhRB+?zBMZ+gI$y1$v@lNIFq z&WV~HA^G$hbL}0`j%bTBaA+j;UTTSFFlTEF7-WS@Z8BZq>)p!!$zddn(8Uy2*E#k+O^uV{`{^~`88z%26>-cm&XWBf^54C?J3gISHO*_}|nbc$R!&O&4nWUiRsiX*BQ9jZwX9aFo z=z)S$fMIM0*<(1Tuaj-|e}8-XM~6G=`1Wi|L>TdBgr%pH`n(`>>Mw83vL6eA`q8s~ zRliiQCC~NNQD$Ep*qFcgkmw~J#v+B=XjC!(N$RF`qA#MgJc|;wTv0iA@an9oE6D>! zbj<%VM4S7k^wvZ(l*(-^+1w)B1=T0)`xGou6?8fQL>=a$&_7oAlkcaPA=7gQvB*u^Ng9 zo2R{-)nEb^1~%Rhs`Qu zuK_~{hbG~TvkNp>*vo;#NpQCU$vw-?E`HD0FNg0>&Sx)^eb*jFmuc?x8$x!+vBF+w zEtI+Y>(hjiJ=$k(WxdEJ`}WRmEw&*(+b?cnUUkS zd?xznv8?=ge3u&D?JHy4bT1UYsb{*HNuL=?uNsp#6OEVb>XfQmNv>8ZENGVG)HnpW z2Fk1TvB*XS6!iKq?OzZ3?4|(mp#0)0Co&hXr`{<_aw++&!5;U%$v~?psLAB48W=yE zxcNFO9Wuu^o`#S*l9tBehutGxBo%oL+jXVl+X$2G^b0SYZ&w3+ z&Nw9$k?(yqPOQp>=u?qnPhauJqW)|x`#!>N*qX)_vz6$Gdb2h?RWHL~a)iVF+F?UH z5jv5Rmz$3%3*&@wl7k$zo=G+9BwX?0)Bnb9=|b>5qYoWhwjJ7!$XJV5pdJvo}H%#B(mG?r4xVlymrghiZl=8TXV0RB)*mqSFKN4 znS{&yT)JxNL_QWw#Vl7`t)MxQOUbk0j>Fp2G%pEH_Z+qCyLt*6EJd6@dQLTIFK9^lcO6OB$4I@wR;dkmcG zFW*A=XWAu#YE2|R?!7=ce^l2DGM&t>_x`!yQTIX?^Jg;m*cH(oPoII(jdVI=MZ&L8 z{=`$kuK^|HUyQoJ2?%0h>l&_cg!{-kCmJs#w>cnvh7Huw|OGI3BQ?AbTr4q zwWIXk=?(E<*Im!!pY(CvtZRz5SPiNUu61ggV-cPd%hEUQiEF+bUG_7PWFI>6vWf4a zSX}tl$TkI(QSKBQ3^9fbC+jx%cxOx)YNQ?<(!&(HfaA-ns_x%;)98|7;g9 zCv`EymHjq>m0mb1F{!Ry*K``g)nfO~u_@ekayy!s?6p}fbqNh}# zlV--RSnZFlZTw;LgsXEmdh3CE8D!qw!dPqc&6DI%jf+QC$lNnRn_&~ldl|ZqDZQtc zgYW)-2EMH=#83NzCNL)z&sCzT2Pd)+{%h}#*uqfEf9Mz_9(T~ZaWt3IQo|PSgRkwn zNUrO!lcw~%0p`L@ud-y5dv2|lD?mR554Bn&V~li>!dz{Va#J7f&bKm`2v>bws5x?d zj`UYupM?GCLQtko8HtTR?}3=m$7eEGrd zc#^}oQY6I0x6Khuzp96)naO>ndwDbOO|p0DYr*@aF5(j~%4>?FPy|Q$9-He@WY5Q+ zn404^$bMnkHOiY;LLmNSot>ZAfb880_!mub@*hPHzUl4t!J8t+8P!&8_$l{osv+k} z3H8>%%hV)aaj`v!_KqKf#1zUtMcbm5oByz2nLK0-#->IM;;`aU$Yl9W2glbvc3TST z!?DKv=b6$lux>wZw*8HMk7rjB=abfn{e8PH|?8t|Cfv3(zF`}ieIwh+7_*xKVjge?U`8@7 zS2EQTg+JfWJN(v$^~K?v;l~4UzG^`DCS^YQxe{^@N1MVt`TDF^e-u9Mkz065_-L&! ze+j&s)zNUpE3q)tL z+3K3p7ZZr|o)?FFJk{~%in;IJ?0CHBm&w>2N&NH7Cp#Lxkp1RrG0CE38hFlopdn~g z4(?43>Y=1A5Xh2Xxt8RK$8B9FDa{GjJLnl(F!7gicRMEjxST}x-DOw|=DDG?YEe9W zl;q2j3R@14{IT~+r04^ZTQ!RpxVD4ng}wA9ybpb~F?Q0)JoXaVUuFEEKZ?}p-?GAN^1ejOXoL(n+D-Fe`HJx7wpWt8Z zS1lW|hYC}QxdXOvWsaLunv=tGCfA}qN_Pm)jVb&|1|(29~&uKWO`!;gG$j z0873uu75IjWl49!a&|&Xa8M);zL!??Be~CPT?=i-(^`0adH%yA3ME`pbI_hA=fy1} z&Vr|U5m@`@RWO(s0M;{9=YLEPKG4Wr=C>JVaMapx%PBA7Bg{KRv4xBHBtnO%Xct|; z%bzeNRqcTtpZNL*TKsUumNzi_s{__9`O;qsc7^wZ=!*4U(szh-HqTrNM4AMbTFE^V z{2Xczj!;R!u^naiPmg&)GQr@=$OU^iM!MGS@N@^)jMyTZrzNgE6nyTWpoP1NJl;3D zqfx#5t*x7j>|f+LHPtH{jIq>5;w(4G9?to$R>k}{9G_?YLr?s<>T9mOy78(I_>b@H zi6k3XFiOZoh5L{kKMstml6r$8>7$NLCfJUJKjb33?Dti!7}F}pyP8Y;{-K(=0j^Yl9tl+<4L7~@eQdyI6iACe&XtwlGcT@O`oAUuCj_%q3a-iv2__e~2| z&Sf8YCU1o9%=B-YlAQ3U)v%KPt0^w+4SXagPB`JT=c{h>JHyWZx}tS^5a=CjYmH-x zo+`r55kjRX4Jo5Elu$}hRH!sGr4kkI_5S^QdeYDRyYKIPUFUfmn21i7 zez`6M*BgJoTdWs`Zq8L!ZgjpKQH<2nme9eXAL(PdtJL5zDST?*m^GFZo4C1hQ|`xg z{SyHZD@c{Bl|3D)jop9GUT)hM0O7wcSvxsWkhkMiuk>*Te9YS^u~c3QWrwtMUQ*5^ z_mT2t{u4UHgJu{xr(=p^KWomNwV?au54Sgx^nFrpDce)3O5W;MOwt0Qo&a?`DmJ8{ zgY(@v6`@S*A5?Q!-=Gh!sefFxG{5>`C9s>9d>}GehdcaF$AKkjB(~&FEQZW2LRYY9 z;j+;{>TG}w+6o@Y2qeVdW7|ybkvVN-WGygjQ-AfbyN6Tpd4KFwEi(_K&xKJb=iEtS z$}vlapQ)hp-zSD`MXc4T!)WHp;u3g3BiFcu`oR&w4FSJZmN`dAC@)5krAzp^~VlBJw;0Kzo=^j;` zVx)X$Pr?L5jPSbWuKGtw4?NBP@9e7U5j0QQ>bZYI7<8o>3^R#qB&t=lQqMXMN8e^$ zd%WEVc6an8Jy~?n6=}4Rk#Z+9;!9tI+;hRgS-%3Q_sJM}?%QcjJpM_1UB9%E<``;T zTjd=>VVBA4`hR}D(f7K2=QZL`>-Oit-OD<3j_^FR=xzqSzV~wa;1r54H?T45hAxDJ zc5SU#qe{6K`#;rpETO$3cwxxh0_v%M)F%hDp!i;`DK)|v951?aqB*rtdg0Qh#TRs- zmG4vSIjW1SE$>U!N<+|WJkFhMM$eU1_tniqBk*YDT;3L04-8+g5_+-767k*~x1+w0 z2X8`cc~+w`r1Qn2A8ty35u;wf<-8!2ubYV3EANPnlYhH|sb64jdSit%gB}<+>K_fR zq1)>y2WcvFx*x0lSO9x73Wh<0`8laq_a-sNJC59OPoKWndDU5+x?dR3)4 z^J#zc+PtNC${dH6e-NI}*MXRoN!%#y0hCM~vYs+gubbO6I_op}()MY$EDT zAnn0DwC)WRSdcI7;QbX*E*MHZ!)5MkLtcUjb#BW*+#Z_wH|`pSQCav7s*=AnQRKlR zStrQ1Zem(Vyv|!6H?sy*4uoh|z^I2`a*?f_?Oks6Y3h8#3vj6vbE4KFMzY5w$laP}qT*VlGUFA*;a zhU|&27Sb8SdAy-qV?;i{g37P8gEZIIXA(WQA;c6Hjz0eH5HrnD3#V4ZzM=E-=ZS&6 zb?VT%C9%q?B^>L|x2$KJq#O{VtWC*~8=4#97wZvsD`DvEn16~A_ObfN9By)jWoL7} zU8@$>9g=(b#MueGE0-w7c)&Y>?_lzFJ{vjsF9S(jl?TO|h^wGK_i+A*Gu(>x^FDtz zz;X7a3QJ2&z|wR0xpR0TPO6?9TB>V=+lLo!)w>hlu$Zy)WquT{4>F_%KhwgF=2tB| z*5Np!dN$$=mmVM_krhw&DEoRK0pCh*JiR@mkBfQBJj@bYu=&4LlCLupVd)>( z^>KeHG!LvBcRgf?4f)1m{UKg>um5tM`(-TF^+}9-WdvhMW~WrXQ4kca9`wygu13kj z_Oj>a4j^jh_Y=-izF^wBl~E>~=FW!9i}u+Tpy9r$l|*g~4Bpusts?J7#>;*M_c&F& zVDJf5k5#~mhYoF3tWn^wzPoDsfEzd-C%kO$rt?~4$Faz5l&AAA@4``gTnfpk87iZk zufO0mcH3YmV&UfVDox7UWW8(lOu$l=+FEzYvrHRC*2}gKhxn9(OgrVzeDxl0@p%x6 z9m{VXQtj78-+uS+i5#x@l`!@(kI5fwJ54-0SEYh^ufWj#({>2BaUkyOK@Die&W@G^ z)A=N2N1D?@8J@POmTwkK$G_k^b=Hxzf8b;P@bp78T9prXs|cClt*^O|UzpDS=Wiwo z&bsJ`D&i8!3q;>aYb`6fkA<9W7hzehf@3`~$3xl4qnMKU<;rLtMjIR72iub8?5@9a z$28^qJMxnS6a&GF!cxr z7%KGR>F>+A$m^hqRK<;1I?d(?fA9Qshl2)QU)n#m?3XEG41P;L=k&veq9-z+E*rqz z+`c?N)fjTu%N7}1o8nE-9lM`T%n_T)&}HY9MR_6hH5MICJx~9oo0o_D4DYgx9(Ud_g`7hB4qe%NvkX>?r=c+S+_-1Lb33? zeQx%Ad>oQFcCJ5sEgJjJOK8QL(tBb3xS{5A4P=G}FqMva;G*n>MNPDy$$96wo11bw z%EF8_*Up6Dr_sJH*CifETC%@(?1K++s(bJKb5X~w=K*EGM~G*I?H=dZBXNOUK4e9P z6K;&nUY>sM3BIjia|wr3vCtE)(*nw~T3tKt8&5pvTT(B+Q@_&iOQ?7Mup)jnu&UHh zuP$FQ*Ivbmc)u+!55LHz;$-C4s)Eo&B)u&$oz-f{2$yR3@zNP3UR}1O zt{UK!_~f*fI1ML_7O8GxOhkXqp9HJ@DfrOQ!mx?X>Bd&V1B3K@8**9caM{8Q)}_t* z1%IqzHn!Y2a3q%Y@*QC!)#8KcPkimBnobSUwjAF>xR#@Gl~ zHp5w89R8l}dw9kN?>}4hy3Be&l*i=GZt9b~XFaFIZ%)0a<#Ef*gX}O<7V$&my%VCG zmP$G9(?EPe$J;Au$_V&+#CN-OGK%J}MH_P$V%oYO@cFV(ILzmKJ0+ZeiQZ3b8s_Hs z%KLaLtEvOWb=(BU{j9*HxjU*{(-iZ;KiwOaJ3{N)Ye~+PPEb9!BzGlu9In4*eOAz= zfw$Qj><^4hsW+m)^={k`HZ2i(l37NWJ~NpAV7)rF)Q-8ge55?5y{hz|Oee~Z-iQ;Y z`=fns{qlRVuJ|uBT(rO|3K^$oy$=l|a~%iBs=kP@4iGQ6*OHrIIKd7okAEyLjLIa= zUu9qgR}7X&R8^i3OToRrOy-ii15w+2lzVt42@C8N$={-#v0_i;>-Brdo1D4thvi9C zIH@SUKYv#VK9wpjJ-#NQyrfXre!>VJYX{aan&7 zy3}FXzDtJ5&j^3RpH~jQRG@yuziiaFqUT0s-aRJb2L7>f8*fxa^T#foTVmP>3Y+{f zo*js9{;TRLJ&d7dwseOtaZ_0Hk~3O2nnPyw`nv_;`gpzLPw57>_i*97Mdc_bTTu@ws?FIe6>K%dd=Z>PY1L<|WR!sJl$pkHw38 z!}m6Ly*B`V>%_dnI?C(m@*dt1V1ixoB7B@oZisuav-Iy*7f9~5%&Cmiz>J;|ONF;T zxat*NUJ`de+qjR^%^MDgSxCFv8W;@0!#}S5w>lXP%lp=Cqx%KFSs7PFvntr4cj+J4 z?S?}61V^EaJ=&Z- z%+!Yru+}xTHaAcSe8Qp83F{T%w&zXFjUP!6(t7f}ATa~FddvQFaF7qcz~zrtcLaGL zhbQ&e)38%GqFC410tGwl6Rmzl!tuGJ$hs$Vzv?b=-l~)au7=*Kt-j%KWxk##?WB!? zBhqsT&%-cd+J3s8yj_Px%~np*{dUzUmvT>jb^Ll#{qH6R^-3bMa*E}s2YFD?e~-Ht zV%B^pJmjB@fQq`e=95Ch->-oe`nLKx5m17;08;dw~zl`B*9)br=tU5D&@on{$ zktaMU5T2+gotlWi{F^6n-}v?5Y9{-i!j z#^3!`SygSWc%#~4*yKY!#H9i}gG$!uMAhSu^n2*Mb&h$Q)fu*J8GWe+VPI=JbJ1~; z7Nj>a+ zs}^X(`IXYOqO2HvxGh%pw>Y1k(0uu(3N_y$ciJmTB6k)u_n|CV~KtugvsA2X9+I_y}Nt(}Qe z5?fYG-XK3r=30ih{2(0t<@K83loK{RUYo_vnux?h3*QdBcER!`R&ox+LGgLF+q3_v zHnx>1ob)}ffWWf7c}@xbSZZ^aLqRfW9LHYm+b=G*k-9>rLr3EsxCGqGfEAzfR=7sF5AKPY2=)9)(``F;(JScB# z5uWfNPt&)}se7o0njK-)#X|d#qlv@ewZtvi^6~B6McM9XDTu3Xex?qwfzko{VdBvG zSaZpIDZr`IdXawD0wHpPk89t*1nNIHE_U5aydnOQ(`tIOPg+qSyNcfj8jV|j)hlRY zeCYF77Udv>!^{c|*V^Lry}O!KTv~W9RbcYYj`sfdt!^Vg2aJ|$FG+WY(>zJBdjs(k zx$N9h{|p;J?n~_|bK;6Hwl>(m6iq?>{w7=HT0>mgqPbKh*O{(|V+NB6)?nl`=vyOi zjw->|F^)q4h+L8bed=eJSBT$Cj+H}X*go_6gn0avy0P)Pm?h4Q2ONs1Hh@576Pwr? zL-3s~oB4f$?)&e?Uc8tH!Dn`_&A0d*;H+V_rLZ&rO!>|&+n4I$@u_3`5;c?YX&;-@ zL*j`SnPezz=CDVCZLhU%uR8d-MW))RA1A5DxaE$zF)B7E-|%u&!M(v6|3{kio;}hY z%G=`r+JW-(XS-Dwybe(oK{Og@+>&p0995S(80&-+-TE-o!C z^%)j-!T0DGA)dXeh}@@{DeaSta(2}cWtuN+ipb~T(~Lys3$x$cm&C+r!s9qcUQhYvvxb$rgifOMRWf-lyZs>yz(LCXBK@kJRg&R_S9QUvd#twu9^E}kmij$|7~vc z4WsKs>A}i3k5ZAad4T7XP$D#5G;aO8GzW(~4l86(ZfGY*@PCn0#0k(F_UtB3k@C^T z7p%(@AYLZpH2ymt%zSya?iHS3`gS$@APaeFyCe@<9#cWq^uE*f&gS?vmocg+>V!2* zB<%!4X`X5*(ZueogRtL~#?RjwVdbW`yf2isajJUG=h$N}lscuq^<&dTJ8QYT(XM=0 z%+<;aikCYo3Oa6{dX(zLkszVZZMe6&@RN1hpRfrCp&3;emTg z`VMJBS$Vq3r&|+}R|FQHp9_FWLB-)+Q#$xPxY~bFuQ|pa?ur(&k45E1w_iD9wkgu1P#LSZ4~$(|jofpU5Lu z5|AGvP>PqP4;yo{LlCFetGqdr`o2XEgr6%1;CrRV#Fj$p!>v|(w|gD=>AJ3bb@=TF z!y`lQpDLwdb$<0=*^g*E)AU|@_*FQ{{eLaqHb_2l1LsJWU`_06`d)fVB>@fc6E*wh zh!YpCZJBt;0bA}0Cr_+O#~ZtF&msE^+Mz5CO$_!mR(D|F!t^ZtICMSbcBI( zt?2du{I=_~UqRfZcMEq-O_eC`{%`Z#tp#Uz-PjzVN8VCfv*fjh{r%8fR-QfN_OsU-ERxdkBOqyp#THV#sQ{zN${WECcp5F-ZPeJY|=`hJfASX zPOP&Q`K8l+i)w>0uJvqSh_Mi(uRL$}wYg(O(SoD#Sy#B6RQDOO)WILwRa$$iC?7ZY z$*y=;Fxa>EysrEfg2=O$8+IiaqFcz}Sy8+utUq>z^;g(pyR2@F+Gcfx1qMD{S-cxw z*GqVd;?yy&p_;-XX#&nH)&|EdQ81A_pD|}`4=tVnPKo*yh#l)Z<48HB*Ba$Adc_{7 z=Hp_L`WOI|g)MUzOT%h`rLohA)~KmuSzW!;1-5LjofCg)!s6ESx%3(j1hB5C~Z7eFx$rdCJv}*XJ!;du4lE5a!HcVXAnFlI$ zO!ViSv_h6vS>XcDJg2j|I@qNLcGWu3brdb3s9ULyr#fZ7jQ%L&&PW;O+niW<>iv?{q6B1A^=#x0x()>BQ)M&R=QhyPHubyv!NO zbuqp2lrQL=df6r!O&8+3iw-4&4oEF-94WI%!Xx%w-(&}za5n0`L+rH#Sb7TFE|*V* zo{rO@HQpvjf4;ajdYvH>oIN(ruL8F9xoJIC)xFkXo#EDQ$dw+;6hkTl$2~3<3 zk@Wf7$q$wpU|_&C{4GNmmR#1{^vb~qHhixy{v7rJ%kdM7`=+A7eqBDAH7o-AciJ6{ z-{y)|{ne7W%H-Ry)IV+6NOKs=>nHeUfWR-I6-7+sN!@R?{{2B8;Pkvrz_}op7#`3_ zzi5X4rc!tgbo!vXEk*k35d(Z@%u*TWP{Zf7eIjj6D(F{g4Klyrjb!$Vg&llas6S$J zQ}2;J9{&36r^lO)Edhy7(&ZDuFP<=(>Pfk&LuPyrj2)59w{z)K9`zygxL<}Rlegpi z_X{=i5m2|w8RBKtMT_O%Ph#d?@XF#UlNPqX6Yu}NeY)p}r502Fal~rTb;&|&$;t;1@R<@W3!uxoRjSVRV zc)Lb>aDLJQZYWC-&jwr9WC>@e25wg0ZX4$2%#EX{O05V>gJ-JS+N2;2$!@a$qFB*un@ z0_pqm<@3*_O?witiou|#p~e#`5otcz6zv3XHfHBfV%s7@GP={KJ0nhB$a9FE{ zllbfeyiSTbt-xJDm5EiT}|=8 z)I?Cb+X2V8yZF|6t3fTPYKN>6y~hiDajL5(9=LPUZ=<6ou!+oiCVj#f%~4f})!AC8 zuzKYse(?xge>InK2|9zJQ!?`*?VB7NR95vqHpGot`QkgorP_adLOqT=bHRC`o`0t)~HLl=J1bk@0bnc!cC%QzpUGufKSYx)ua0Y@W|$&kHjB!gf8ej@}XYMI@2S@vy!gFEk3h1 zu3HsvqnB)bGLi@>`4Gt};b44i(Gw=LQ+OA}WDCMVkn?7m}{f)=L)Z@>K3?Ywmc!!y+o(l1@emB^Ky1-E2gW)pd z#j~vrM%6{>u+Hm`_>+DJTTj2uS#RnM%d-xH3>yr{+fvMWncgR(?yEiI_?C+r?W2q> zt=R~hy4RTVz!cKPbLU2^^e6{)``^8nL1-6md{_0{g!uTmrPb`Jlee29j`xy#Yv3EA)#w{qxduWTqD0UjGr-s2y>yh!b&4t zscq9Mk)XU^|LjJ|yTmh63~cdHB2=m&cU426=mZln5;fD$?jY-Rk9 z@57@<+sBrj*TwI`hI0{Tl;EmrUdVjk1W#`li8iZhK*3w>bW4FEB$yP^1$XITezu`{ zN+loXW3Dh2RoFw{X2rJg2gL{w$k(+ZUeB`5$$OHI9B}wVO7g1@?)bUV;J?M}x>#p( z#&G(eH}0-pzu_;tEqLx9UVrD35xy;#?4CI4374lFH!k&2zmEOKme6uf{Cl^?vHuqF z@2)C;7NhsxG`?C_mOumY6P3I5Uyj6@H$#s&{N0eX_R2{q6(#Bu-u@kw50A=^%`DF1bUg=JNi8k#tJPePWMc~FlpJ4 z$D4&`YvRvV*X7Xpc^iW{_0BGdp|NMJ9+(Ski_T2D;b{Gj>u;#SJaMnI~0P@GA z#ZRs=!qpp;8&+3^QEs5ZsnP;mi{I=&D?~ljvufTerxUQdd%h`$_T@a6rkq{c6S3H5 z$9%_(8+rzRPevApp?xXS%QcA}sA;QZy;@*EdskTPXkBUrEE`4{Z0huOiq3f2kq{1%kdTHj&@#ghN7-G3HB zz1t(_+TZFzZbf`<zUs(<(X0D4aKrPcSQhmfJ46GyglotX z9MNan6Rb~jL*6g>i_GY`w5gTXgt+0oTUVZWX^zo5f35w5%+VSx%I$T|46Me>c{?e$ zyryBzDF@n5SVc$r#nT>TN?Y=1aA_zKKOTr|&v1sns@}#v@{B0_Q(XP-NdhcBbe4W2 zKb=pykEl>yEH)_*sBC>61NX=l%Q*7=$T_Z7(CCT6G4I5#xz;%7X0E(`io+cxQG(YW z23S))~bcX5Kmw6 zcolg}*3MnSmtb?0KUj3&Ce78emhbqmX`LdZS<9!L+x0O1aoc#$d*U1o)^_cjQAhu9 z{<+m_HNj|ZzjyH+4VvGwEDrWjL;v{_N4q#BXhldIv83lU^E0;Y|HO=NEdD~{dU+EF z-??9ux5o)ruOB|UZ$=R(SKR4;z?A}Vw`@l9HEwt?`j2_?u@P#EMVxG850Ka4?tYGJ zJH(o`)c%~*#A&gDhS-rH;uj1iiLCa9T$gSB8^!?KNwK*WW2^@@P2;{(bp1S4mYliq z(I4e3-tNPHTyWqh&)Rk5m6kPv@OgJSKXI8J)zBi&U`V2kINgWqQzV9!Zfd&`jitr4|8|v#V&x2n_JVtmm;re9D zD`PGnniEG_dqt$ z9)!lf3_#B={a*nORCCr|VrhdT+Y2}uHan*wGbTFUA(P& z74a_JSI;O(g@WB<^U1Rk2jM!lD7jhH22UL|`QBV}#6!0m{tcC0_|~0$^U6C#@*(d$ zrgX#*&xSm=?>#_sjf1BGhiZV?`%_cX?orsWWJ6f%-xO%AG1%i!YXrN`Is?~4G$3`# zJ!8*uWw`B1JhpH;2{xyhLT?j~-Ep3yW%6wn%$K%o`mx>!+or2lmYt@Y_r421y}yuW zb9uyDljGFGVH9!RRh0|{{?ThMYEAHDRZsTsCqZC#9JOiG)u3F?aQ^Z2K5*>Tf_@wA z`(HP)?{+_es~I<6gX&hxd(Ia4W?gXTKK17}4h`pOQtl(g z#kwL)O$(wdwXfJ&v@q^+zvxDl8up3EMQy9g#F=BQ|351-oPlrMbKA|wxGQI3#N&{NOR6GCsVCYuNiGHPnPVxs<99 zf4lsszxNh>ytVY?-a29g!{Z-&CE2_o#k?_RLRuB){tNse-sTMlsmC2nrS|Y$6=S_N zlk)9r)#5Ev$&0vT{q?dtzKx6?l;6IAIpfTj88*o z{Le+F0*Rv~qjCSyO9y=ayGNlVITCX57h{6R6EgR!VUXdJ8GdAUtm0ek4*}&Ofw?qQ zsKkw5Ez`?{_4z-oP1`+is;a=&mwHQ?UnN|xr$#|0AU-<5&Jy<(EpI)_(?RUq-^H5C zDfhvnp7rgPCe$Z)M{&3Y;Q4W}3Wvcwh!yl0e7xa;J`JT_4SGM`n^t?gu~QYJS+-}t zDf%Mzx9&bM#Vlm3xxkW6{0DZ4@f@Z_bl(_@IpwvT{F-@Z_o_^K!jbh%V~e&rbjw!c zcJio#PnE+!Khyx{mk(#Qr)pzt#^-L}DHqy5$KSh6{g8vlc`VOsC8PfAy{{7!zEEGC zb8n&22*=GN?fM5OU%J3$Qqe=4SFOR!BztRA+dAtC|5O9h14HJWsV>m?FM3hxFg-sq zUuj!PYh!~0N1k1YE0S7I>a$jBVQ5^WFXX5Xir>m8zNbC_C!bXA_cdl1F4heb9#zDM z*w#7^;sj2!aK=3kh`ai~H>e7m;iYTrYBu7|u^$p({^gs9=A@`DUSCZ>DfFjMk~Ws_vunx} z25x(%pXpvszC8a{-5T=3o!)a#+Px_feCJJ$sCb5Avr_ByihY{&J?Xm7nXd<*6Fy<% z#GC%{Ud=A&67g4(Z-nU75Vy$h-r%GreGYwT58WTC;-7@c{EV1B<|ek9=tt;bM62#c z^Fb-BpLnkN@PZ!rzvQuLT(UyKH_7#-XPodcY|FQUl=J8~&wEqtggTxG#q?HCPx|=g z2tg(46CMA0DBGY#x^Hj2s)2*q}`YgYEL`e>=s?=eiK_i_7GhO<#-;EgQPv@0{i zFQu0D$ztN7Dh>U1C<`WUeZlT4lhgyL7`<@jD)rDsCn~-4{IR=akzq`~9;B}IIJdaR zV;leIc`VMb6%I zUiq$4Vz_CWHU#A2N*9RF*n1=S-c4&4{CQY+G||li2{qj<>h$m1UHLC~NJtCWTf@2* zS(>5udAUMcT_D1vBbNL}{HJrl^08ry$uE)0mVrW7ToI{xf8;^}>O6Ouo;l})QtiQy z{u9IjNi#q3TZ`s)V?kE4Cjz1UgDJ$7a@q;&SJc}`B|yqqDwlDR&MhnBr(GhQ@H)b1 z)z;ob+$<35IPxL})8g9$AuC%)|WXOmVmo+ z%Uwb!e;3`C)_t)u01sm2dW%Ykdl2N=^ET@+<~B?J(zK?zIoqrnq>+>nnep zJnA-h(b< zr#v87qo%e*nsQHZ%LO^fqY~oIXW*kk?<*T)Om=k|!f(fwJg+_KI5`;^{V~pka%3-0 zIMd#ecSVgLV`LPBrw4yKxKR!*!)ge!0=)dGUV*$UTv9edenfK0RRm<;-;u z8EYTzJDTB%7m*M3m)01fa(m&TJ1*oKTXs88f;bGTuig^7O4sqeDD^5q>X-LaT)nfD zzNg(mJoUG!kNEDq_)+ow!Pw3d96oW1=4om+%ISN(af#b0cl4<)Cf{{N?_V4TZt<1J zSL}Dh#6>n0mVDx=H9ffWdV>1BpYP0-kvGA8;BH$oor^nrzXT{t=t8;wQHAmF2c(7#;AO!Y(SDKo-+Kl)*Oe#Z`VOQ1i5HPr z{~vGqb()7TF0_nHCX*+8ytzMGI1qjtzFjJz{7&~gUz~rbE3#B~DiuEohkgCrtFK|^ zc-(9yuR!Pc9QLc$`NYBBv}K8w+!j65ym}=f!mfi8Ihkv8&QiWfQ{<4yTI$c-nb7$1 zIvMWW&5~+=T=D&uyI>pTs2$vzuCgt4hpCW;b*U)z6d9~ynbn*j^eO4deL2e88(U9I zPZ57zHv7-h#&~2hT=<|G<_K@W(<|eTg(Ij@dRZV(F-*B_o=zXmL)r(%2qp3;W?vb5 zC-Kh!y%k*d4ey8JMfJLZV98K;Ek3%$Yqcv_)DE_-r=Ifqnsc#-Kf7b}nqSTN9y_#l z%{=?nLSEXjbvq}{s-pC&vBy^{>f65O5Z*jQ-smo?1{Qv*RjmIOwdA%olFYU#+4d;o zT;;VNdlco-J-JfCn0i2^!$txp#XK>W7n}a8Ne`h_32#ezXwUD#dnU8N5%IkKXUFJ! zSXk?HZ2Ua;x4(_GTC(z@)A^rT;e~;Hx#hTE$Nk=1$(G%N|z2PG9 zfNV?*9t^Dr!lS`R%YVAMkam@4AM>yU-!-15)7|^wo^m_mj;J4^cSb4_su2eT1amd` z($Kz&KPirKCJu|%Z*M(C93?}Aj;u&!j2QVpvZY07oai$ZH!)X;rTzBKyB$gV*nfq> zIzd=BsW5#1wmb1PUQOyc&^w6l+(@XL88qrykEz5*ValreZ=amDd9xZPQx^JjAG~?{*8}}~1zImDo8w>ETJ?fw&UofCQkzqtiq?#G0R>b3Bsei3FyCfWnyX$=0M;#&6VdEO5AB66>&3->_`y!mP zz2{#|6vn!bDzmhtVoKp*!!rX{B)(1_=6XeWv6XXo{=GHC?QY4rQH}_RZd#o1;f@!W zm_*Nu{YgXKy16w`%X092aCG*GtR<@VOx;W$OaiM;royGnI4o)U9#wkH2V6D}CA5zk zgW-?yN{{W7KVIk?^ZOZqjp}mar>C?rsn0E*n5zQTu#YAC|I5W>lapt%gblb~CUutN=;HCw zOO;b9l=I4Lxp;;@3=_?l9VP~h5Ziic>Z-FXK6IOW9OG35v-~EJf8u)3IOD*7)kpzi zE91M)yo;jVQ_NQ5JCTU$`hESchzTw%dmkmhp@iYzUw(OS3I*5jGX=e^cKE`xU8V4= zCu$w+B16rbad_}xrN&7sC^crZ0Zbk(IunumPz2i*%zu(iBecPz_TEh{_OFiiA;&;D91vx_QgVNo`dE%vp|2>!K z5sk5*9qafdgK_1x&Z->(ewgg*t{WB3fKZD6hdY;&;Gy`$^+lmIeA|wHc~tI%|8fu5 zi-#IQ`n};cEgMq=oEg0QAdjA-4c`g`ZW?1yYspBovp%*Aa=IR+zKGND_H#dyskbN0 zv{O{p41C<9*Iv_o7dvJT?4rDc(2dr%DqbwgiHD(S;$o zb|-POv#Z5f?IFg8eJJYnJ zjQyCbC`&>Vach#(r9p@5q$3{<1MN zuE<>+cK(0AeX8DTGyoF`3ik?sg`rb(SKxc5SWG8BJ=%WE6bo<7W~y}#BHQ7WU4OJ0 zjugsA`3{<6``Uv=tQ_QtkJp-eNB-!T^vRbycyf>>q*_oBOZQQWXQu}f4}x`OXZPeZ z_1wn2BZp6iLRRJd>9ICPtc|nS(XgHJu3JxKuOxo;KarRzxog3=G{)EBOk67#$!!58 zb<__oLAMRp_}XrjoxLu`x-7e*B6oxBc-5JM1=f zTHOr?jo#ZTH7A1gdcAyPYY?2P(+5A$)LBR4efC!xW0@;V27Q;E{WClxFLP8 z%&XcCDL>0%25Anv%ygL-(_Z?$C>bW(rkP>ks@N6l?-qDrke~SA2lbnKH_NQ0q4c{| zFkj;trbl|MZP-U`5Is7DtLb1$YOy}#mEIOZOZq(RlPEGRTt%)M^&SS>>=)}-lV$M z01ZxVw(QxC|Mw&uE*m30Ymw-yVUcL;^Ah@g{W0+*ckE^1qxp#ZnB3uxMl14;Ii2>C zSA^;v?gKkDZP6D!V0yNYIG^gfG*@;IcUf*wmj8e!lm>3qGm$r&C0+kkWyKN1iFLkO z*k=v)HRYBkPXwT&yl~I+a|vkQ7rY>})d(k+xPQ+*eh~E+bLw}s5(ny-nC-itX5iS? zt(-xeZ|fE>ZSN#4c=`BD{{P?cE^q97izwQ;S?CBegJFU z>V;PFKW}@+(O?!H1^3UzS=8E1FdPm;fl;{c;+$v% zuiSQ}g7znW`~N;*tk%HDaOR=C-G12dvT9}eSTu$Y$-Gc|odAZ*oIGVW?ZDw0oUQXI z7;GH;ejDn%(J65Kt(dwM!cQ*$?MvTZ^?(1Jy0p$7#q1JX=T^jFWOwTs-f3f;mT0$> zC^Uk?D9dBh+#q}@DR!zdkHfd}o{-zujd7VhbmyX+V3dly(cfi}g2IVbrCsa>xNNC1 zC`EfV%RNDB4h6bmlygLwY`_M~XONnZGV_YYGmYM@>WU< zk9muf@XCgGF1!tqTca>-loEe9MFSsQ%2+(&lQIA6vFVwq2t;@#+%ur}yiL=~6}yAg z;cc#C>wPMPVEolwxyQ5+#kF=K@H*vPwo1JaeW{4YGqa{2HEF)zGtYea5#|05Ft+s2 zeLJgN{&+jx@Bh2IZ016i3&Oe|u93PzoQH^)GEuD5a}o-h`P%A7&rQir261OR5RkBX z`_u((j+Wn_G&(|GxOyS@h%<&0i~b811*yQ!%KW#Z9ywY)6)d5Gp6z5Aarn4uz#E7-dE%qpT1~g(OPf{r&$o-u!;gbKlo>p2y)6-QPjx3+MN8GdpmYpqtiFx&Y%ZBxgIji92g@dl2CT zGQRU8`*-I(u8+>`aDvQv3)QAjd(2WFlz6YIgxQ_vS6)peqVExlY{O%7sn*Ek^FmzM!)$<@z4>ClxHrS^7Tht(bnLQ;4rZN6E2+=Cmi?q z$h~9La^O7b=DVwic-pcl)_HAGkiS%6Q=zH>iX%odYQ)$1QiARgH=7bv?{C(Qt8+ok zY9B*D8u7CTBo(!G2t(BRZ@Of^E?yc|4{bNpKyedK+w(_y*gWl{we>?1EPkr@e6Q4i zmuL0A6eTM}ZrPgK`v;)fRVH)A)CxUkMnoPPi4xw+nYfYD#FO&!PKHXJHGX9mt*_iy zhWSYKoBT>2Fuh6mn?>$V`wvuW)W0TN=4L@pLy;&b@;a)Uv8jON`HMG3V@4PixO(bs zkut6w5@GeV@kD1A%WWEgWb|uX)yt?OJ($Vc?3B?Yhk03CnU8omR8)eNf9sj!%Aq>9 zb5{so;{Ii^=@YW3&g0l>FrNhO3@(%9Kr+9SaqFNYJhJoJr3VFw&*|L#8Rls+=cIhX z8vgL8G-i`m+!t(BLD6`w^oxHco*M7kzq~mFDYPp0MOjp#nO8nH@+KS)Ovk2N`xH=G zUCp1wLAcY|qE-f59idk{L_d9B9(gDCo)w?=M*k(zQv=MNcy-}!8Q(E=%q}!7+C+Gv zl{>!HS1=NLd$gL@%3X14H^-0GFcUa01l61Il0I=G-^%M|6|AurOnsfu0Of}C#j-(i zFA09|vhp4|rsxv5@hM_j|b;xTmY{{a0s(rgh4&9YOLqsVQg3Bx{dRId8V{0Fvi=nMkC*#EYj> zsIS)L52lNWHW4?S@x(OTppDzL09~% z@A|v=(-Y548%K2+Np3yy7cFD2Ee66T_kI#396OphzYM~ke9b;8&qk|*#GQ6CJ5-d= zt0;7#B2f$-&m0$4w*&j1-dpP>bM$*-lDvhVuHe2f=WZQxPqbf>&N{YE{57fCb)u(K z(O&({hV263ZvKf;XjmgXWARhsnVxR=^EK7AgPhyUtVe#*o+o~XJ8L_i%@RLoW?^%> zxErExE(|ee1Y*2RZO3>~2!c6`qVzT0AS2(NyGMZV?e_}WgkK{b+TBCl#Wm(wIy;!R zz1tWgFQV9%31>Ug%yqUyM;wz12UQj0B{29tuB%s51HY=+ohVD(vD9y6aHz-|TsAdF za|Gis5^zo7jJO^GJLF|vCA%X~-Tu3wybcchF<^MCO!Bt@Li;X@n_%B3e*K1vg!})Z zTsYoK3u3H%9210_ID7b%udJUArVntPt|$8@n=l_LQ9*C?-p$exP1Xf-UMBltk`wpq zqy2Z!L?4E1SynnnyrAas?PRr`8Eo!_s`JH*q4P@vWh5iHx0FUY&r$ke?HhejHR1kK zc}a*o5eUIWa`Bla{Vw^hvvE&CO`sEb{P4whzF_uEpr9`x{fPe8EPH>b;HtcyzCoKS zl-Uos{~cyKULBAT!e5POw@Ms!qqjbifJS{w)=U{lbWTZ;9&kcOr|9WE%NcM3d{0HBWd1dafsXbKteu@@Yy#}jcjW2!T|yEq;YXM)Y6I; zUpcFfpD*1aC4x<1dDpdwgPP3O`+3#=leUKG-hf-To~59y^xABshaSR&Ui7mB$U-E} zo;H@`YD5GWsK--{@wMuo<3$D?!ZG?l?HcHZg6W@(^nB!A@lT3fDJTvR&S^(;2`BE9 zp@{-durWTj^c3_E2u1%TfJn1~4&F=A0sP2(>lGo6DwPsTX}1ykUxX zoq*OmJPMf4YDs%d=3ti#HERXP{rN#xSnov-lD8?)bSO3^JjL!EGZk`zFul`$z(Xnq zU#~hF-6wPTwmhrTc_eQVOTR7tbh#PA=FIchLPa6=h9S_0UJ-Q1;!VumMN$3oLi6>5 zgzNf&k5{@&56=5WUQ_=xgGT=9YuVsX^iRJ#{+;l+j1ozp6Z|Q5SA6|}aPTa?w14cEMgK#dN%vxHOsMTqP%k9@;L5)pcYaC2{6i0c-(q^5_%*z z*eo6+dv}a@^yW$f3b&De_pJFrCUPD!_Xs+=>!LK6pR{K0Fx5u>W;+?-d{a0W$Iids zX9&t?3%ZKQgcq3dMU_og1)`6C2>S;oKz`elNPw~n>|B<7+KjyM&D}uPhcbkm50)%c zx7k5VwSZYO(F8x_YVNk`$-_mXuQoqj1LDJWbKm!op0<6tY-^P{URQ4xY-w=Cjf^g9#fyb zzAkq&#QPJB2JV;LP~+feS-IB$_U|ZCj4Ld`Qe=>Da91c)uBZ+eeDu55DqvZvJpms_-y(PjU7G;qR{KXp`8fO!f&gJKTj0AQNhEVO+xvZ)vwU{Mtv(Va;nVSK|o3Sayx# z`d{Lo%P+4jOVS0Knd#+~6nB_Krqmu2m4*4Y6fW*0@_bzEICIa0UmpC+QGK}E z(JZu+M;^VVdYcyV#c}G8#23E>BT(4ca2XsTe$Za4@b+?fD5dW7G-D$F9q*ev9NS1v zIqSb`+vrrWbalB(>;d_leUlRpVvR!fUrRIEuL1B|VRLI)wME%H#Ra_t!iz2u{t?6I zjydx?{C{}0FmkHLJV(m}4wIKeE`x9y&Ic&6$9v;(x0nOPAukk<2`z2uii7!=__(AE zTU0XYYPh(PexrcmqV!uuoXz}^LvzRj|J1zQA2%D~(ahLivvuO1cx7M1PTq6!dsx(K zbx8m0T#51ZHWMhs&n&xAI}l#&mY%EX?)ZIPpzUuc;m;~(y%(<`y?f7il@#Lb5nwG3 zSI$o+o^%Oj&Jf~HxGu@$6e9 z`N2w)S*bPhobi1hcW@{e0WTaMwiC|;#bWTXA(;#RJNLj;WxF$YpSV_BC^o@BU75+B zfpE~4eoh=W(#ErxR-Q_GQ*`or3!K!kgN4Xy>C@}PpK-UN_Cu>CtXB#>9=rLWGLlJ} zMS^%eTkc73-Oxd)MS)M+xEXfK+!xIN%4D}15#@LMV@z61k0R( zX4?;SJd!+iX+Nh4{Pp9*rw>TuPtIE=jeH0E%^%%N?JWy_-@gYwT@J->rOm%$do;0i zR`xx?pApUyr)5>ACqA9O=_~sx7Ze!muE!^!8xqLoH z416VCRSs8(Z-i;9uL_x?=sIL29f^~`qh4K|e$vZ}?!LP=>L(AIyz^UR^K|z5maxTh=I6rVkdKGKb(iLs=Z{b(n zA>aoVc@wsyBwwfTEz$9VV+{K3YD)CN9dPsM$;(g5`A`3WZj2nQ4>m6cYHIS(8rz5yJcL22k>LTjys;4&n>EPl$)v5dv-{a+vV1*qZvtEo-j%Jh)( zmXRZX-wLi-F4a3{3=!&TwE56F=~-XtZFY7H$E6ppZ2I4P!JpZ+zh_Vnn`)w+-Ud4( zq32ahxx5BSf4{6mc{oJUpIrDFACAv%%6buIhM+ylb}s{O{>=_Y2`86} zEqx%W^)$XQs%s%zCG;8#9r4oyn)jYkv;nv8=R);cq{k*HHg@2eCQb=a@XqabLs-I` z>LVXJgnOCvB{_!R;i2{Vj9uiuz2g;Azr78f&OMsi#cBcN=tGCO+vK7CD>tq#(igK^ zhABEY6rl4!X!JJO>*aJTWDJ&(9#*Zk2L<(cm9;$(noU7^5AQoz2*WL zpX!%u!Vy@caN5Rc6N0c6vY@;ffOpkdOh(s;|Kqea&5SGIm^`|b-1?s$t|+7ny;mfC z+`};*4RNx>uOLzAqGgEqg4r_q3*>twybaBpf=i+VRU*y!2-m%Po z1>DE%nv>!^am1&sQI4P7Tj^M;B8aDGmp@IA`)%T@Wu`Zk5KzIhT?^`;mqXB)+V>~k zSqsm4s=6tMh)4Cdhj+5CHN=9CX%!_pf;fhSmNdOFdqp7dQ9sG^?dy$|toMfBK_R{n zdL4{JdkwSDY2f9^yhxvz44w$-cOI1TME+x~NAhfe*qg{?-bdye!KG6_Ix;=s?V#hk z($J~LVc2%zyZXp~>R@{$n*6eke19_jraw&& zL(-$3<~uE>*gqv&Jgh?UWLja7+1yrymnTDenllw`RrRx73gkV1KH=$aIYXGU9L*G5 zC7-*<8)vg_7~$OHV(PG?4N_gd^pyVy!PCe0!?uJFE|!z~n(8Xy`~=ozQvSLEtyZDA zXGXSgtmqdACUc_ds1OdD%s_|@F&6EwG=!bX3-@t#eI(DkW_uv-gOP-^mF=V#H}m(M z)#MxW&mT@GuigM?GW<}_cu5UGif_Z~_cbvBsW7wRqKk&j2~G@`!=S9gP7 z@~TITCm7YaawzsemciG=ORu=~ZB&>Bw6E@{y%ePkqH_r0AicPKzVA={xa)zZ ze4jGzeN=?#fm?T|=Y5fUP{E6_%?MvVZ92jfu0_L&f|!>ibhWAgZw-_W2p&zNi-rCf54Nk@{7EjKXk5wuAmj#J*PBUlpEhqk5`4)H64I9g zA_FBMJ)Qpj@3J;{xT7V`Qc7c}{WRNDd<3#&SS)$g{18~X(_pYp1(d;EY|Ayacu?q? zxZOqyyT*9hRkEG1MJ$|Gh+PFHQ#8T?(vc(|tleV&k@x_MH`BHkD?nZ;e4@EY8BQAt zn!l#a;l;f9TywiF@eNYOMv?xR)W}`gDRR$qcQ|F-O#U8p`|0xkszyPNILe|E^O2%; zCexEQ4c#t_s$r#J5c_^~sBKyeqL)<4PG0wiXwtx?0sTPG+X*S=bZJ3+&PpdhMj3qc zI%e@jUZ~L)8s8qM176L(0CPnZBm_M@{rs9V_W!6^|H-b2b!W-Ywku+|I~tW!A4&Kk zK~HL39Ch%nMf9hF5Al$i9UDu1NxVV<9sWOh)JP7Ovr3Jl!jC}sNtT!vflsn<6Q}g&f!o{6v+4flQ zkpY%Jo|5@7p@c6RH&SA^bG!MwClp!6 zF7!#e;zjLKr5)9NaGQT_jsbZvfD2f!*18?XVo$%)PV&X>$l7Eiw(aG6t3Z7&20yhc&nIYZa>Z-j0 zUauW%W^{KzLvUjJ3ZFmNR8RlsZ0e3%-v0YP$GgIXPU7XKp(v=Q2;JH*5sbsbo$8zC z5->XQ%t0zO9iRTZh}gJm4Vz5Y+83gsxcemMKk@JSP;pv?~ z74rm)og8D^O8OP!l?Iezcb!q0obtAVa4rnOf~KSzrEvOp7H!&PF_?5)GCd|e6~oy2 z;rmUh*c!?ocZ`SluPAj&>3Jluy8S}+U(!<;T-$QhM<)ud4cs5YNuJJQq{TR8CZF(L zXnD;t%)r5v5TNnR4KoRE)EZbx&fjwF`o+_-IQOAgA}c``!k@oUu>X?7D7Wp{cfx77 zq`+)^>lJf90#_@w z{Fioy;ZOKO%~&cmMBP8Yw?cAA`AUYz+-+qL$fmJ!^PL3F3!M?W1`1s zZ`Hqak#Ot3(wH`Hm&ZZ*{1?|alwcUr>1{HkinV%X4*U0)p&Zy2{!7{ru-VsUeLez@ z0=j~C*C}J6@+9}vIc=O@Nk}+D@+G!e!EuK=&EOSq_qW+%1US(L|kmm@hfrdNo!VPrSZL0}5bjMvID0dJ_+W{E00SiTkfNx$v^6{9Zh zeI&=2nE5tPjzb=D{+d+nq%W0mnW`j|R|YIQUJK-s`D3$kQ+4XN6l!j5-YzO5g+Cl& z63m1{xJ1Dxa^gf3&QdX&mXFK8xi5C+S|{m!hO4W3rkWsHt}Cd5_%tulmw4?c5{8rY zx=jV)>#fsKog8f#Z{`(L7hJaMts;HS z%WI!6x#p1k+m={1(v$W2X4;xdxFczoGnGW8$sA%=_U%M6w=Eq$K$qj=fU`gILg+1I zu&d^eMRnaJY&zTfuQ=5JFROYv)4o{aBv)Jd-$SHFSyeh(cQ*o=bG8(Mi^NxLkv^_k z;e(oY*LRj2HH6Zp<$uE$C2**=v@+sHAhsPj-F9Ts9b98!kM_%}z`R(AN6Oy}LKo6z zg}IdR@POQ=;xA!f+IUl`DQ$q~RFY;v4zbV~c`M#;Y6i0h!rFRKWG>=)O!fbJlufcv z+piq4#&O0+o77{};66XKS4LMH{0^;NhGo%UP>Vcvo7Nbm%s;sK{)mHtfp4y|zsH&QFERy$YgxRyrWFO1BN9SGb z<3zNV9r1i_6oSCwLm9QN3CEdju+;fI*`JMX|I%<=6Vj$1KBuM|AWUSKj>VdAZo5oh zISvvZ(oHp2+b`Pi><`!OH6=NRkYC3?P6VJ;TaDiA); z+4tORYVcJR$sB&;jJL1920z}Wg2TBbi`$tKG3FXNv5@Wpr=JnbQ-t%CY-82?YljB5 zym+7Rl971duicsCdZL1d(Onyi9vWo7s>OMQ+#_xce%mJ@AP=W=)(1>JdqY9PeT$o( zBIYB73rC}aF(N?H7ThYhzo|TxhWx%q^g}=TzLCJC?4`VM8qyo&bv@|OE`zxpZ;!Ox z4ZyjTr7x4CsTervTmR&tHojZVHc%7aS@GAflD#7SU|)04Jx&#jdo8=i98Q|Tm+Ifu zhr(p9B>(nYz+bt zyFXSTJ_)>>ynP&>HOL;wW?s$K2V9JkjLA|i*y-{5yDsspP@C(Rahl4&Z0BQ(7o(m? z<>(l-<_N;z_ICQ6JIVb@T-);(+23nVUcYm<*#iMPE)L!yJUh+(LY{G9WX^N)-A5YY zfm1zsZ|vUg7%Z%Qr{T78Cj5)UOWk|okV(}se!$lmR&O>%1WucPi9J&zFhmZoHqA&h z?RLd%bkIiUK^?pvntHFNrvf!U3WvB~5eR<%?M^kBJKwN4&Bv0d0CB6HjL|1oVQtc8 zX1v)RPWKuqjtu$2BE50@+*6Xr@6)bPCH=5ax#>W`CQB66rzPqptD%LruAG`v11tdz zF$yHdUZ;F^y5Ke0N9J`4vbqwF9M!gqiao*DS4bP^#~yZa|4qV zpL=wBW3cMb<^O)dhJ)3pOus{CM%+{JHkiB5#g^tQx3Cy%yEwenLf?bVVAE&$hFn^|=>J3u_ z;p(aW?sd~bL~2ca0ZkIl(fAc~k@>Hkj#hA@0pUR_#HD?H=81N9OH;NqU&OfG_Mq<3 z#^8UYr*7+5LTq!-NtZfr&|qn(RM!wG2kZ8y9egZ9z0u2&95gGau0uv63;myQWP)8bJ>!B4ir3zg2GN#?%T zwbc}7d<2@yBS`-G>?`}aeg|~0S~@)|bilAOePw!{AwDgSayxY!U?QBUZPEP-;SUf0 zV3@SR@WZrCX&K&l)OGLhE^=QeoYJG6yRM1v_B2A#+k){QIeoQ^dVt|)bV2nSKb#YJ zGJT}f45QlmGg}HsFHW~TZFj#JBHml%W)vzC=j3#I`#-|Zs;v8+8SRFQ*2u55M-(vl z_f`mBls-zRc1<5sR)t1;>dk~|c_>;RV}EQJO1L_UAyFk^c=+GC`p-{fAEwYUdqSJ= zq2#4Z*_4PUDtO_bjIAc&eWX@5ehxs~t|zP&B=4I&aa$slGZ->FgZ|Z_GE)6J^pJFW9M1>}|X2fJq6Ocz)uE>gi3*m}MimJdWg3H1p1|lq-pe zDzQSvt_M4o$@li`p7HNp-Kh}R6yUz2mh`p*6qEL>#A1e_J7Sc{4O7>zDcl@XhJ~`+ zFD?;Lwl$)$^rIfd}Ro(#sYu#?y*)DL7 zPvvsmri4b;KR<685r5L!fr`ypK&jM~jJQSVbP78I@Ka?RzDn{*TWb<~ZPX0G zE^WF=o%D0=P@RqbLiTot*4|5nJ@bdp<8usiD;{_gLr1rD!Vc!w_c-nLaf113>N%rM zMRG4l)E!DCdAu8)eRfp}cxCc!-|tP*(7nQ3TtoK0bfdq0)23taK2m9~<&rygL-%sdVK%4;=c}@{no$Xj)#w0D(n;O zA)Q)tta+~|nkJ1pkC6G=Z)-Hz7JDF4p!MJ}!ZSEO`fr!nnmBZW=wjer8sXV(PnvH1NKEh!yxKjg z16loS;VpcMI5D37reQz}b)(5@7WGDWef>jYJ;_^%+-1orFf~WZ{5!WI(nFGLc3Y8T zzlt494vjDTNq?j%P_5fG3J3G<@DGcH!Z6Y>`4;(J`{yiCBqb?=`po5-pqEO7E9p_7 z7Uv1CA=RmnUUfJm?WaiMPr$q+`^@oC1GpxBm*ScZM2hfOTwAF!a=sKb{v^)<^>nvC zBmD~G*WbQ+Bi#>Dm6f5VGu;p=?$OJ0-wuoK%x`jhGJ>|rXy{#ab+{?5-x(*lC%MfX z-9wAy&kJ(6eka=-wc7;^6Ef5>e6VlxpKQ{HH+x}cPa%tEhu$?V{vhZ0RTq6q_50K|4cu2cG2m8ighx|!d!q(K5##wU-JO-})u#r}^y3%AAxsm(m!54f-zLFek!V%kdQP)WiqL@zIBwn8Ite&&Z zY*8XTQr|l-$v#xWvc8D(lQLp<&j#fekUiknfoV6zII5bMrx052qtFM zi>{5vBJ}dZv?~K@m{QD_weKixj*vN zxV0tlC_=B@p)oI(%q7{Zi_UwRVtrU5E!~IAcj-k-IxC~GXHr~Q_@XxGZ89%DYb3oJ z^)A&pNe^^0Q9Wt_>FNCA%eQ?O05hgI=eAle3>+Do_f!pm(5*v-{iJt!*49EfZN(U) zhg;fv(Z*x84s_`#>;Yb>-ixE3#9lFWo16?E|v27EimS zah)aovB4k0Q6XO z8TwI7d?o#2q1Rz-eVRDw!g})q*?0bXD|-DzfC;op$t$+?|c$8fW+D3Yr5rf!Oh z;;dEr$^0@i*xu&Z+{C917eVcS{Ulf8r_(-dIb};c_P%VY(PZu**GYfnvJ(2onQ~H; zd@$y@t9PF^$*oMtUg%j-g74wiS}~;_$f=p;n%xtGlJBl%Kgs9+1p|%BJ+T+y0WIyv2jD2EcpBkoF_VuO<{=Dm-NXT=+xWked^KKcWmi%KKS0;I5ZiYWxLY|nj zVA9o?4#Ak$wb^Mu4Ycgns=7673oi~$70Cr9+(1pRDEWDwy+r+LdxZyRU(kM5HC960 zO8kv^@_yekkv(_MHwf;3`LA*9A@{x}LGxICIqdw8w)mc>CXP1!R;IQv1z+MR>rbI7 z$aop4!h2F5S=X=51YWbp_k4q{MinV=%S7JbF%3gSMq~ZO0dM4F+CKkglZ25unAWx# zgSGawqy%dS?p!!kdefbF$JM4B2K3Z$vnQyaU&q)}a?Veq5P6-m z;EhO|&eT6firC}4@l9(Z7~48_&iWou$3asTITO-z6zn_}UY@3c<=;OB6?klMz(w}O zoR%%|f~o%vb`Qn3TJeE!29kd&9E_i3i^4vZ7fa_r_9>mKm#;)dAz8-2d(%5_{8ptI z_94$Trj(6?ox#aqjM!{{NYx1{FCOs^m#Cm~RBnIjlnXfiUAx$j77OD`T)$;a{h&gT zUfv~5a=o=Q+E0%IJolw2r)k6>`YLn&hKM1w_FbK!pa_Nmr_ql>cMaI)W)|9$-^ZOh zNH);W9O3DjeEWJ5kbB;(pmj?qDz~&$#EKiEbKAs|FSGLaS|g_HC+~_>5w(Opydm%v zxV3xVbSTW<&xmFzs3PaG!ahZLW7vCdRq(tN3GFqjDzUU^Y&mz2E#Z?E$iU9;MPV>F zdNk7+Zfc|b$UEsl5sWuX%S0Ox zuV{_m@qKf|kNNq8^zyg|68UnbW6DbL2s= zb4A92@H96g*>jgO1;z;zzg#){S!Ln@Jmo6vSx-0x!%;iC61D~cit=?+?{zUI6!fj3 z&Ie5&c1GWP7mcnWZY_buK+x@=a>*k4gHB~B3v$64ZE{S@_3iM zRzK?NiGSa_znCY4qtT7eA*sUzAJs>XUD#;|+jTd7B|#M=y>$N-B|v&$e6wmQB9ahV zpQGpDH$~>B$hBtDH`QHFJjf{Eh7}*4pe;M(&_*xwEPlT>3bP-Lrrwdnb_&rv!JXEi zicT4!FpI`pLyM)(cH(>MR2f^BR0Z3$1Cy0xPGb|o_0Z4H0T*8}+whX}OtTb+R&T#N zcy`qmw|(-(!L{CVgQrNoz~-*a%UKZ|IdeE{e48Hbtmf*lba~;f=GE1iI(yjWkDjQt zG{QR5AE)DC+Ayx^T-6+u1+zW{li*>(b&g|eRLNC@5M^#{%Yq^l?V7f3Z@rB5KX2dI zvjoHDgGNc31K}?nH@{0oawsVqRI=^?c1V14vO&+1}qe4gpyR3k-|^H!#R`Qb3@J=tZ`qyE1V@cD&CW;r93KlGMO$aSFQ&S zKS}CeJ4_%y->AK)9SO@1>Fo!-%|Ic_VLbX$84XiO&9Nhdqi;3mon)?xC+{l191}K% z-;qNO>XD{MbX#%!*lrB-$Dc0_3z9tB={4Jj310ZSP{E{=FNN;`>-; z!WS6v*+uzH9Ny&=J2MZfq1EO3Z+bs>XpS(vr&uJM`#aaouDOyuOm5S+8A=~;o_+g0 zqb&l@a~z7a$h{}%tYU>QzY}O;4W?iJBsm$SNrw3x!nfgA8+mP}4zD9oVVZ;!aw8P$wBIh;272s?fK0i1rj+}yT_Ou7XFyYqspG+z7s#eh46j$`YY4=KJBTI75 zuXJoRArqj5D%nbCLqQF5#5a%e9jP(K5ck`Z`9OID2vG->klsFZy#5W< z98rYK_{_0%S>RTNns6-1UtOC!da6Uj68@#XU6XwYXR0di7E7fb_-?!v7SQ&F){&pT z85?ZK|NqXDU4<&>+WzPFfI>LVxZbE=*rS6U^#Nr;4#FiY6tXW}x&+;O9tL|3kUmTP z#ktI251iZCof{Ay45g!P?5C?puhxjAZ$HVO$T$dxidvER$6B9d`YUgYcNQld|7V9U zo?+r#{ANg*=l8wqN&H^6IJuwj>Yo=fspigVfXOV;Wm};IPQDoP-eqo$Ro8Lr@JJI# zZA8)z{?$UjGxJY$I^5^o^bJ7bvPBLA=HVo#QZ) z#J8VfGlS8xAa&4JUn^9hz#(6LiPZ>&7e7cx{r192N!|C|L4-rY-_I;_PXhZm{C!2D zrJ*XSY!a74?weJ`_BYa0!SH#FH9geJAgIAbR8)Dr`LbNJP!e3vAA_EV0Z{WcM=-vgZ?>;6B)gFSV_ z?cF^~obL77=SO@v??Tu2dRn=XP-K$TfpNri?N>mo#>64Q6&0R>b2z z7}#+zrv+r2P( z_g;74H4V(4%S?W{PYjVe=w51X7=yW8YRL7vkCq4zKHOHcGb#=9JGb>Lk?1z(tvFKGlk8M2uJ(d+SO}Y zVzBf6B7=NB;lkwU(K%YF!{@GZt5KB(s*iSZjjP0gtIz+Hd@{)i^N8?|GfN>gGA>78 zmo!|h72O00*J8O~56uVC2W)4#`6Q)H0Sr&v?+JTLphCIrjF7J`uKKFF2PTpJ^S{@6 z`h+{V(5b&N*RDx?T;aQZlAevZ;& z6?{dP^lpFsJG-B7z9O7-;<^aJE_C;l{iw-Z);HY^Vz*=b=2ljKi8p+m;i}A!v#F zdm6WKjE5j3t@ZLBeottgDeFvB<|6Ue~ zCVSDbg?qgJ=L_s_)tcX}g(t7#RaQ#F$y{!Qap{C7iay+zN@i9k{6AbWBYwK4w;XJ~ zlYE=bfR0S{nMjf|(AalPi}a5;brn)>DB!c*jPv#ecfvF6 zgV)&`d!IyKwySWUOdU20=dIdTT*3;eDtiIA6VOkMWNe^^E?S5)RenUs>|L_%NMxR)L9d z@($c^ppiAf@nsVmNkMYH>zmwHdC?xBKQ)h^d>(}Jm-xcBJodsdi*=s&o-!zC%;ncF z@B!7381IO&4k~EuV`)B101}H$>trsJ_2hFvmOaT==xWlwmasr<^;&vsoCd3;+3wU`-o$=s0Hjga6))gn5ZMY&aD~2f~Pl=!MY*PvYefvt&N# z%B}?USKs(vA2LjPh#e1 zc=5OD$Fw+ZJqYq~qmjgaV-wwag2FJatP45#m*k1jb1y)}1)1K>o4bfdaVXq_Ip0A7 zh4 z#1AC$w3@R+2Ez1KyuuGeptTRzL2pA`WwK8_GC zRukB_8hkSNEdy2V$7X_g+K4ZcQ8;L-0^UE=+(W&g=$uiz|B=rdDoQSoC56a)g5zmz zA*VG8jRd=@Dikp+C_Pub*$K}&=kC`ksbcSE)kbFvEj+p&H`@7E75pLR|GmBB2;W@Y zbrWA1EdLbpDYQevI$aEg9?L;P5z=38Vk+MhR%U>7 z(K*ovm67QFZ;OVEnG-$=j9h&CTLPx$E3#`0Nf30s`RXaTk4$nE?6)I1*IRrtHGWc7 zcx&pxSJ12n8=)Yd+=Y0gF|iL!-0~y7LAA!NY7b=Wo4>@Ca}`>ucF9ra=k{3Y zd~-N&uOYeT={o$Hki)dW+8KJQt2pzlXK=TxDtT{helEJ7^tTgM?75N*AX^YXl^>&r zi@)-}Owtoy)ywuf+l5>)Q)FiKJJ1|4=7$Rle~Tku=jHG@!mVSBao(V4Cw&*$tLrB# zLP2R7Z(^hw0i(Abv?HAnc$e?|JcZ2PWju#(g_#ha(IyqnO@yuDwAl%M*1tH!pA50s){bE4w zx8nAZle)xf;Cr|rzMMiDJwJYi)d_rUNlPYg-hRl>4F!LQq&D3Zysq~|7oyI3sm*HGcM^n!FQr{tH${xSX)VU z^$7YwKsYh&(mhu!Pd8?A?$;;t%No@cD?4(YaB}IJbp)T%YN27W1uB2a>57ruK|~Q} zxzrhMvY5%8TM@A$ zeQ}CKyUv@QP|oQnbY{_nd&gd-E1Y()3BF+M#BBiMsK-xM)?6TTQ`m3kadXmZXnuS` zh|E8$dJZ=jh7dn|vr}%MGy1kh?Hy&dK+HL5nZvSVF4eV~x0OvF+l^$|Xtk6et;Tfu zYb|Ry0MtVM6*yM&=OSC~J&AoA~`egrMFITds{&4&? z*-OPAVmjU8>4)qsRe>Fclpv_$#IGsff&hP+z#(ojpPJO*oF_dI^P)`KKQk_vIdQR( zevdjh+PNEYHKTEAP_ibEUme#~2CM7~NRFiJRb)85E22W5%RT+*j$XfePbzXku{Zl- z9Xshoe5yWE!+yXGCzAuavJ+I{ekRmSiac*P(LQv#NCLxst6RG4O)y_SXcJ9xM2CIc z;sk{ipm&JsG_59i?>$Ut_IplxQ`$zZiL-=D|A0krh2(VQgiQ7-(MKcwX|%OXf(5R& zf4|SZVgND4D8;aJagBAN&2Ur1Ora`hVZJ(X^UEDLzybMHG>7+sAh7 zz4zXG@4dHB$tW3-L`D&!kdP6|Rw_wlWc=>$wUlO>$;v-zZuYdb2HG% z(S)==muVse!r)(xW~yeNGt3;zo_Lcw(U1+;iafWkS9qx}l$+R}XJnSFE zb4nLgp@)p^dep21^ln)G;_)m2#-jqJ$7qtlk!DWv+kzqN|81OWTpt1)Z09~4&cOc4 zEn)qJdSwt++jX|MYybr+4$TTm`1h+HPLe*6g@OMZ1uxywM*oB9+hfQv)%r&NKBL|Y zSg7w#@%;6M0{OR+VW?ZzdHy_A$m$K}bEqrTe_^h%e4*^X81}S2Ve}L=ih!rlTNIOl zYRD5%{L+HG2u`03G9@2txY~zyeyP^Vo}Nnr>CN zBMkwn6bI-a4g_UORJJ~wf!vbyg85G?7z(*!w@^~7hh_v7d#K%tnqMBo z8`S@W7-El6;aHei3&52vs_B^)V~}|#wZ(J|b2d)ib;Gzn(rZzRJM-BeLRBLtErf}X z%)3d@X|sWIvB##Cj!Q$lnlq< z;55_!tlx7XGVeHp->ZMsWbF=MQbFc&RZj{2FpDuz-mp3;Nam-SNSR}0`96B!WEu8Lxe~t?5 zpE}C!%Mb&P^e-$^REL7W%Fa*oG)I^nf9bZ4I;i^cn3XE+3$aD9_ z!@~EO!mhja;PY{-ViGwHfgj_jSbU^FkNc3`#v=kOav4f~ipAbZ=6`k172IK8;AX%2 zGX)5LPrEML5(U{@O?ha+fc`CeGyapVke+R`&A@L7qlVQKE?)@1vah2#R7Djwitla;6n_>XZPFHDYtSgpq`H7Z&_OIa64^v1k(%q@1EWjUY`?~$_^ z3V;-sAH~ea@%&moz9U#f1Y-v}iJU^r6ZSWGCI7^J%DCAJqLi5T96j8Z_QwH6yGDn_ zp2&fkW!=9b7R*aYbF(Cn1;F4&AbrX}D8RtEu=jz;3myF7CKIQIdt#nzL{AIYKB_fL z;^_usUc^V6jUhmCONdIL%?%u5Su%Q%H&ylWvDkIYBdqx{(#3CM-jylr#;Z@pKsN62 zb-Y{+rqV=rIcJN|f4Y~J%IyPVy5xn2v4@Y8t2F9TrUDpr9XBHza0br(R!xvSVxhS3b{yYp*=eX-2$Oa zi{W29(H4XZUoT~qX@gT<)QSh;JhkxwtPJbl4TsaCl`?i>u74;_z-=} zYy0hAWnoYFcj=w_dm6w)(BU4d#k|kLfezzZ%>Vv4vC>DT4%vY$m&3lHZZX{8AaY0? zBnujv&qP?m)h8kCybcl&oTY8Zb|eByPj+@!3hF^MC4H4&z9MY%^}pPD>khM0OI$nb z=x1gW^$qThfTrD1m6zzdKhk-Jb=MnvY9@8ITSJgrZvJ_*74NIO4@Pn~?a^Q3e3a^I zIr0G4-Yl8n{z>S)-^NI{1qjAvJ+va|fQ8U>^$FC)IKH+&_X`OC4^TM5d|4k<*J#%m zinU<80OUmIi>%doL5HMqzqY&=EWCd4D{R;XBpAr_w>F7z=E^aK>Jb^} zwrN>Qd6x+NPZBvyN!-9Fr-#*=MjjG!dnNh?ebXG8xs@?E~ zDh=i&vXr}plVIy^S6qyu29Qe}oVUqWf(Xl|1vQFv*zO zuOJ9g97*|(^jbl_SYt1%KKjWIbBXle+}=DcduemT60RR^vvxf03Y4ZtN>&o|z({}P z;%_!0tO|zNwB|~IBJJ~ZS|&}n-fH}k&q*HYe`;ROdM1y2R=S4FTos6W=SUW<=m*+0 zjt7rlK_Ad!yh!Dm3tWBueJf%JbyiO{?HRnka+tSD4#$f^*E^BVolPOYOXRzI@0KOp z_-0s{G7<QdZ?%d5ZsmvvAaIt4i2H3-Rry6CVZ@mwS{{?*@-$Guh6?!UE`sxggjVzMcwv7F6>Gn<;bm?NJtEpC8a7x-avoo)q(P8SZ5IWT=hf?NSjwm90w4O zEGiVwbWajqN3V}}P!J%ZFEJ%=*clQxRCdgj;=r;{J$>~8akp-qT?{w?WJ$ zEx(&_=e}(Xdm$P584TDH5-2Jw%cu=^PKdE3&YD3-#gWItf%f2a&fN6gU8JMqHl$nxx zP=6Pq(sq?tAK^Ysh0t{db9DP#cPmejX~FNm?B4O%CtU8CS=LC4xrL1T`P8HC5PX=H zFOk{=WZ6mQs5!L2UGL?qTgTBK<*M*c-T}|UlO7KhS_9#>5VbBPzdw9kF-G3=xt>WqA@4!xnBYhSO*f^^-P@wn%4P)8G1;JY3MJO(1>g`qC+ zg8CfgDDne48{$tC^+p4~a|_$hS4lAT4=&mtC=Gg*k{J|a-r#9}VT;w40JNl?OWD{H zGya=nAwmy(?YjB(wX(HAU7puE#7zbMDrlt8vzfr&{pL$)nCp;x=Bjwo%o#S1onE?| zO@PznBL=6CN5*r0=V!B=3Y_}VsNQ-r9BO|GtbhFz4)i>!ckFJcfz(@_eHjDZK;0*y zP5agh)HWr$=F6mkPS7fuEDG})(-fV2%IHH$-TuC}j{HUQu;LU^dtgiWFnil68uv>x zw`qK}fT!m|W|giVuo22=NQyKdBcPEb-pl*azOm zNHV2f@`Vzeo%)BRF33lYkPAW%w#5}+=Q!Ly58Y#-8F9tAwb-~V>X8rF42H!#nZi6A z!@eX7D;wOObPcAIq5e{Gwn*}^DF|@?${Rr6g^iXQ3C%lqu-!`WI*MHQ^{)Qx-W&4p zVeng?Oh^QrY^}?A9_Rs$u3eXfklRJG#rvPj%{U;p*!zi{^h^h^e zac|%qOiOY6f`|gVOS`tlyyOk@&vIHm_1Hn~A>+i*PYxjW+x%Zkx+>0fm5n&!N0dZ`9i(+q6ae zVRrFDMiv9^{Up8Rg{19Z+0fSs$0x4ZIfLq~S(#i*UxZ8A5jq)o%LqD0RAbo6f6Nd}(m7u-S{u28HUAj&AJ3hdhp!b8g*V3t_8 z>fjp!55%@gRC1Jnaa@o&?2Qy0Vt+yHjyluC=(v~on+)jKOLl!(iu=Tpa5G*WD@Y#I z4or(8!1;*<(h>A!<>`z}>Ee8`K28LWWALwi$2B?AJ z-wtUSFMpU$T3xsD*Mp?J7o7JnAJv@7{BJWl1ZruHSH8BhMI}fOTc@c!MaEC3XqmV-@X4V5twY(7asK) z!*?}t1~%GY=-5Ah>&Rgz*!0S2sj5o_4J-B=li$!ko66vjYJk1w^_Ml<%hln?@JMN9 zFmga%&&PYbmIX@DB2`7qc{1#ld$%?q$9b|lr{2pI7=ON~l+MSTS2Rmv0)9@#-*<;D zc&b3mE*JNfSPHzaA;x}5$6mO@Cw(f=_wmT~&@om#cQ&ejzWH=79He#yrS1&tLWJoi5U>%t9%clxy4;d)IURsNpUv-#y}nJR)ec^A*tZ;I2zqT+L%|sF<<7NmB#h z$+>815`#StE?I^3c<(BvIP=gAvM^_JxFeXc~z?Q3B#Z3y1*&XAE~cp z58M1f{0baIa4OVT3w9y^gDeN%an$8yE|lhT$`N5EF^#I7#1!XuO-|=TBE&A#^G;5i zLu;`3QTdewn3#_gWNp-gUms0p1C6l{I?`~0-qsHK)j3qtLh<)~rCj-LoHCFr1j!}o zctHV09kYCs4II02u^_Ni26}t)86L5i09D-n*qcNz@V7fAT`wMi{3h3Y3wG>>F;p_3 zLVfw!{4V{MrZ71hup`*^mn}^&>Wf*$w#Ts;@l5ePGU+5meD)!EL4?y2dRQ#O(GI zsZ)t?B`D%f>ycDg{qQ{^vD5;{7jIbeMMeSh=5O;L%rkhKJbJe1hItr;bJf0eG4QVO zc+hfHAc$+1rnv17hUr>6^J4UwDas!_s5lq|Bsy~?P1@MM^5ol|Iw{^~F0uqfzb8PS zNHk%A#S`irWIBB|0)cBu`p@$Pci=b>fAM&t6Z+Gl#Ge@A`6@v(^TCz~MhsC}E~s-5 z&%O;(qH=;ycRGfi-Bg2S=c=CgyOAI|{b7LnIrc%k;f<=vl!E*OeVc4!Y0UA|I5O5t z!36KnOD1LH7sxE$kUbO*`ERE+Zqu58>-^27*%(=%D5U(cWoZM;^+J5&o$_#c((lD@ zLIe=>j5@NfD#NN2clW0!=CIiMC$+2)dASGVoi}BSft_Q*ZVd0gZ}UY$R`I!g)-3jY zKEnZ;iaExsPN)LS?ErPH78f|}XCHo_#RbZxq)#!WBL~TZdVy2~xkJbHmUqw}b2IMN zp5+quZ*D!3xkN99T<;b0oLNmMPOjyUPZI;%1roF7S~u7@nLAs{r3yEX3@zGf0JJea zbt4@1f|UK-x%KgCaP2K!cOUv4&o{jsbwYn$t~L*uYg7Pqi6!2rZVQ00V~2Lw95EM^ z#GfI7d&F=J{i;${yf>8=JGb2QfHM#C+wLIe;I(S}i_@jX@Ng+}!Sn^@*v$CWHjr;Y zU^y}nW~>Pd4-foDqic@+GUsB??Dqz`$p1LRFbC~+!tP9DIlg}x_Lx*#U4YIddW610 z4TOFs_PC(G>9~sJ`Kub3bI+Veks9(9g21WQ}19$gF6ExY^E3jP|x?(M5y zkpcMrj#%i9Z14h^@QkvDBl!GTj?yxuD}&o<(lSD)GN`;6U@K5_26uDqBZIbqK%zc+ zsdX6rp-!68hSyb-Y&<3UbIy@QX^IP$Wwy~L+T9Rf%GYz?s(!+3i5 z=?#28$bUGNXC*8PR4!~+LLPa;G>4dT%@NGkZtt2leOH4^W5k}8I9U`$Z%_LE3;^=F z6&{+un(*PHocT>wBIv!J{_@um`_rk~*^9$5M@|1go`Z(~d4xBcm(>+v&y4@v6J|+p z6q#bW(xmtJ*6K0stHC#*R##VjDbdHA*pav47}T1hqh6N z62CeAbnXK7+hixx@>)i~95aZ1!29Hct_Ra;`&;ngPW0KmZWrVedvR}~PWZF(z`q{! z9h>^c^(uLVfKgP*s`4KqWJ&gCeD{@rHuB6bL+P076Bl^=zx&X?;F&dPyf61YeR6ho zUK^&YUY{p0*g(7x&4Oy13z$@2eScFGxemilfz_XVpf{3*Oiw}&L^EG=IGpo?&2L)w zK0Yx7qh{UezQ=M9D=!>2XN=tTgJwIW`1yHjok}}xp#`dJ?>49}>cYZKhe53AeHc;L zQSionON6(e;rl^1sFY6Z{$ei&Sx=9?N$yI4quVsThbWO#VG;0N-~sNNFAn#45~+>1>=AQa4p89AM3lf~Dw^7AQ^+lhiKRfPwsikm(fm zMQe4?v~K%D)lMqou?iyOl<~9tYYD}i$o8dP?NCtqR8Kd8dhXon-#3OshEUL5%9w9t z1;@v_gqIQn;JlY<6`7PO7+i4HHC9)G84(S!6T|59c(v4I*6aa0huM-1UIas@W2j(X z7y1!rjNbJ&>%o^-nI5#F$a|We`}@(?5E2oIS{Zv z-Zgo}94sCrsh;IV-&cixgX0ku_|x6ex@jZ_KDSCqG*by66{beXYM2R&mYiI!64@aC z)yHBEd#rv7@}_byW{@?6}{hTpA5uJ%xGg11v)72VK#hJ;F#VNc9Y#uk_MSWEJ&-3lcIu&9 zVS5y4%CBBjPl1N)thL&t$&=w-v++FV#iX$iZz1fkXlHs|D#=oBKyOqql$f zU5$(cBu{m{r;Jg7aO>eHk6Ki9$%;P;C8T8w$L!pNDI3?{MmO6dp_i_tqCu$J!_m> za0AARJk1TrJ@kvZ&(}{>gt(n!{<*k^;j$?p58ifw-#&`Q;bwa9)w9C%glY`7zgLAT6!@}cAqM?fDMr*5uwrroVpg4KxoRHc}AzCU`; z_G*VB(9OH9JWaI&=BEQEK6S~0dwoms7?~n4XFiSh`i0;7{a!kjrC?AT`#IQl+Y>x; zpHWl)jsX5p|L5P3J8&W@VSU2R3mUlBPG==L!&Ca9V}hv5Qc?C*7?J9N%;)_)*M3QZ zpz!{@VlQXNzarI_L~n%LSM_TVrP#|6mr|`*;sI7cTFkLHUwJY6tITj>?)l^UP!fJT zuQ>F1f1;k^U=N&=_`TJQcn_;*i^24NtkI48U14VI*W>G+s<33A?QPK-43?z21v@X1 zL(k>H%x)$JPnA>`Q{pusU`~;~Ttpfk>fA8uR}}+Y<1pTUD+bVN5}f`1PB>iaEIplP z;0lEl3}egu(qKIy{fUfE0p>@SIUMrTfL!B@_7DGfkPm7M6yn4_hr`3Pj^RX*p*-ju zaX6`+~lmOC+;@Q}FYPS^Y*QrwOjU#ho+m$mLmC z2;Ck-pRhLBJEI%MKwgvoJIR9g;~8M6y&qO^ejM2+7*9yb^* zxqFku8s|cG-BVn4VsMI)z-uAt4%_kP$d0fnzyRfn;FsE!fyGGKE-UA6F6yV<@W~enW4f3Qq!ua(6~i zr#kb+AZ>IcG z6k+x?`RDuIJ|IcA{Fw7`C;0|cCVw5F#4SDUS4AkkheRW)9t7bzw;tGwz1(c{ ze}v}NB;n|-S@GIFS9tTaW@>m6`!h~Qo*-@W1pBZp4_a#ts4`J(`Et??3PYCtq$nd` zs3d{w&pvz(IGO)_K%O1ht>ze(t8PHKC7WnjV+p4mJ5%-lM!~U-tqF-v^sgLB^buAmY+GZDb8&${(BW8s z7uA~1`7Qx4VVin~)5RLPrO&DGM(cr7=pb|QFDZC3bxi2oDRtzUY>HHmsKNuAn4HJW zxTosVFC8QL!`(ZdYhPpTE8f~RBB&XAf^JiY{=qqZET>FVKgk7X#}Y;=^g=-kEQt$y zc#nAKJXnnTe`+oR2ul%zy3N>;-$BaAm!+=R-)soYQRnr$E(L)1@5RXnqv-eN**B(g z2J>d)-_7`cf0v$8`^(f9YOEN=&&`?||^101|6@X09Y zz$$VidTM@~7NXC1?j85zlXFDiF8tLwqbdi&3=7ZC2us1MQH|_d6WCiEFFtn2Ujqt` z`<-ghV#XetD+0BeM+u-toctFbY7GBfh(@^ z*VY Date: Thu, 22 Aug 2024 13:03:36 -0400 Subject: [PATCH 2/8] Formatting fixes --- chandra_aca/dark_model.py | 2 -- chandra_aca/tests/test_dark_model.py | 27 ++++++++++++++------------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/chandra_aca/dark_model.py b/chandra_aca/dark_model.py index 0d0453a..9fff944 100644 --- a/chandra_aca/dark_model.py +++ b/chandra_aca/dark_model.py @@ -391,6 +391,4 @@ def get_dc_exponent(dc): y = log_dc - e * t return (a + b * t + c * y + d * y**2).reshape(shape) - return img * np.exp(get_dc_exponent(img) * (t_ref - t_ccd)) - diff --git a/chandra_aca/tests/test_dark_model.py b/chandra_aca/tests/test_dark_model.py index e33c813..5467b1b 100644 --- a/chandra_aca/tests/test_dark_model.py +++ b/chandra_aca/tests/test_dark_model.py @@ -79,27 +79,28 @@ def test_get_warm_fracs_2017185(): assert np.allclose(wps, exp, rtol=0.001, atol=2) -@pytest.mark.parametrize("img, t_ccd, t_ref, expected", [ - (np.array([100, 1000, 500]), -10.0, -5.0, np.array([ 171.47, 1668.62, 852.79])), - (np.array([100, 1000, 500]), -10.0, -15.0, np.array([ 58.31, 599.29, 293.15])), - (np.array([100, 1000, 500]), -15.0, -15.0, np.array([100, 1000, 500])), - (np.array([200, 2000, 600]), -6.0, 2.0, np.array([ 478.16, 4299.02, 1399.40])), -]) +@pytest.mark.parametrize( + "img, t_ccd, t_ref, expected", + [ + (np.array([100, 1000, 500]), -10.0, -5.0, np.array([171.47, 1668.62, 852.79])), + (np.array([100, 1000, 500]), -10.0, -15.0, np.array([58.31, 599.29, 293.15])), + (np.array([100, 1000, 500]), -15.0, -15.0, np.array([100, 1000, 500])), + (np.array([200, 2000, 600]), -6.0, 2.0, np.array([478.16, 4299.02, 1399.40])), + ], +) def test_get_img_scaled(img, t_ccd, t_ref, expected): scaled_img = get_img_scaled(img, t_ccd, t_ref) - assert np.allclose(scaled_img, expected, atol=.1, rtol=0) + assert np.allclose(scaled_img, expected, atol=0.1, rtol=0) @pytest.mark.skipif("not HAS_MICA") def test_get_img_scaled_real_dc(): test_data = {} - with open( - (Path(__file__).parent / "data" / "dark_scaled_img.pkl"), "rb" - ) as f: + with open((Path(__file__).parent / "data" / "dark_scaled_img.pkl"), "rb") as f: test_data.update(pickle.load(f)) dc = get_dark_cal_props("2024:001", include_image=True) # just use 100 square pixels instead of 1024x1024 - img = dc['image'][100:200, 100:200] - scaled_img = get_img_scaled(img, dc['t_ccd'], -3.0) - assert np.allclose(scaled_img, test_data['scaled_img'], atol=0.1, rtol=0) + img = dc["image"][100:200, 100:200] + scaled_img = get_img_scaled(img, dc["t_ccd"], -3.0) + assert np.allclose(scaled_img, test_data["scaled_img"], atol=0.1, rtol=0) From 183e431e59af8cb909b68a1a6cf860ffc13f3631 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Tue, 27 Aug 2024 09:45:35 -0400 Subject: [PATCH 3/8] Support scalar input, other tweaks --- chandra_aca/dark_model.py | 43 ++++++++++++++-------------- chandra_aca/tests/test_dark_model.py | 15 ++++++---- 2 files changed, 32 insertions(+), 26 deletions(-) diff --git a/chandra_aca/dark_model.py b/chandra_aca/dark_model.py index 9fff944..d247b94 100644 --- a/chandra_aca/dark_model.py +++ b/chandra_aca/dark_model.py @@ -18,6 +18,7 @@ import warnings import numpy as np +import numpy.typing as npt from Chandra.Time import DateTime # Define a common fixed binning of dark current distribution @@ -344,39 +345,37 @@ def synthetic_dark_image(date, t_ccd_ref=None): return dark -def get_img_scaled(img: np.ndarray, t_ccd: float, t_ref: float): +def dark_temp_scale_img(img: float | npt.ArrayLike, t_ccd: float, t_ccd_ref: float): """ - Get img taken at ``t_ccd`` scaled to reference temperature ``t_ref`` + Get dark current taken at ``t_ccd`` scaled to reference temperature ``t_ccd_ref`` - This uses the more modern approach to dark current scaling which does - per-pixel scaling instead of using dark_temp_scale(). + This scales the dark current based on an exponential scaling factor that depends on + the dark current value of each pixel. This is a more accurate way to scale dark + current images than using a global scaling factor as in dark_temp_scale(). Parameters ---------- - img : ndarray - Dark current image - t_ccd : int or float - CCD temperature of the image - t_ref : int or float - Get the image scaled to this temperature + img : float, ArrayLike + Dark current image or value in e-/sec + t_ccd : float + CCD temperature (degC) of the input image + t_ccd_ref : float + CCD temperature (degC) of the scaled output image Returns ------- - ndarray - Dark current image scaled to ``t_ref`` + float, np.ndarray + Dark current image scaled to ``t_ccd_ref`` """ # Confirm t_ccd and t_ref are just floats t_ccd = float(t_ccd) - t_ref = float(t_ref) + t_ccd_ref = float(t_ccd_ref) - # this comes from the simple fit to DC averages, with fixed T_CCD=265.15 - def get_dc_exponent(dc): + # this comes from the simple fit to DC averages, with fixed T_CCD=265.15 (-8.0 C) + def get_dc_exponent(dc_in): t = 265.15 - dc, t = np.broadcast_arrays(dc, t) - shape = dc.shape - t = np.atleast_1d(t) - dc = np.atleast_1d(dc).copy() + dc = np.atleast_1d(dc_in).copy() dc[np.isnan(dc)] = 20 dc[(dc < 20)] = 20 dc[(dc > 1e4)] = 1e4 @@ -389,6 +388,8 @@ def get_dc_exponent(dc): 1.90718453e-01, ] y = log_dc - e * t - return (a + b * t + c * y + d * y**2).reshape(shape) + out = a + b * t + c * y + d * y**2 - return img * np.exp(get_dc_exponent(img) * (t_ref - t_ccd)) + return out[0] if np.shape(dc_in) == () else out + + return img * np.exp(get_dc_exponent(img) * (t_ccd_ref - t_ccd)) diff --git a/chandra_aca/tests/test_dark_model.py b/chandra_aca/tests/test_dark_model.py index 5467b1b..d0e0f7f 100644 --- a/chandra_aca/tests/test_dark_model.py +++ b/chandra_aca/tests/test_dark_model.py @@ -8,7 +8,7 @@ from mica.archive.aca_dark import get_dark_cal_props from mica.common import MICA_ARCHIVE -from chandra_aca.dark_model import get_img_scaled +from chandra_aca.dark_model import dark_temp_scale_img from ..dark_model import dark_temp_scale, get_warm_fracs, synthetic_dark_image @@ -84,13 +84,18 @@ def test_get_warm_fracs_2017185(): [ (np.array([100, 1000, 500]), -10.0, -5.0, np.array([171.47, 1668.62, 852.79])), (np.array([100, 1000, 500]), -10.0, -15.0, np.array([58.31, 599.29, 293.15])), - (np.array([100, 1000, 500]), -15.0, -15.0, np.array([100, 1000, 500])), - (np.array([200, 2000, 600]), -6.0, 2.0, np.array([478.16, 4299.02, 1399.40])), + ([[100, 1000], [500, 600]], -15.0, -15.0, [[100, 1000], [500, 600]]), + ([200, 2000, 600], -6.0, 2.0, np.array([478.16, 4299.02, 1399.40])), + (2000, -6, 2, 4299.02) ], ) def test_get_img_scaled(img, t_ccd, t_ref, expected): - scaled_img = get_img_scaled(img, t_ccd, t_ref) + scaled_img = dark_temp_scale_img(img, t_ccd, t_ref) assert np.allclose(scaled_img, expected, atol=0.1, rtol=0) + if np.shape(img): + assert isinstance(scaled_img, np.ndarray) + else: + assert isinstance(scaled_img, float) @pytest.mark.skipif("not HAS_MICA") @@ -102,5 +107,5 @@ def test_get_img_scaled_real_dc(): dc = get_dark_cal_props("2024:001", include_image=True) # just use 100 square pixels instead of 1024x1024 img = dc["image"][100:200, 100:200] - scaled_img = get_img_scaled(img, dc["t_ccd"], -3.0) + scaled_img = dark_temp_scale_img(img, dc["t_ccd"], -3.0) assert np.allclose(scaled_img, test_data["scaled_img"], atol=0.1, rtol=0) From d22b7fa42ac6f4795fe7130541279934039427a9 Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Tue, 27 Aug 2024 09:48:53 -0400 Subject: [PATCH 4/8] Fix black format --- chandra_aca/tests/test_dark_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/chandra_aca/tests/test_dark_model.py b/chandra_aca/tests/test_dark_model.py index d0e0f7f..7916c29 100644 --- a/chandra_aca/tests/test_dark_model.py +++ b/chandra_aca/tests/test_dark_model.py @@ -86,7 +86,7 @@ def test_get_warm_fracs_2017185(): (np.array([100, 1000, 500]), -10.0, -15.0, np.array([58.31, 599.29, 293.15])), ([[100, 1000], [500, 600]], -15.0, -15.0, [[100, 1000], [500, 600]]), ([200, 2000, 600], -6.0, 2.0, np.array([478.16, 4299.02, 1399.40])), - (2000, -6, 2, 4299.02) + (2000, -6, 2, 4299.02), ], ) def test_get_img_scaled(img, t_ccd, t_ref, expected): From 11cd253f7d8ab89e0d87f4e9bd891bcc88abd8ba Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Wed, 28 Aug 2024 05:32:35 -0400 Subject: [PATCH 5/8] Streamline code --- chandra_aca/dark_model.py | 32 ++++++++++++-------------------- 1 file changed, 12 insertions(+), 20 deletions(-) diff --git a/chandra_aca/dark_model.py b/chandra_aca/dark_model.py index d247b94..f8f591e 100644 --- a/chandra_aca/dark_model.py +++ b/chandra_aca/dark_model.py @@ -373,23 +373,15 @@ def dark_temp_scale_img(img: float | npt.ArrayLike, t_ccd: float, t_ccd_ref: flo t_ccd_ref = float(t_ccd_ref) # this comes from the simple fit to DC averages, with fixed T_CCD=265.15 (-8.0 C) - def get_dc_exponent(dc_in): - t = 265.15 - dc = np.atleast_1d(dc_in).copy() - dc[np.isnan(dc)] = 20 - dc[(dc < 20)] = 20 - dc[(dc > 1e4)] = 1e4 - log_dc = np.log(dc) - a, b, c, d, e = [ - -4.88802057e00, - -1.66791619e-04, - -2.22596103e-01, - -2.45720364e-03, - 1.90718453e-01, - ] - y = log_dc - e * t - out = a + b * t + c * y + d * y**2 - - return out[0] if np.shape(dc_in) == () else out - - return img * np.exp(get_dc_exponent(img) * (t_ccd_ref - t_ccd)) + a, b, c, d, e = [ + -4.88802057e00, + -1.66791619e-04, + -2.22596103e-01, + -2.45720364e-03, + 1.90718453e-01, + ] + t = 265.15 + y = np.log(np.clip(img, 20, 1e4)) - e * t + scale = a + b * t + c * y + d * y**2 + + return img * np.exp(scale * (t_ccd_ref - t_ccd)) From 5d0adc6a505dc251c440add4b5d6dfd90394ed08 Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Wed, 28 Aug 2024 13:11:42 -0400 Subject: [PATCH 6/8] Add reference to modern notebook --- chandra_aca/dark_model.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/chandra_aca/dark_model.py b/chandra_aca/dark_model.py index f8f591e..9984d8a 100644 --- a/chandra_aca/dark_model.py +++ b/chandra_aca/dark_model.py @@ -12,7 +12,10 @@ Alternatively: http://nbviewer.ipython.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model.ipynb -The get_img_scaled method in this file uses a more recent approach with per-pixel scaling. +The dark_temp_scale_img method in this file uses a more recent 2023 approach with per-pixel scaling. + +https://nbviewer.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model-2023.ipynb + """ import warnings @@ -353,6 +356,9 @@ def dark_temp_scale_img(img: float | npt.ArrayLike, t_ccd: float, t_ccd_ref: flo the dark current value of each pixel. This is a more accurate way to scale dark current images than using a global scaling factor as in dark_temp_scale(). + See the reference notebook at: + https://nbviewer.org/url/asc.harvard.edu/mta/ASPECT/analysis/dark_current_model/dark_model-2023.ipynb + Parameters ---------- img : float, ArrayLike From 84eb5e191c47e69fb02f476695ada65c52d4b67e Mon Sep 17 00:00:00 2001 From: Jean Connelly Date: Thu, 29 Aug 2024 07:17:32 -0400 Subject: [PATCH 7/8] Update test names for function name --- chandra_aca/tests/test_dark_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/chandra_aca/tests/test_dark_model.py b/chandra_aca/tests/test_dark_model.py index 7916c29..2d2e09b 100644 --- a/chandra_aca/tests/test_dark_model.py +++ b/chandra_aca/tests/test_dark_model.py @@ -89,7 +89,7 @@ def test_get_warm_fracs_2017185(): (2000, -6, 2, 4299.02), ], ) -def test_get_img_scaled(img, t_ccd, t_ref, expected): +def test_dark_temp_scale_img(img, t_ccd, t_ref, expected): scaled_img = dark_temp_scale_img(img, t_ccd, t_ref) assert np.allclose(scaled_img, expected, atol=0.1, rtol=0) if np.shape(img): @@ -99,7 +99,7 @@ def test_get_img_scaled(img, t_ccd, t_ref, expected): @pytest.mark.skipif("not HAS_MICA") -def test_get_img_scaled_real_dc(): +def test_dark_temp_scale_img_real_dc(): test_data = {} with open((Path(__file__).parent / "data" / "dark_scaled_img.pkl"), "rb") as f: test_data.update(pickle.load(f)) From 2b3c94c1d649b1cdcacd3a70400c2e689bcaeeba Mon Sep 17 00:00:00 2001 From: Tom Aldcroft Date: Sat, 31 Aug 2024 09:35:13 -0400 Subject: [PATCH 8/8] Add test of NaN values --- chandra_aca/tests/test_dark_model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/chandra_aca/tests/test_dark_model.py b/chandra_aca/tests/test_dark_model.py index 2d2e09b..c9a4d95 100644 --- a/chandra_aca/tests/test_dark_model.py +++ b/chandra_aca/tests/test_dark_model.py @@ -87,11 +87,13 @@ def test_get_warm_fracs_2017185(): ([[100, 1000], [500, 600]], -15.0, -15.0, [[100, 1000], [500, 600]]), ([200, 2000, 600], -6.0, 2.0, np.array([478.16, 4299.02, 1399.40])), (2000, -6, 2, 4299.02), + ([2000, np.nan], -6, 2, [4299.02, np.nan]), + (np.nan, -6, 2, np.nan), ], ) def test_dark_temp_scale_img(img, t_ccd, t_ref, expected): scaled_img = dark_temp_scale_img(img, t_ccd, t_ref) - assert np.allclose(scaled_img, expected, atol=0.1, rtol=0) + assert np.allclose(scaled_img, expected, atol=0.1, rtol=0, equal_nan=True) if np.shape(img): assert isinstance(scaled_img, np.ndarray) else: