From 659b030e3cc78475a9125f75ec27e1422844197b Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 29 May 2026 13:52:25 +0100 Subject: [PATCH] feat: ingrain extra-galaxies noise scaling into imaging examples Make the __Extra Galaxies Noise Scaling__ step a visible, runnable part of the core imaging workflow scripts, so users cannot speed-run to modeling without handling contaminating nearby galaxies. - imaging/simulator.py: add a faint light-only extra galaxy to the `simple` dataset and write a `mask_extra_galaxies.fits` covering it. Also normalize the file's line endings to LF (it was CRLF, unlike the rest of the workspace). - imaging/modeling.py, likelihood_function.py: add the __Extra Galaxies Noise Scaling__ section (load mask, apply_noise_scaling) before the mask. - fit.py is intentionally left unchanged: it deliberately fits both galaxies of the `sersic_x2` dataset as the science target, so there is no contaminant to scale and a noise-scaling step there would be misleading. - Regenerate the `simple` dataset. Verified end-to-end with PYAUTO_TEST_MODE=1 (256 px scaled). Co-Authored-By: Claude Opus 4.8 (1M context) --- dataset/imaging/simple/data.fits | Bin 5760 -> 83520 bytes dataset/imaging/simple/galaxies.json | 45 +- .../imaging/simple/mask_extra_galaxies.fits | Bin 0 -> 83520 bytes dataset/imaging/simple/noise_map.fits | Bin 5760 -> 83520 bytes dataset/imaging/simple/psf.fits | Bin 5760 -> 5760 bytes scripts/imaging/likelihood_function.py | 28 + scripts/imaging/modeling.py | 32 +- scripts/imaging/simulator.py | 499 ++++++++++-------- 8 files changed, 367 insertions(+), 237 deletions(-) create mode 100644 dataset/imaging/simple/mask_extra_galaxies.fits diff --git a/dataset/imaging/simple/data.fits b/dataset/imaging/simple/data.fits index 850afc2d4b7092f41b14e217317eac8937ed56e9..42951337fa4f51d1b75992aea327d9ae04f68cb1 100644 GIT binary patch literal 83520 zcmeI*d#tC&c^`J(sZh1NN1;Qe2$5!AT z=wgXQ7#<>}gH6%kG))f8$fY(CtawqUKB=2SuR<|v^N93{3{6qylFN%oky>7sT;A_j zky6iZ=KFb$^&Sa`ivanD=bzb`cjlRAp6B~qX6Ai=`%CZtA9n5d+h?8izO%l^-??YK z>#V=~iI0Eu6QBCnr_Z|J6KDPR??3mK&ieEPANlyNe&kcX`aSw@H2t^lKX=#r&;OqE ze(JqHvqA48pZe5C{=r$l`jHDha@Ga^;FG`hy{T-r?)>@l-~US+)$hsQ-`Jq{@&EZ_ zzw&FJI_nc3Mful1_G=&i^sj&7AmNB(tFRl-~H!0{zv|^)B79$>F8|; z@0b3&pa0-TJ|_kaBHAOGe5NZ#h%2EFtDlhHf-7vKL2@89`-_S5bl)%h!RZ&SkMiR7g@yUwwzky!*SPk7Uq1i)KimKHt{Z;kwa>h}p8xiP zYya@?obh1I7dq>EGaj7@-I@=d|LQBRyz|4CUf1>ve$ZbZ^-Iv=3n6AfjWlO&;@{7D`-pr3)tp~r>VO``yp4zAN z=&|_7x=X+0n6i0}?@7JZMPBrm^?0saA9PMP-?Y8tOCRKbAN)_&Rpao9A97k8xNiG8 zZZ4+h!h?C>OF#3D^1_A={;Vg@OHZCp_+VG;0KJj}K8#E0d5zOA^ILz>$3Hp~{)uMueH?6i<0sV7H=dsVAao>pmb0@Lidal3u+2dIEuEI(9Am?@VK+fvu;2WPQOP^V1)zKOAZpSD7OE3H@ zoZ*wa&4(X!@kL&l2mPsD3vb3v`)fCkGVO3t`qtn4_zZv6ho2fh^nI~*eDY8H3cXbO znCDote&ycpj$2#gKf%`)oQIBaH7~tQ^-M1CIOK48;c9XCr&?d;H*dv@^iwy#lAYAJ ziG0iN!wDSeXWmi2q^)+oc9guW&p3J>^xj(f8+`6IuJ|K|*!z^vvM2b^pWSLtow_-*f-)s`3ZP6E@kzHlL@Zy)yOl?*N6VEt91@JRla2%p4U8)1G-z< z$s2#>(Qcji#25PrAMn%E$xUh81po5B(zEq0$IkYq44twU{IUo0k~{gX28Z}ZPd_-( zf6TkRktRg8}rU;{A3;G z+ma_v*Ae}1(a+%jShdGLsn@*LiSHrji}~k8M?LZj8|{LoAKIf-o?(} z?x}pLKE_`{E>(_kGjYQ_=vI4h1Ggnd{+gVW6MjnHQ#nuV(>k?_CvXxz${$59_#w}# zXWjIwB)_6F^t7G4$Qd0a`ZYe|&|Qar>&N#L|LoelO6yk|H^tkK%X{!?e(S0F_p{r~ zPoLGFJ%|f~uU)lXeBdu7xuBze(GMNZ^(WWjdoegBA9kdjzRgSCqkf5b(c_QHFOcVi z|I&NpReH2uehe9Q;BmPZ&Nv=6~BtPu3lJ z^f0y8iaYxAqp{E8H|^GCT=64*uxxeB*y+q|Rw68PI9 zkMPw}mVPF5Yh6=4;75P>FmL$bcgfv;qWHqc)Sef++-Ai$`J^nm_#g9~Uhye>6#w|N z4&%l=yUPyDOYZxIe=0io4L$Z~UU^6LM-Tsd2Oh}1`h~B;CHdhGKgO$vkD8yq8uMLT zxDDSmAAZa$o}^y=(BI$Xa8c=FZDIab=kA4{=#d|*`OU|!rg9Nq@FTA5``-JY(wFBe zgP*Bg7S7CzEfnFXjo|l6&gaKm676=%3t`=AoyO z2R{=!`WLFN7PoMaTS0v2*x9U!5NGQyNcS>c(XqTjXecSDyC1 zPyX&5r4|=zQr!b+An$z zZ>4wR2cO5PAGzjvOWpJGQ{<5ODo&A`GB~VtqC1rr`Ap}>=a^?R?(7X-B9~e}`A_*9 z{9J53ruG$lR=#Ha#;co0Io366@^*5O*I1|ZprfC5{qU7`>nr`C z8@>u(|1S5wYYX$gZE#rjOON!RU*;`3@X3DhwG@3>x98gR!@qvyL%vG-AdieMy=hNB za*%IKao|2lX&+IzU73qzVOy;eANbF6xI(}9(GOo4Uv^`?nGc=JS9ru1KJm3Lc?`PhQ@Zx?aFhLxct{WY zZQGx`4BaVz=(9)k@bCGQPjxs@;`@rc|3I%(zokDq)eicbebC=ka?1SdMM-|@6;DD39_Y2~1poAIzO>i;6MVPDtP?*OSL+?mx8zOyulcvw zuW{s&dGTHMN9G^%Zijni^`~FBfJgEj^e?IYtFb5a;o++O{>^96uYCZ0@eA;XFZQ#x zF#p?Db_(bCF^~KU&OFDz=j>y;zAv6P|2wyt4?p3f^l7|#@TtU4%{ z`w|z~E&9fr7r*4=IXuCQ`N&)Qw5~lRf3`m89g{CiaiyIraRX@-9GdN7Y@=JrinRC@c6hOVc!NXz#3}g<{;dZeaPv$R9n=&8Oc|*W0d&C)$nAJf-i&-hcCa7CrN1T=ftAvKREtOaEEN z`k;4R*`I#+N&hXl#D{+5KCP1jKJ-)S@45A69KM=<)=i&qfFC8iSP%Wd)0D4$@sDt8 zy!rMg|JTnv=$QvU@VmA!|Jzpkht=4vcKx)YXMQEznP1(y@Qn_6ng<=@ugv?G_AA!s zTnkRD(>%tTXHS>Qh2_^gH$VBUY)D~44uzxCD{o%za@k&a)D(B*8+zy{`PuUO>CXkPkTM`^cjXy#9ejznc|Zs|`R_%c7ci9G42@T?zs zlb?FK*dKl?pCm{AQ$KR`Tpcd-rw=$JXYK0fz|GX( z>Th4H9~@i1`hnPue%JK(L>_DC;uHSSB{y`~h5GuX{eY=(r=4BEKYsC{Za(v{3(v`E zsq??9_W$gNT#dIrb}_{*xsda6r(5@F_HBg={pmwIQ#XG~{uI6|p6nU+sJ-Y<>A+L^ z!*$VN2baaZ$t8TPkLPbSuJDhZ`OHIJ--llLpFA&rL4Hbp13ma5ANnCrdZu6cyk?>Q z)x88g!V`VL8Jv*ImB~x^Q}nf?2gjaI@$dOH^RLj>{$aa%d?(WQf4Ig6{;YegFJoJc z(+?f=(fe}V3$q?}fsXZ%Cw|G5{OH$wD4&;Udj6O%ecZv@Wmd= zKV;s*fApyzzq|#vaBhBhHqT|<9xg7t;yd%z`p}v3k#&eaA4b~|2W@I zi8p_I&;0M+hW|h5?+t#T9qjsBs`Dpyf=~QS?NuH~PWs`KJg$siD0$~SN%A9Kbb_lb zc`ALQt4`18Ehny^tKXjZZSPMc|W8H~2)ql0GV@KJqbxh+Wdm<+ za-qNIr~I&X>rs*mzX4DD$J)SU?X&3Db9#kSdNUtfps!!_T=-e+`kgg;Wk2Ri9sX*c zu{7+dDbXGMH@`nIFS>9zoi}y+qKXgtSs!^@7d&SkeAam&xhl~y9?sFlANkU^d8Xsg zrB~0_7UqB3%CDopI_$U66TI`!{JD4m$8f*D+rdSRzyI>a%T4e5*H>doz7-{${7^ZaW)54)-zJa*;=#}Bu!_I>|fZ-3g)13ulZ|NH%~nIC`V z(|%2t!$SG@>AJ|7JUl0NCA-ugzWEP0)UKo-_C}B9TkiVU-u8@mQTRzazg~En^sju% z{GPA$b=P~>@I0ko`On18YF~lggdg*YPZiJM+xo3{HMqkkJTCV3P5Xy`azjT+-qxp1 zFUDynclw+1$G>_`F9$kbI}7LV319F&?FT<_%pUmjKb`+v#zptFX~{kgCeE-6@#4k< z{jc40PycJToj3pcK)bb@)6cX+(f6NSIu8`TRd#a2<^8YSu<`2{UwP%H8v7Q1%!glm z$`9GQ^|E{Vuzvd0u1+uX!LQQ;eXHxYXV^pG%sBm*V>jxVr|bm|_V#|8@1KM&d#ijN z9e51B#HorqW8PiPi^`siUyB~hgU|3^{jG<5$X9=U$2!#MgP$j_rL2>^(E~d}2Vdym z7a#CzoVcJa?_e+V&)>2S`OLM+Ti{T?Dc<34I_~qQcH7(ffOgybXJNkbYP+kyhWm5vu8ohEKK8S3YKOAlzjYqIe)sl@FE{ml z%!!qsu@`ZWzaS^;*A54(T^~Ct{u)oup3@ioisSI1-xs^RpIhhCp#z`ghgTvGWpGe- z1YhuS*|4*=*!i7RJ_Rp(2dfp0;^*Lt zzL)!ao9{Ql4IBkmg=_uU(V%;Y{R(?1`$H%5vzv-D_za%QewT)wRX%3k^b>E)Q}Y^U z9sD?X@*Ct}e*MU6Drb5UC)kftf8*r^>hOTh)!iL`sm@iu(*4C(8jNtZHT%Z?!~e8jSpYC5utXsBiWZ9?eXH3_Qb}m;QBlM zX?tp(t^LCP+n&z;)KKBG7tYjAN-KGpp@^^4}H?VlHToO)#X87 z>G{l`v_I|n%BjtMKpxNU-q`s)yIq?+i(RS1|IPig?T=R;1)n$H)i>0wr$HOvDczp$ zG{lMfg7;&+-`(x$;QGb+&jarK&Gy2^ulc>}hQHfh>h^b9dp`NsGr|4QZb$!nyYKS( z-)Gz59`7z{hjKn=z2YCf`HusAeLKr3M$=*KCt~$r?oc`%G z_%1xc({dlbz1oAz^6&5ydh3JE>F`?k4_?`4jmMw))!D<+;2&M_sq7P<)&=KEa`0T9 zWE{DZYt~oixZ2qVJzgEZPkz@d^uHR%pX+Db_3^Lz%OBwWVDcvRqaEJa$)C+%?^u6F zL+@+-n(yzluWq3BS6=%Uja{&Z+Gq6b#6Pcny8XnRXU+xvqwO0#9{gl`Fnl~U4 z8S$DO&;vcuH#?$7xIhoiJQo*~=y`6O=dlm?EWaMQ+RIOyZ!)e0ch!%*uxsPk$>)PV zrS*yr`YjFnkzZB5g+712H~q}Fzw7U!+TW~ocyGR6t;L z&Wn^cc7D$;yvvvD6O=dS9G5*fSAsw1O1H!=?89%J|7_xk|5NP~zTV@-JKHzsKlA#o zKW#_vng89}o|?Q!Aq*(vAD1Hv!PPKbTJGXHA!eV=Sk z_x$X`ZJ5s$4gJG^tN)e#-jsM_|6$+8PRLKb#GjC(cKH=O(;t25C*HFYc%Wyv(XP&~ z!l8NC%~J3cycW(Ed;6@d1Yhb(_+IY)wwM2(+RKD~OI;kQ^96k3D{-LsnLf8Z`vP*y z{B{1WonGL5>Zh)b-O?L<(wlMi4fF!{%4@p37Mk-@__j}wcku75tQ!?E+*&YS;z zpxxB*`5VoC!uiwJM+OeLO{mpi7r+0eUg*a#Z&augpeEAu1{y@jk z&blWeXZ_j97rQ*pt#}R>^r^nu>Fp^0%Ac|WxYZ6v+V}SHo8PsM_FYq)!?Du7K-?&N zuXJ2&euuo+>Cb8$++No8e&LpTP|{0a;4jEmz6Jl{GJVh+J@Nla z`lU~K1G_oc`CO`VB{*geQ@^jBzlUf2|L8Awem38Chd21y z*X4b2i32Bw|o13`yaIDhrpRn)p z6Mg)@sQZ^2M*dd*1n$LMdf?B*1tmSgkNtwY15V%!jy#8V_5<&5=Q%v*dF{6+h~Wtmoh7+{FG&yME+izVb`S6O7YNf8qfDMvv@= zeRzIjm&4gj9FSMYlf()3aA&`7_;Q`M-kbdnyo(R~@i%hbba(FQl|un}_Z1&zt{! zpxO7ZPwV6l*vIvKdvmnj|Kd-k`4_)+<=FT){*0dB7LM7AIH6t19^en2v>V61_|@gE z*X?bwmoxi-!m0M9*nu+i^;_SY=SuU$9&6vhp3SSjI=;-a+Uf47?`7Dx^;qA5-fw5U zkIkQG7ca$Y`#W`U!a1w_m_FIVP1%=D&z0B#JizyD!H;(DXNVKGcfbA7I&Za)zPIN) zJK8-&|kr(C{DKEAoIo~=mG*bHvc6a2=xqyDO@~->3e_Ch{mk)Jwu;76`Yeyem| zP9ImtZ?Fq_Ex*W*!P^(Ryw8QpCjWvLxPU8m#%|z5JO7+|;TS$YKll?jw(Q5@8lAnd z6a1l9t^$VuCHj~$^-q}>$$4qOa8WTIK1(9Z<7z+nYef;@kbnSj&XDBn7xoEeZl)w zkNm3rx;%mY`AzNmiyv?ZH*lrY51rtx?(aN@Yj|Ajc$<0e+*hDuJNd+K)OYy&hUeD?NAg(t5<3w$;GUiEck%#v68|TUx+V1CL;R4B zAL;qr#`$OB0za&szrH{D&HV$w4bJ)RhcBH^^jqzLZjWDRN3$P4c22*as5k=e4~+BK zcC_y&exubqdB1AfkNt=%_vJlA_w4dBaZtQtAI{N~?16uAUd{hZ<1+un4(LTZVMpu- zj+AhsE`OF!EyaJY#NWaj|Ei`csNnNZ(zdfzAL#p*e_Zb=9_jC|e5~>_cJ6&&en37Y-#nB!;(cHD zJI=NFBlksL>vG!o?j+}i{M;WWPN~!DTIjPYdcQjIW*_jjKk=RYO#LrB^1sR7x17H& zj<`{I5Icd(f7Nkva7%oEbNrZ3XpfTP_S5Xce)yZ+&tKM#^apph`-59{@KE1({JqA%=YK+K zkH)V*l6^M2cqD%Ov7hdLO}ua~_eAjhNZ!Zw-p6Cfj~-9{#J=2nJhZVZIdkVfs<^@~ zy|??#5A}g{KJ5J1IkEkba{_rJzd=soh<%%URQ#Y1afRP<9x{#p_T%v2Ih?Q`c;iRa z;Z1+I(=Yh;of4cY+fw2XTUoyI2b!_8;USezFVg;-M#)~vkzt;@^tx){Dpn+=j?%BeZ}3Uh8k4&0g%U z4<(PSc#(5Yal!kj@@#SAzKpZ4xIg}1+~9ZN|M4Cu9%)a<&YpQ&2YQ41XObs9m-nTf z&wGB)&!fycKir->qXV%$pYhN2=1bdi{XXJ`_G0eso(o>RhyHBhrSCY{XXR@HG0lAk zd$09%y)2X;IFj>*l27jQ#noxNo$gDVd$14surHa$6?US({7T-yuTJBM_`y!#2d>xy zKMs$3dc3@_-t!PY#0B^k7tk@^G@eZ3i1Fq_|BB!TKlnC}e(Z_8SU>+Ej##IBl^yXf z>_fZ|zn!P>V`Ufi8QE`tJ@329Bg6}S)B8+Elh@c!Jd%C+W8F`j-5w9lp9+qj$+^?B zkq5lLko%1n6CYmAzTxHI_vPeKuO?r6G5h%!L+8bQzHmx?w}I|6eP8kmm2ctW@y~R_ zsl3bn!@k5mWlXcZ*|#6b`GEb0eTuljKI})FXO%ySAH1pONgLnehJSHHz6Jl{zIXu# z2fBQB*6&F0t8lay{{^>j(4|s=X<5D`mt$j)CrQ(Zt#7-vTE3d*A zKH(eQ*@<<@o5;m}f_${2<9w>lMY^9or{3fHTF2d0?Z)6%{J%B#7xtC(YCk&CoD1sw zhF@hD@)7%E`|Kx^*Vq?46MOUh1bn}oJnE&s4|qp=Eq3sF?C$lO$NAK2$*;bZ^QG7F zz0ymmzmk40_V+|yYs!if)35s9k=jM8`hGddY)d;b*jyGkM-{dH-w95BV!`oW5@FZuH&tev|Wj z`j>Y*hkD?Qc>?pp%{i2NiDSuoo<6e+pgj@)_;iQ&dz*ae#q;}yqP^1Py3k$?p1*Zz zccio9d*1b9v*UX_yJL3z=r=pA?;C!9c6{Hre{y!*pgTLR)B9k1wVwySuf5vO%l=t= zIrihdLhs$OmnU`BSq#2fwa1NZC( z?%6H;TZewi107yF>-_-ngx?@H^0;;qpKt8r&#vzz+}Da1^d#QP3zXi|mPgdRUFLD$ z#?LB!FZ5L6z}QBd-=6R1QETmmVUT4HuZ?~8?6_f9vl9l+XD5#H{@F>x@y<>jab$MV z$opm|p4b0oCw#6O;Ow|z53}Q9XW#Di|D){h<=CgV^?dgG-Z%9AuKm$t*?-l3Bzc^8 z;(U&L+~;`zQG5_jyjSU*lDA{wtq4n0TcwZ&tz&{lf|Sk(a6Br^~0{>FL~SIfvq>pO5_PBVPM(4;1a& z9d}23mp40U;CpuJ*r&`+>-o=@W^W$;V)o{-Kbf66&aGyrj_>|vCyw`|X2%alKRf2`QeosZckJ(l}<@$lg>Pn~Cdv$MZ(zth*X@m@gg zQAW5st=@094`C1dk$7@b;w!yQ^8omkS9|Wfg`MaJXK(~R!GGoX@CV=1x^{5@?&Uve zuXC;RuX{xJhxaKT;u1c^Kke|&{?+0B>Vb>;JsbU$_8I(-eTefzd93eEy>IS5+_=Mu z)4tE}{*&)d`A7LS`(O|56`mY^yzU{MkDb^zzncC2>tp}nJG|Km!T%|V3#aw(i=LYO zVB*6<--mpB_LtA-`P^)w!}^bAZyA~I?6izKb>M7v@FC@yz&{D0v& zeorZ$O#LTZzz!+Wb9uaY3A5N7}_J^w=G{z!yI7 z#eT#meqw)c>mK9U=<%xPfqjS{;s!erKkU=(yX9Bn{2jyJHs>$m2D=b9?u$N;W#246 zhbQ+8?&+RCvkS1jnEk-F`nC5@&yE}Cud@@!dBE(H(Qo$Vf!o=(aqczyp~Qo~n)1wk z?)$H2XCz+y-HJmud*+{uI~rmRqvUI7Y9Zh zt@pI~Zg)`^DK0j(d;UUrAi} zv7Sf$yV={v`O)m{|8(BrEwi)!up8R!hhql|*{8gDaqg!!-(!h9#0~!AzS!0D{Sv>oXx|}z4DtL<{jQPolpDI= zyS(vx^l6_-|L_7I@Fbqt=knJpoxg4M-YVR~Cw~vW@W%h^2fuJ{oc`vE{nWYc-flOW zzbmBv+tGu#fsXaz6J6_-SHZt~4SDAu^|bmYHsu(<#x9inll)fxOJ85_`uJ>ppXxnt z`IY_7aD z4?FSu3{Q7?UQ~7=zjD4iji=tr_r9imIXhq%?r8_wX7ABW4u_<#1mu9TlozAt{T6VKsZiOyw-Pv*f7{^f=E*Y7~*XJ`G+f_CwN zT*SdIXMa15FU~vUSLE$IIPqP{E^g0x?w#j!|6S*n{3E-t4}U2ACXZ<)Pw#I4!eyTE^F z_9LH}GxXcDg^rUiRDO0c`xs@bbF3a0?rAUP+}Am8*+tHQpU8gMxr6<={L*^^_UGHQwB{kQ)7G=I-7*crQld-cSV ziUV*SJE`x?(dF;i2YZ8i^w|fz8)shqJQo+(0lU)9pE#!wSIEgZx_#gEU5-a9e)DJa zBz`FGj-7t9?`Qv7y=VJC&Nrv;f6F@_OFR)jydUDdke9}OqrEyDY0K}c)%nsXr+3FP zJ8k4GvxVO8@a&Alga7ge=Yi#)&))IleZw((`@wzzKl`!4*X*w(UTlk<{H5$$P8{p; zyCdzj-0!`beEx;p+dJQNPxEx{=O0h}6E}QMJK%XmeUHLUzR~CZO#N=+^xngH$+Wa@ zv_BN*#SJ)KzqD^W>;GBgAK?f7;Kn%ip)UVYvJd!#b2t?*E{p${Z!LBC?yBD@S?=}i zb&hMixPq>7t-HN#^}ad&*o}Ft55N3@bsbE+B_DQj!^j(&yi}a;=x@&vSkrj^%H`d)*w{q+3{naQulk~b*WotZyxV`%-)jy_(IO9-a7KS*+OrBWcGvdsq8OK z9qXujljFxc&I6j?Ieb3vr93m{uiu+~bmVjOzP)p+qdCWVFz=1I5A`0VcuTI{bB2GV zxWOKrPu%3@p|m{<=BJwvgCd8s3&@! z_20C|M*OMYC#kq`#_%7`?c9%huiy8>@>csh_kaAsJuitTF2k3|k#-W36e6sH^cKMuJ=g{h& zvlDi6ZSFbv1@TwA{Pd>mBjvB+3V$Vji2L@B^y$0MJCkP}$@yjNH~N0|-sXJ@{`6Sz z{#3p*^gfCGx%&iov`>Gf`{Qfd>kIv_eaC(GsLxJJ{5U!L^^>y?@w>O$PfT2iowe^| zpY-kIBkni6$MRhAEBT*$OmXnhe*f`v?ctHH)%UgR!o5(H*{6xa-n$Yn_#x-(-jAjy z=P2R>zbC(9ANCpSfPKIZ+`u{fdCorI7rx+^U)IiU*vFocKh^#muEh!S>OZB+KF~Lx ze#UvO6xZZU_9xn{)BN%+@{?zg_l=3u^4&unZ`+#lL~;MN^X7jaXz~I3OY!1J?pwr> z2ebcqu>1G7Y~GippYQOV>ifsvYL5>;-+XVuK3?F@lZU;M{nM+7Bi<|gPU439!{cKg z$LC)B_?&kgcV=I3V-K(9-sPnp2M#p(o$o(*WX3$)w~~=_j}t*IY0BA z2s^ajHvdx}njc90a@ntY8|yjR@%DlG-7EVG@q=9CQ{x!mg4%zJ7vi;ZMENcKIIm(K z_SJM*$vcFE4Imz?BX2fDrOY}W)|Aa3sI7kA3{Z^8m7Y^F8*pzR_&Yp5B4o_;KIwz4*Cqh|PB%Rff-3 z#yM?s@BVV=z0l9GUMYL?y`b|l-$6Myz!y7q-|xKPp)Q!i_4^R+H}38BzN795 zoI8?-@h{}b9=zu!k2(~*sZagg)!BC{#TE7hrK+S!49iF)u^?-Q~MCH%sF z?4rJtiapftAj0!X)(7Y4u>xe{`gDwK%V@McJZ2js<@tW5OH6g z<$j7C&@;Qar~lnszfa{nv%Y_b9Xy!)U%Fuw1o??s+X{PW!J;XWPv@!h#|TJME{(e zp4bWda9{G+8RPuN_Y_Zby#GdXA6MU<=UhO1d3L;K(w-gPLACmQ=@rSl2jiQN(W$&>gQ-*F7H-MAl*KAbmQpY@C9 z;sQS={_}701A6BFmF!_Hehlv61@4x+Jhs>W52hd7PU1#=zqDr(AK0T3UHA>&>t39F z><@m}A6$EmU-q!C_rJJ_3+S^S^vH$W)SXwb7yd@cAF&($%Kn7>onw71`vm)T`w;rR zV+4XNzauSv+|%uA<6JA}y2^XAj}=e&=?8Njbu|0w+Mo3Myq{`i4>xy1toMQK51fyc zeT9GR?qTZwGWXo~Z|s6L-q&i{@yAZwQ-~Y#Nc(L5!9GGf5I6Xl>k_}5%fde#uoL>S zU*zBDS3Us$;sktXXJ_yQM{sQ40*~@{cx5MW$u43?b>76T_;LM%cX6ZgG0)8d=j_%x zw4;l^%VMAUk)M9n#qV1`|1q@}`6{{Em(UYA-WdPNui1Clzsn2cS+|dS>;@n5Ece9l zQ}^t_IX`_Q`4qn_53ApOjvqhL6PLT2xMLr1blh{)y&U^;E>O>VeBJoHw$Qu({OqSg zmpwS2mN!`sd$_yHXJ^Gv?`huJ>HbXJ8Z}w^AtX-X7VjujH^Ii6$E>CbzbxZOB-%E)z&MCzgzf1G=%jbWet^dF5 zuAXS!U+2Q`=KB%f;oASIAAQ?A5d2t!>j!d=Gs4#u^}AB`JMQP)=j+c7{eFb|O#A7w z>u$&ItaDoHyK{q!jqkA<|KR-`af7`08S#Kz4=v39w$<+%IuDUo(Jy_e!~3)jchizx z$U{~#FC6ZP|Ao`!8)YB*FLwOTs$GB2;XU?L{XEB4>_;4_`;oPU`QNtse?(VO7au~m z?laBDKJ;ft>;&%F2fsp&;sp8gNBorh%KL8gXkRGLn(7@M4#ke)L0oZOdFR|cpa13h ze_-LMyX(6yXz~?#e8rb8pN0BeYR}=`?>hUA_x{+cyyr;Hvz#kAM-vy==QrlpejobD z`hS}5NS-L(*-gGM{~xLIKKK_u_%rdFoxr(d)>D@l`#%}LtSG@3^#MDmU(>;Jvym5bKKcj!epS%yq z-sEB4t9Op%evVzR2l=0V;)K$A?7Ql`J8?=rX&(XSQ+bga|3rTDaUg!5-1$HI4ECWP zT+=_?=(iTVz=eJ*X@?Vb0N2Lx|M)-$Zi8p}(-xc?5BJM`eVe}@t=+uvES`if^qcWh zK8(lL)DN&X>yTf{v)Bhe#7@L#@(?%Z_h7ffrTTv#ug`g-d`iEY`?1kx9%0`pp42`y zc#^jtik|rEDgN0B+}F!KHa>WiJPb1UwPEAUKj^a~$w4;RkEf``J} zYG23Z{W-eui5?v5&tGd_4(=CwIcs}5ZZ33hSn(py*-@Qap);Kq?%DlPr@yQA9r&?c z{rC6zFRJfX*b~3OuaF=AV;uRbU!8aj@AT#zjC|RLI)5pTV1N2~4p+`^*bO^#k7oY@ zr~J0@>W8!bX}f(ge{LS<0{nl`&-pETuwRl-;)flW4}I%nKjJWX-q`hdc72!3-_Qs7 z%5T{NdlCQP+&CrN@YD1TziSKgzika2_}3p@b^u2!8ONW)9XuMZ|7EcU_%)6_EPpqy z*T=Y{EzjL2R(vpijN4w%@rfVvuZ3Uj)+tWl*LePfT-X=?L@(k8eWK^Nc6H}h?Bklw z*Fwb?b|rt{FAHDUzl%58*~!H5S~=N;t3Uf`Ku(hkqpCccorlE0B}$s6Dw9^g~HK=1m)1ABlMbl^~l?qZkY zto-~+a40^QH~6Z4#%I1Q?=Kq14$xBu_u`CwO6}+Ifq&~^7wkrfA95fU{)Jy6KXyTm z>dsH-MOk(gdyq%lufvmexD!|SNBGljKVsivzkgHi=k0UZ$1TB;yo zy?Ofwc7uQWY<8>N{BX<;tb;u$J?B^01^Mq!{JS!F)q#$~o#p2}r(gP$AJ8wqH}%_a zGR4i__;G$zX}tM@-^z>NSih88@SAzrL)nY*@UO%VyyG`^S?~Ey_>#BPIhN<>kq`W< zqpy7`XYwR(`r*gu)$^;nKmStwzk=%e%g5~J;X@oL{B%5Q?qlGFUt|Zy!6o|>7x-Q0 zLG1DR_;0xOTzp|C`k^CkINuf*#6|SYFCH3yFzZq$Z}wn){0sT&M^F49|Hg08D|xHy z7d$q&g?H`fPVo*WN_52!cnyEGUrGO}qt9;PI(RKR8ST64JEqu6(ZP4@rtTlv)oQov z9rgRvYu!(8t2+BgoG@?gQ>+Uec0m5x#d&sNzo!&O=+C<0#wW0Dbnw|BH`sFK(fS@9T1YZGHT{@#MtL@you*L!QLn zsIwRT1U|_feRe|L@Qfe$q+k9MuJnf|b_5@>o4O~N($hbotyFRiX<&@Ib%x&#vSdaBCj8XIJn!E#VB_?Qh_g-Jy>!{mahI>woqAq&QP{(%*x8 zxbB6-5%N0F^|`a&1F>Gu$x**4{^^sR#SOTEV|rE6hj_pa;S4=?Wt`{mBtF!9IWO^C zxhMEfeT&~WkA85TI($#n6Wxd)BpOJp95d`!JsV_V@K(RPX=N zi}CPme)RZv`lmm3qJ&4dhmUI(`d|H?t9=1`QObk(ZMf7Q?&MQ9T;Bg0{_t_IuXCx| z)zP&M{qe^>;aOaOL-~;L;v+ed7e4`q2SL5k$-Pi?tSnTs`?k_TL(J>FZG#-CSd{6tqGdUQ)H|yg!%tOBW zIvy^reZKghzj?@A`~KJq`s|2(z#}}!AH*BDF&~`67yZJKe(egZX#;^8uZvH=cOVKlZ8n4t& zU-&WK1XuE;^}*kD_6l|p{fAe?N;%CdAvd+Rge&Qd*mn!i&>}J;BP`mlG2iM|B zi(Sggwv=$DpLwQ!fW54B``%W0Kb-=lYZD z-q_Ec?icF+M=rdg_uDxyHr}}4hTkt7!jE}_vvuRpTUnU@ZL58}c666Jf7`d%oA#ye z=HL8ff7XScNxkfuyvS95{F#rx@LW8Z%9)*jyal?T8%I_mnZ&BeJP zPx$rR{7Um5h<$2@JM`7{gL8P?pLmErb@N&uynD_bjKjBaQ{2NL`IDFN>eiuu=<=^{ zTKil4!QIM+^ycz|Q)O@rSL~zeOTnR%9?=b6;IrmuAI8C1`s26wM`z&plHzl*&%gP8 ziMsj255A}Ln*PCsbrt>9PJc(`dHjHN;ve4Azs`ZJhhE^1Und8+q6auo?v0%FUz@9I zLto}&ckrv-eCV(b{OS)!%4<6Pg`&?M*ymb*zODX$Bkfn`eukf&;@!N)8;3vZz(2bn zFXODsJnVzs_s2f@2X%Cn#*?pp-%nrq8Bf3Fg$w1t>*nttWL(*Uan`dExh{S;t}A|+ zf7)+p-~+zlxNu88!7qELeq+7YZOOaLmvtCdbSKa2KAl|1fn4d4UaV7p{?hs)AAHvN zu;=_B`$3>Tf6y< zH_kfCKjfST|M;553G|JpFL*=`j_FM~m2>b=_<-+;Ug@*)5dD<+4PQkMe<|VJ{=e)5 zzf0fSF5sxfr62smxj0~cd~d-yd7-nmF#p?(r{ts`{jQ|nV(gt>$vw|&UUDISdP9$V zjE4_(`d0GyQ@O8&FC{y{|Fj)F^T8|J)4%!9I}m%a9_u5Ay+eQX-hp+R4#t_uISL~ zYM*z9=apyD2YU2nJ~$v}xG|1hsM9|@uWfKY=ihkpN0;8M%Rb`D#DyvU=$J=*L&tpL zmAZJN-~Qx3_$WJ^v%A5s{?>;N@|mtzU4Q-fo9MT2DDLP-e=C_EKImnFlhWU0eBl&7 z#-(1kKtJN#o6nK03zZBkL_cyApcQPx)%b7o!?z8{m5gBJ$XLH?PgDft0^C$ zUvdk7JnWylpeNXa6d|1zN>{a`ulphwq#zkHgN9e~q=#s;<+#9=#JjyQM zko@V@c=LEpZ`SEKy{-m7gb~@kLOq6MxB1q*B`x= zjalaNpl9E#j!$$WuhJKNm0k51hk35dzJMO#)p+Z`pEB#BmwL|*uGpRR(X;XR z7f00PORKSGj_=^=%Nc}=)l*Nu@iP zMX$@@8;;UmbnyLsN^zw6n@5RXb@Q8-9LD^oJGZU#2)NXLaje^V%TCE*Jm0;=U*I2} z@WW4|pR(*1ze|x99I2zD4o}uS@iSZA8(8ghca+`Y!@8#>y5>;=R;>-;$`M3zo1<|bjfwGuXlGNS9-AycCMX#$G92)E}tm5;oCT6_~7>{ z|Fk~xOkADzt9deiJrDoo@64Zd=vU9fPvI)^sQ5r0X|HjUc?vQ`%5Nw9u(=%-)d5zfXsE-889dnx`3U-W=) ze0t7~Gmc;1g3rl1ieL0HPss_r!Pn{hM3wk79{*GO67R`xvYx`jgr0VG054nY13vLN z=Gjf3wVv=*>&6d!nV(#hd%B%oSbD;L);pzN=jiZ5KE|16+HRh`(f<@D`ootv5qc%B z)joa)f7{j;=6~DTVsD?-I?RX8%IH_?K_~RJm;CWVzvLDEiyycp5A$j#|IAc80g?e>1b7Q2VDiQSM##mzxyW_-Pe4BtxpSf6p^qpqCFLqF>j*U*n#i$CL6`#N{f zPucf!a7Qlq3;pV^KKR;Quuln&@ za1vb9I_TMYBEK3B@8qRFIi`Qfmt9y7zpFG3oqgZEe~{;GSr0y}8^57b{zbcf_z)M; zUiuyKIlb_zzxj+KM`h;WKg(a3*E-eFRrfsWuldQ@Jo<+pa;|wJf9vu4XeB>5o6^fT zcx!Mp*0sCZC-nLI!Vfz1i{GqweenHO>nlDz&wSNBNor(Vx9o}8xs7$1di zdVt3%J?-?8^_F~%XJ6K@UH_~TzR7z_X&nD;ex>K;gDZ9OkNGy^kGgfh6S}FFKhqxm zt9^CEjb@(UrsR)5B|7F0PO3kBdu|+l_QgNsdFkCa{lmu;x7D9~mPdSBA9A~x9^nOF zo)^8$Q_n~JlBy5>>O2G8^6%D>vc}=Vco?YSAwJo0jWf^QE~i~(Z}?gs`Zr!#ej?+o zvv8krB}aNm$seOf|H{g9*mc%Xd}}8katmJY!~fPi6TQKk=fywy#s1K1)(bcK;SWBw zlLI-YeSL4_tv)S>oOW;FCwjpzIX3deU-+o@jH~|P*LrFma)l#t3jfB19{N+d_@n=- zTQ_?2HsN1BSbin*79HzEkDMb%;|oXA=jKoUqMLlodP?t8xv*DsC;S`V@*ECZ=Bqk+ zvhUG<7r7c=>&ZODFTIZT&GW(ZW2e?p<2;|}nO~B)HY@Kt_*o{SsM>l}d`%P&MPH6Go4BQCN#>nwbcN9HRz(Tnwv*JM4V*D?O|lFw?F zOMOpeUj6I5%DUkL9d+~Lb8qAw{o_|jF7#OPK@a^wZ_9W4H7?IjnzC-ScOtS|f& zoy>1|ZjP)HWc}K1_j^52P>KB^FIC|G!^JQGkJLGb4t#e=UMRExJ_2u9z>r$_A zp;I^>{kK=%g-_%0PcKUIreCcap94qNl^-(C)c=Z8o}179dhb{dy(`JDO5=jx7XHYu zem@ic)`=f}!Mexu&3gvx%(&89=A}3IH7~ps-{vzvz9;j;gZN%_@WsCkzRItbJ)=8W z&pO=hja;Lz!UKKKFFGk}9m}z6?d(|z2kdh6ulQKFG(PKE-xq(czBhcLi~r1Ty}s+f YU+E|M*3Mr|^`-wCT4^<8p{Qv*} literal 5760 zcmeHITS!z<6pfSu!zc*)u*d!BK~$R61c9weY9YTi)9y;uf+3GD-4jTUu-{mkYIiy9?ud9}Y!eJTPaxn~ zucdoJ9xb#!DuT@8~+n#E?M$9kK!1Jx~F3$*P{*{cJ>CNy3GyT5cFOc7U5#G4S zc;hDiwy}#p2`};g;Th(oPD@UUq&H!#&Ct7e#E7L`%uQ@F@!bCnZ{7m?EW6_``6YBe zKNIgS`TYal49E0txMOX#9Iza)9Iza)9Iza)9O%9S*u3G*tkVP5^^g(Hc5Etcj&tHr z>HL;0LrTH#a-aIX$2=D{*K@J6zleS~d-MC~WxFaFm(TyPPwwxGhgHO{tMexgia4Oo zpO$%Xeiz;qZJaak9g0M@^FKiDSNWVIycDJ9DCJk7c=M#*U``h7XcQ@v@Icy&*|!^{zdvfo zzZZDc<`zoPWSAOTn2gUtc6b=x8*UoJXrb?4@q8Zbn~i_Yw1QA61HK z+4AwQg9@d5ZZxK{Wy?@}qWFwYqK$jiqESpyY$2YBv2EFk8mXK6gXe+go%?GW&FgG6 zcdWab*5`riR5zUhuS6BMFToSiKL`)0tUFL+`hd%H&Y^texy_?>C%CW7du|=L&YU05 zDbHgiwP`%=kPUlLL4DPPYtm=-lg~3A>tdelPW+t58Cj1Kno}3}e&RV~J(sA>Lg4kz z{VMUoCa;t4e|2v%Ki^k8$0sST6!1FNrSrV;y~EF)=DnxxUD-$Ko%!T`l@Hg0&zt__ ms(V$gt9QmVpE>IPAW!OjRbIK@m)Buc9?NjN{=KeGxqmwwFN=BgSli!OmACz~`MKVQ-Q`{K)n>gt zT%%13P`2YRg<~`A0MBZe%n2o2CYTm(gJpX(m7@U4$ugUB~Y^gc|1PBlyKp<0r zXD$K+2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXFWG?Xb@9lt1|1;+>Ue5)Zeg{p@^_`t38j+&79kK zi3!B+aaVPT@o6`6fv)yEHaBx_<0U2#yT@JCCB~=S%mupI^Vr*|TTw?=DlW-@RY^J@(%H4{h4~zRiF6 z*P+dZZH_$U9zZ>^xAE-(f`r#|MI{8rnm9`J9=}3x8IK2|M|bo zccTs3&6)2(Ypr$2|0Qq#yV-x}b^L!uZ?ApZ@6vwHpWtt!|2cpEq4yK~{ePhMhdp=t zA4uotraw#IX9@f)fuAMtvjl#Yz|Ru+pG%=xoNKjU;m-w&DrnP zx5qto%mK5?t)0GEXWD01xA8x^8>Z5+7u{nyz;E>^9}@!+-z9 zxN=jE5%2Hz>s_^=L~>khg<=huCS zN0vWh{k4krHyk?j&p)(XDCwJb%BNk|fg9f3x773TUU_r0tLKSN^SEu`=g)lb;5BKd+MCyH$6WvQfvulNzme}_9>#}X@aX9m z@6hERZ~T!Sdt&G8kNl%dea1ZVnLmY(c*Z}y;1AC5L*C@#J-nEo`J`UM-_?5cUO&(H zQEwcjdbpr(^*KJd#5~Opef{ftd>W4(q#XG3@Z09Eq(^jW`PcY@d;FlIpK-~D{JqB~ zyx}+H00-hu#9QO0UfD&&U-Xscr5`@dq zd&;YiUYdS2$2YHt1LVNovt8hTJ&}v{>{h)z15W0WXJu#lnWuS>13vi8953uR`=QHU zgr3MZ;>Ct#u@`GUq+SUA#6MF1QJ$r%^e6SOx{^Xg@)_?qx zQrzL!>UAJL7k(q@1&-*69L+oR81u)s{_qbs`e#YcQD(pJZ~ApzW6;){EcJB8U*nTU zK1ZI6U-jtFXV!)`--jLD_Eux(Kem25_apE@&v2nXd3!I<@E(2q-hTMb7wx{@?{a;i zoqp(%L&nJ(*Wy&vvoCzWNsjx{S@S*7aYU)x8IQck$N1=wfASl7q<-~#&r0$aU+~4Q z&@~=C%0uZH9<;|VKaXDODR7xT~cCtvn6>5<>&zrz2fyxCi{LkIqp+QR|;<2#?_M^zky4|M2- z-z9f>hItu3`5LG`*cYs;-dE1pFfNJqil~} z#PwSK*{;$L`Jqo<`g_KQdBYPw1V`j-JaQ!8+TUx3{vE&GXwn7yosxEyb{6s=Z+c`W z#y6gN{Ho9R5^@k%*lo@?-s3CZ2ma|5ANm;=UGh**4*1cZJ<>z+&EHlf{m=tHs9#OT zIQd@w9k|K%F>iV?FLGC+V?HU5@HgtUM~^(oCEtfVk^{R#C)>CB_rKil%=ck;#!=!Q zK9Wx4&-99a?a0@Ac49p8^j>@8(XaV>5C8g+qxa$zy}}LLYG-`&!5{jb@oRqC>u(&- z`k6;9NA-6+u*P9CUV1OjSI1%pY{kRRX(r3 z_DLu3MvmqwPt#86J-*bNztVUar>gZB9OgWOZ}K1~&&J8|%`@Vt`WipAerx`6eDR2X zu%3zdMepb(zJm|)FfV+NtNu#;J)5tR{BDoDogDZN^vwqy?bK(x$amDkAGxB7FZjCS z+WFV-{_Qs<<5c+594F|ar=4-pflv19J-L~;_v*DHSMA7$e6tLFq&yB0N- zR)%(e*qR-B`!NjrET9#-ks38M@|!KlIrndiDF7 zPtV3t;^X#SYyWA;UzW`Cb8;d_?a?tGav^8-AdbSRd|rR^*59-FD#au3^=B_i{CPGW zI_5#IX(yH6%<-Z=pQ}98Jn^HQ`QoFt3wC4t+J3WN9S0(gvOnXKpLr_HlYZgGx=uTM zXb(r|k`SZ&z zmT>|<_@F2KjpMx%KlEpO9yoLZlr z&BM6rjc-0m^MaRruJRxAMju^r@|^ADJ24-)RI*?6*tZf6(IvN(SFGc;qfdC(p8W7- z-YLi6v!=`cB>n2#(0g*dbLMgPP8-?27(dInM}i+EoN0%S(tNaMZ{+2@dhN`wrmvmx z=-E8jq4~lO`$)SDd)Gh5jeRY0_FkFzvJYl|Nj!krXDVemvt>Nq_jw z`LTn*v+?Ky4zy3bhMy+~@@Kc?X}oL~c0x{BAO6cY^aeLc-~G7YTYa_<`s_j811IKV z9`22>yQCBRXg?P&gAe1Ad;Oep;1>e-b^DZq{oRl+zR)G#)Mwx({x+Grp3(jGxb8$M{BvJdJ~n{&)Oq+MVBwI z;LpAmy_+Zb>2Dl*F%R#JPk!tijy%&NJm_b9?Nh#id-Y1=^84De1M|?{c*Zd<{Hw3+ zC+UX&M^8U|uv_CBk6$5A{k10-_#{{T@ugqVug*WTSJwK+r+WI#`mjIs@C8qBn*7H3 z4Ls31{+05pJBKvyvg}Q(6!kq?j~uk#w26gW*b#rPozlLOcJOQbJ2!uG^9Kg)Qs~`0 zjRSY+m|xAexTSr{Bi7|Pe(VPuhkoEnJMYQQvv%f}c2(_D)bfElda*Brp7|JGdpI&a ze(6)bRF~ul$LQ!UPuI`5_&{HO{k-?A9lL-_CeWd7)=hcDxk zC;QQlzUfgr{on(>*aJBm7yju}JAPSzc91x?4hVbFPTasBzKv&o=%~ksd7$GxI^;yY z_+l62K@RyGxMmmhX*^}t$NbQtU-N_@sY&^9*-Mhd+<~RXApc?56hH{JzqA zB|P`o^OBbLR@keEH$5(Y@5;ZPut~vxulPRmt~K^LyK}oM3S6@je6lMz(oP;oPT5Z0 z6YCrENO?!RBR}#*2d?OWKIxmBJ=gXJ*XU{oH|WD7yD`4;*aQ5kN7wp7fAcdh{luLd zH}d@)zgmxIN1ojSeVld!zvvpj z*TK)++HL6hg?~YhTU#6Zl3lo08S7>KqxLiO>peZgDZSwXF8C$w=ox+qyL z+sCW}@c{Gadr`B(4hMLd9W_>{-Pd)*Fx;Fvw|=Y7sPe3{M9|7WR7_90%d3-O{) zzf(JRJnHb$x7EM@mHua4+-IG;pD6tzp7(uc@8@55cFWQ?zE5H&eKwl0XP>Q3FMXE3 z{~y==W<+7%;*at0i%dpBS5pbzt>U-kKHe=zbiejQHe+cUkvZQ{&%M&uDSUg#e_ z;Q&724qYW2;R7Dki&NrB(hq*pUV<-i6aD0~Ztwgd9@nLhaWlvnI)u(#Ah zth>>{7drSQKlo*z;(}*+2YaD^ewcm8gXAr6s9lYB_^bQ%9(>#ZM=bMf5fA$9vBy4F ze!o{?k8s~V)`g;oAF=0dvrm{-`i~mCTemO!m4Pd6Gvb&%S1JQz-PJGd z`E zfloMxTO}N@Bluzu{65^PhZl71H;FUW88!aVF`g2>;MIK5CpRUzWEuHYKKp)yUx+*_ z@$C18!q3wK`PF>Tvpj+x)RQ;-!?|a8fZwc-dm?ZmAAv9axA*f$9`^bl9?$)!KG%NP zXYiHHg`Mbk*L&j*={b9z(kF15W#mVFFPOdX(|=pBu#3LmF8AR**Bw&^meZQJfLy2El=KKLK(1s$N2`g z->dYyYuV|~9JoV~ci}^R2H*TeEgyEzPQ+#SChv?p?h8lWq@DNr!$ZaozvB}A0KVuO z{^34x9yn4D*V?lq_NJaa7>C{Rv+O4AG4c=3Nw4zX=xc|rag1BrD?3Yh#(s}}Cui1@wBG1WuDeM3~o_AaAl z3|#q&x7RA8SAXpINjtn$hCes)_{~O5F2grjqw|&@O)rDO&kP9vWS!Et`&RdSyXWk* zQ}f~%$b-KjKc#u|`~0H*;xWII_y`=Q9}NGeJsiWie(WgiA@Uh`@m~JJ&*Mu8ueJTa zse0qUF}mu}_a5EkTOJs8jz96nv-uepKFPs&a1Q_EWnTOTeWV|X^Ih>lzOO%h(!Y}4 z;oq}7=&temd~n6e<4d1ip8oiv_P;5K@5tl%U2#I5#jZTVf8W7JKk?u%+GZXFpM9S? z_3{;V9ajb{df|N?`W;mIPr3foCl;QWabobT>)v|hvu~H-BTqT+_J>a@qc@s5X{{Be zm!`0z#&N&9yj%P0OVa`8&AxDpM@!=oTi@KW)Dva&8?9S)>Uwn<8uv#AKiBGs17=)M z`mep?8yjr%UFo;Q8{ggC?KipJ7ys~`{v+Zi|G+NNzeSu7U+I&7%XV>In(xCd*a`ng zkBN`K8NDVB!`|Vk_SgJB{|-O!l=)Mv-xL2;9MT_u#?NQpdxSqj7eB_sr+FpLm36%E0mO zfBVCZgUZm8rfq*!=Y7klC$C%mwbpx=#>lr?LN8;!Uvts3zx%F?IpV~r=Ug$Xw0!^j z=XO1R%hK@3`F)m}J)n#{?y^HZUgg6wWbfHq9J=*y)9%C*{PP3k#INuZ{1-c*uZ-96 z{TDsZ6Mge%^avm5@w?uO3(BMy>skH0he!BKzaDbR`q z_mc&#Eu%)P*Qe|BNu~LO_n*CY)zM|lC&zzw_wS~a2~T|Y{68AMDihjtSaqFCKPzLS zU-M?qf4Jd|dzaB2uGypY8Y`8d;pfrszs$6k*17A(v|sBW>ksQTc0zvgCH{mQ)yuEw znf~ZaJ8^=Yz#ly)e(e7Tj`&qLG!FbE4guVyU8m<-`yBV z^rdI!KGcvnpBc7u>v8Lj>{&)XI{dJfF7K4)N58vwwNLLXV;|pkn=fY0SH{1(>1U&7 zFIOh4c-Q8uo^xRtyVM@neL8H*(!9jX@!zg>S{Zdp+p{+N$NXjZ!E1CIw9(aN&^~{B zc!P6y{n0KW&Y8b`GxB5?{EWDf{wUUuz;2_O2Zqqf^TR|N4^$yC7eA8vKjP^g(a*#Q!Vl zmp_*;X9RkPfq4xV(AN#KGterewecH*_GrFvI^0T>bWWAmHv4Mkpjy!4b=R1FQ z#i}QiAscUT+tz^|{%$?p|EfN%XZHTGHKJiX>V*+ z#E((i&U?{^T{bGC@BMhu6*u0wG)11)9OE``chy=8bR1S1<9x88Q*bksDdQAK{pLH%4aT31aNWKK`><8ZA&U<+F zUi_&0sb_bY?*u>UKpKa9rr=PIY zrOT8Nd+hSd?XPSo!{YlB`@{W%59cf9#ZRyc>vE-dEY4@?zH!{2H5qp!Tple;OOurxNI`weLzTjMZ@*cjG?3VwBclG+I&wk+-*pd10_tt0X zwIdhfrC*AB*!robzl;NckJMw_Ut%BppS(dFkXK}$6z6B`VL|}x~$zR zdz|}Bf%gFsANXVYO#{O(WEtmr_E!hI^ZjSreA2EA{BWI5Z}9W`m^f00XSwa z;)Hr7dw?_eQm-F7;8*!;_|^|j)n`12^>p3N{ayBA9QJ74!Jdt)z31dB#({JAM$dfA zkNo7L{E2$;GUK)TtZ|QyzU;HgkLfe*A#Q@_-JmbNIlS-N zyJy}txXS@k_G;R-42gZyA=~%d`TY4i6#a&7{L&@sT-#8F$9==0hy8lM9e>=e7|;00 zXZT_JR?Z{)N4_mTweP@>$ZwVQ%jrX2#ed0b`9*#V-oyv^m4CqtTqLgi-fzeoPBLyp zT<{FX>_MCoH~4dL0@ z9xFXrPh@cVr!k^N6^pa2Kee62+dCi;OFfRG<8|J}pc%S)g!nOBeo@A6Cy_4_SJ@G4Cz@av9R? zUkm(eQnv!<{P*yM-rey}+Yc>6zC3EuKR>#m41MPIZ(GjVtPDSHYR`Y3dU3{)#Cz;Z z4E^LUi_H9LddYFzuVNR`j{S%$gC98dm7DJBT%5DZ&%{CTj(y}lTEu_$kn8i_E{hM zT>9ZpoH0Io5iiUOUrPPuqv)hOBJW@);)#69I+|U;f671Zxu&0UpF7TzGfspZ@Qb;B z8tdo*vG2xC29IgJZov947I_Yw@x$e)x{cNAf`eJzR&+cDp=Hf+s7`9mJ?GEnrN!~vk^4@8O zZ#(}vc~5b0_^mpx5+AM0+Zyd%yb>6O2iAJDtJ zjD5gc#ueuumH&lD{#X8<>+sO8xRGUb4@y2)#|JpakMWc_ZmhS>5B}L3xg_2retH%c z@XxN;k97k%*|)O)M4s#<*Bg;%>7V&@tSj7eRc}3QKPlJ44=nNOl^={LgJYj)Nbt`N z>=z9Uzdd}r7Z$wo>~Up8tUpIAvcg9Po^)BRlZSu*c*k{?Z!ROAc<$hDcfKs+g>$)) z@qN?qH&pyskylO}J0kB{IrkXWxmVYfPy8z53cGY~cTmf)ou}Xb zcCqhe|JXjU^%1$qBl!(-5=X4ttWT`#=tErLx9lg@@xMo`i{Qa~cw#^B#*cc2H|^n0 zJNQ-)=gM4vhy(6b#d-zavy6Ob!;U!a{=$4BzOs%BX%}(MWBoxM;wQUMFFt2p72iehNA5Yf_rNZ! zFO<&fvii6xUctD;9hB+ z0FQ72_qh%Tf33Z^A-=$|^@H^WJloHLUv%JH|N5Do<#@4w#@@sM{qc>yxM{xp20rX# zkehv7_Mu+9*4{dsf6Dc>JUP}E&b!5hJkO5xwz!e|I&n`t*SE1>>V7QzkL+z?n2bIMAX1(S+?#azdW8Bki+~t4=H{5?_X^eY_4RQY1xckIeOMiKF zX`E-!Ye#LmS7}_f%O7uR|7vNPXP1q;UVKz(T=j&jR~vFx-a~IVWz7EHPoAFNaio39 z*W`C`-jVk1e8hfm?1$M`NcqHhjksz*TRf@PCH6ho2YpzV)NzHKXfMB#H}I=OPb^JguT4|tz3k~>OP8j&H{BHb$jz6% z^zFz&+mxnVH+}uQRo5y_zyD$EPTh|yjbFUK_Q8wKD!$u5w_)-19)9wRFY}%*K1QCs z{-_6DeJ-NjpgE%7Jf`4(p&Y$dK z(LcWmN11=bekdGy52tXL>wwTVTqfS5UO%P%TXw)+#24|%d;OH~jxT(|H@xG=yyQ*f zVm%?wsz=BElzkEYCgXf{PODy?Y+Y%cX+6qb?B9qRX%}&hY<+BrLHn zd-J-|bouv<>#TS~Y3?}k+o>mRTblMgs^b$KM`kI`7&q-D?%!Ic|M_fP9s5#n z4BwvNR-U6?o+d9sUzu?r_Gyw%><6Ng_8D>$AJxkbl}RVo@#ae&{0cv#Uc4Yrd93|I z{z@EAeZ_Zm?l;-br+?>R_MwK1Xut20pSCLYp`1%NhZr5d8#QuBpBb0DyIy{mCtqq> z=F`)jn)PgHj_)LzWBleX+FkssKm4tXjrH8v>+f8($EhclvCDSZYS6-+%9!V0>->1T zzm&1D4j9wz`6){rF{i%&3Cmv6{>9RK(NZnj{Njbu9Q%z;Zw_x}Aw$ z`+oV*q5u4GT_4}o<~k(K&7B)MPt5(^_+H(2y6y$|4$FC|{atY+&qMrfkTU!SyRaVN z58)>BsOoo589P3@r?9Tobdw_rTq|{#iq#b_Xp1r|6yM=%A z&`v3CdM}>v8{|eFbH!)h9m%6&zeK#CC-Gihpma}L9+BsEaqrc6n|)cO?}bM7`NidD zZLoA{xOnvV)ozVq|7p9F@gIM<{)fHxDC1jie8D{pyOwdMU-ZdsH(yf5 z#ku2{2OhY3`MtL)V;6eq(Q7{GR9ZeiaPMAQ&X@Mr6yHgSAC2eqdHTd5pBL+O_jj$2 zM!Y|~^|~V`=lv_^mEwv0Ir4FyV;@+25Ks6W`yBluPV!II7xu5@Q}*LM%bS((LH}^V ze&i+a40l;NX9#lu_NSu=-{v)oVYrjeU1pn+poDf&gVL$j~ujYY{bq)JrXX?cd z_P|c;Pmy273HIT4oqUHEc4ECqf8vDorF_Bpp7pBl?d0j1PX(SvFZ$F!?|5vLVjqg1 zZruLU@7B3<>s&{)Z1up@*^`ziV;29q*V=E~li%fyKWh0U&v@;DGU5HrcfDz;8_UEv z7oXU&&TB_*J*!N7=ge#G>-B7z@a_6%^|^70GNI4R=7Y~(xr~eRuCcehGGp3ubH3Zy z>Apw5++=cT{@}%Rhi(5!v2JPl>)%({>fY~59{Xd~Nh9a% z63p2u$a9uq@twH+tU=Cq+zW{N^ZA{7{63I%I(y)c#1rwHUh6ynzU9^4+izhf+QAtd z!H+zv#vgpcqvtH02gW|F^||<9UrT%Ei144|#rh}th2mCqj;l;7cvi|;We#C^btaZWL*ce^J}o?{P_H<^6r z-nXn??%Vs=WxK3)K$#rpHIwGqWXnYx-zgJ!8M)@E zxcxr=b!my;;b}Q)>9(uvxM69WGJc)0f4?-pAGbds-i`cQ+4TFPpUm?Y_jE^`+0biT z%hmavxqS(7)q2f(ot*6hvIqAR-Oq9ko9o4}Q}Kk~Q;H|G|AY(pNnFa4!vBj4@XIdX zQa`13o|V?s;)S@7aU#|cSzpJWuse2vFMQxD?I+e5`~*MX97FtMN9;rV5I5LKp7({H zwC#7Vf=EpZ@2S z=D4Rj=Hi#8?tfPMGVZQcU1qMjdYSNxTaJ04{jz1^^pWqa_QTHQ-d)O4t=HPAOgZ?q zHAmgiP^P{%@ug?qe7Z~v`f>@>@BWK?wC|f`(srAC+0biRVOQh39DUI7 z*G($p7XJOI{a;wPwA?vh+xd5%w=|zR=BX!6y0$dO_g9UVw(o!Je2bNa*w6F*;HY_4 z-tFE6=lqUO++Sh;_NnuHCvaySV!zkD8uFIMvKRa7&Yi3?+^@1Ow{E9@@j~7W2jW9r zpLt1~b0!XBT`bjPEeL-Em*tU+&)M!DUW+@fT&% z{WF$+=#ExpO5}m}fBf$1^UglHJUs8L&UY?yLwV?LTTD3Uu*UMh`6n+hYw%uW>i(NH zygT~-GWnJXH@x)2=VfA?%TH)^-ug%G)2oO}V>TQ5!4|Yhnw@q3H8 z?(lnyu|G|3_EYSm@O$)Wok{=vDty3`cw(JveIf7RpWz-p`Fr?PKL=AVo5%&Qb*648WtaX+wjj?~+^yQ(O?(on9>9?)F$HaF_Jm8#T7U^+gnYPO#Gw)pU(6oo?ANPBw*G@Z?$CjOV#XBPwE)Ru2o7U~}u+H}4H&c431{Q99c zt{(eV$@NFvyJRPRpJDWgdoKLX`@hX~hy2QYtU8{$m!J1FV?WM*uDpu;Q$N19jPDZZ zkzVCfiHpcP;irxp*2i$DvR{O@r!4LlBh4|M_ zJk(CT_&_e=p!4uLzS!@OU-3)s!HMrmc2VDhNc-^n1@05Zd8GTV{HQ#{?+W@J)$ga- zSM~k&n0ZHSaP5+F_Gufg|JCaI9aqK=x#ZJ(uKKV{Yd+_yxZ9d2Fhpfc^x zZ|vH<;-}@|7=X=f{cUoy`d-?h^H@_l(Um@)x);0E{?K{}7 zmS4IzV14fVRO!1H--V^!M7++tD(=ZxC#p~V$9-h$cKfq%09ToxgudWcdwyD+U>EF+ z-RS38JP`-1$JvSR%+clV#Rv8V_ep={|NL$_yz6IN?YtKk*a3S{pZ+B7p^7WyWFOtS z&w7uYE&_B+8ZoV7dF89B1aMM>Ewl7n#lH@;HtJ!{w3-@5kga`&Wx zdv3DhYh~QYSA6m8l<)I*M_S_d;F|xs+jgVRJ+d^$?+e-Abxt#Sp2w$d@Z>_py%TZ6 z?>9LA@I4AU85n-jK1Y4-VZWp_;h2!$1DS6?t90$9LT8H!ty)eAr3msqy=@;)MG= z;?hlcc#A2u>O!&)#rYBUp)3d`Ct2_&i|Ts+w`)vkG&=HtCsjad+euw zf9lYEmMIfDEcVjN+b&n`UVpPkHvG6#ne^vDj~vwY@Z6`G+G^>2XDzd0nf9k8`?id~ zuS|`6ZOWTtSH0o()@AZGEq{IGnpx$Z*S{Y${ocFF`1t+8v4=cyPs>RU6z3t$p+CPz z(%5gM6KDVFmLK=8;+~7&W#9+wU&||Pp|T5{x8?( z{9^q5D)B;|ZQYUil=F;O&%!zU!6|&hv3;S88?j!3-#jOea}W7dO~6g-06eoR&+LSKiW};! zFO1K=^-q2x9$A;O3-M1qxshkuhwsAU{to}eKV|%g`+b=Q#QwAVz&W6E5BpZurS>^} zrza2bJ5bKa8(!~sAk{;JDbb6PmZ|ZwAF7f<7T|^ z(c0Y~D-+^-zq@CAI$`F^pOuL#U$^-?AN;;dJma%({`l?ea(8^kG~u+9_Frs=QpU$Q z{@8C8+G5?C)-3FzWk|d42i*NZXtA*wzBjtG+|Tu6=+k$hnP+{aNdC=Lztb z>-7C*{OaxgCzqD9_iO0<*&oZ8xCiQ9)wnos86S8Y-~9ak5AFV8{_d^%adFRWOxW3& ztry?%Pw&1^TH^N+@_Zxiu{3_SfByPCfbE)mVS?d-5yxA)c@U_L2C3o6sBlc`shTHGILZb%uI&!#>0r>vOoy zI1%G&U(;nD=o>HFS*M3SwG-FmP4YYS=E=V0TjVFtN_j_|&U5nkZrpyNxX+Kt2dpo} zi`>tT@5P;m`5so@q4XV|-#z!ep5OgvAB|zJP4C^g?LIs7FU@~A^pkV1nOK~syH_}7 z_gimn+jMwwemFMn&y9`WLmd03bFY1=-Fs!s1(%I_xXn&K+C%K~Hr=tz1#3>|UaT*C z|KWQPcI^I(djRgS=XZj4{O^Bv=I@U>|8t(cv z;+nX}4)kXq+OZq*;a@U-L|$N@N-2K(ev`h~hwnS+mH%^}&b^4in{}UX{MhS%{2iSO z{&IZB?H(+?2OoBSuX(RJ@s48uD)05h_pGC1-P~~F?akfyKCCp}cI8_?d@wEVCCaz* zetNaP)Dn50e5@sK*gSrvz7GzWbKfxT_crx!J@u7KFUjvj*kOKG665FhxWQkZ+r;UC+Hf5~UGlaIq4 zyt60xggbcCPCqzzZ-iZlBk*s0_F{daG!A?4%wE`?anWIy?3^8tlf29NMV>{T?v3zA z@-2C-xWYcHM{*xN?pOGI8ox*7{Fi^U4)yyva3sIVI1=ah&J9NQ=>6f<*Uw)Xj{5Vr zA6$NZe#htgJ-`3l6!%xNynd0tU3OKU(j4a+&fS~-c}C~IUo+?Xz)!}#waTU=f3&yw zUeJD-@1UH^jKr6HXy^X1p7q^t?oY*cg7)w1d%LghzKHKA`4|3&J-FxQJSx9~ioB}! zch(t7afQ9e8~Ay76Wqg*dUhb6_Y9xz6S51X_wdgy;F&%6{T_H`PjHT&e&~6&zknW` zE7<|N&_18z_dM7GdGbH%*$4k5uCt4b`+-BhKf?~_nceuEGx)GR^}AD9M*cDE&gPSk zD+lL2tX#LpdY9jpfBD^K=heQ4_1$b^+$(8p{Pg9uKbo&J#_y=Xzk7t<`#oIt=(}_K zwC;uZeNg)g&N0l#K9t{)uwU!<7Q_wraGbN|z1;Xtki6^lhW!-rUmnHp@%#KDe6SPi z4|x*5pY{>!b9hq^r|_#kTziIJ@qz#5_x0!RwO9IX(0J_3IO@}WV%;fkQg2@R=Y5*E z=VqUm{MeIt%zns?U*V@RuEe=4KbGenaqnTS{VRW$HST@*E-Ke2vF~dg;@l9P*a`b^ zUShx0_Y{7IB+uh~cOK{8;!DH5hqQh5kLQ+#xF3_h4;{Y)?Ed5E>!1G3bB{LU@08fD z_4_6MPN3gS@m;Jk&oA6lk9#%tF|0T2?~$80l6jTyn*8o+?BhDW;>YNXetSm##2(bk z2kdXw_5n}uC2oipaLiud&N>}lwX47P@C^U@u@iAjJYWY(^-6rOH}ODx%zFC`aev6X zw72eH|N4_FemK{F@%!ukZl?2^ zwCC_c;-YzqANKReLtM45;XbliqhBSxCe8yF?1CQ`U#wr@1zzEZ z-M~LwvIFno-rpOCfAPb6W#T^W31qwK93=C+*w?atN&cMA$G(esieu~v?(r{**@F`Pt?Qg0%D3zTv!DF?t#KcSUzA6Anxidl&INc8$yYz4Nid zAL-YsUH&ec??(L2o8O)C`-pyTCx15~_9Lu&{N980PsWE>Kjgj2iJcbP@66@$_ojWP zVZYA&>`(Yk%spm#5;nGPJHN;{P~9hFmr1w6yK`~&$-dbiTzijS_K@)+&b!zLJ4BCM$j!6;Joduh zDET9H!(UmSkiUH_zw>6@ZXH73@@oFm{YCM^{!Q-x`knOf1L6rko#!R-_ttzTlKD>j z4qN79?gPjDRQqFTU$M@~`Z(8ie(7AseHit*jtV>ReYm(GkF?I_AFLz9195|&aqq*v zEd0X(JE1S@MgEO`YaZ~AF8;(P?Z{6%^Wyi-pZ}=sMZQX|)+O{rj@BXkn)@ZC8{~`B<>+eV6J6!8uzk}oV zqTxD!e?9IIu@~zd-&r`%*Pb2b??=RWX7U&J%5t3@{F$HMpWq+z{+Cd_b4I7>|8upLP^}6zn-ax{=z*q_2idj+*9C}$ltSdBR%r(^vPa4%O|Wu(oW1WCqn~&J@7`w|VGnW7GVvbgx9Ha6n>V}6xZykBn3wtDn|-mHx*vN< zdvp#G=f3bSe(-1f7(0P?@j||Ay}{4H0e|S8mV5$^(q0z5@V*ZH4l3I5*Y2x0M`I_> z8Ra=}te)MV%>Ao&s<1EHC#A&z|+K>&2a{kM(0+5AWiWJWGAr zQ|Kl0tMGT)u~T^y+_PKrvi@Kv>dk{($q5~CSR4^&=#6}Px`z<=-Pi*?vLp5Eft^@~ z$|uAN_atgNflucEO7X_|nf0#r8Gquu#XgZd%)NU1NO?XN-y5)p%>Ux|9>fWy`Q*CG z-!F;%Ci!HpBjWeqYk85I^9b@weZ=qmh#UN$bq4zoAK;q)6F2qm-b6m59lXH_JAiBc zl>f&EeDU+_LjF|aTz~E0yRJ7bJc}#%LSO&6_=q@;ui6i=H}jBR!aw`qhuDeuOdjF} zKW4pX-67t{ueGc52 z*2(<2ar|B_|DW`Ij~DAH>m~Ule%OKW&^JH!lW{nHAA%i-8~hD@q{%RPAM%2%k8?NqgM9$_O1#DIl8 zd9D)q3;bn$eCN&o!zp_~2QJmaJNsm3aLta=|Hu6f{3_9{^Cteleh2xay~H^TzoZ_X z#TW8d@;CA=c?0~z1HYU3g5On$bvZnwJ;c5`I&i5(7k>G9b|XF*7rwO9-+1-?Mg7ta z;`Jlj8^7iHR2+*9Qbt=Hi=>*Jo2 zxRUFCe+MD%Jz0>%eCe}|1-c+an}3-TAoGp~wu2Yq-?zw{?RNd3nBk=k#=NsXJi{BW#) z^*0_I%8L@m@javVHGYlD9@tCvuj~Ln@Q+`13D-4W;sCq!9=(*0eXZ~l=%-yo{-i&8 zP7`=KguJm;8E#isx@-gq>Lw*K7@(OVVZumuZkp0}xi|<3l1%B6l5PRgm;o5ug zg`H@Jj<{ieTU-zq(Ko($sK0r6hHLgu)Td*M&ulibl~C-g<%_<>LQQn>93vfwS#~7)Gz6YM-f-De>}6(w4eBUF4`HNJ&}XFj$E~?>ubH} zuiiNDtG?C~eaVyfNA?5v?8AF_N62fuKsAG_2Z4voV; z^jEJxe3Og!{5YJaJmPmB)6YblH7LvURKbY_19%uTy>b^ET!U6rU6ZXM> z8(%;8gEMxS_=aO8JjXglf8!~;#`yr8yI10w-LNlqlX;zU!#Eegzj?B2^FYV5l6-6V znJ0gd_>XfV_QW62pZV&?uPV`nA9)xYkUzVUXTYs-;GSKyjxepZQ5b{eiy#9 z7ypd!S)55b3Hv1vafDrnpW2&mzK{J=?GoS4tKIjA{gw1@aUKK5^qOU?6T|~{2xoA^ zuJrRh@f7EG@Jo;0E8*00?e~qNUE)0K2EJ?iVju4NTAyf_W#nUUXx?cLl|00~x;oPX57pn>f*a$(IegJC9BBtXO8L9~>_|WH2R*okd$`7r zIFtAef9^T+viOc%|N4FM5&QJwEIApEJmH-DJe!C4>o4AG4^Q-|zfwJY;Ya@(SMsE8 zvEObzp!6Q@#F?ZY`#$W+JmDQb{H*w*^gX%p;2nMAsfTOvgk5Gk-%UjcXWAKu9pD>X zc2Az?!=?V> z2Rd`L3p;{uIEGvDfnWBZoxBe}#!=!!KXmK&DG&QxVJGBDkMv@m`tg_MhaY^_a^eTs z5Blbbu6lY?!ds0$`A?1Cj0dqFO0Mh!Zt%mtYMhgY{&37L@rQ5q#@Amz^GyE`>o@%4 zi@kbAUw^nIXY}Zw-jubR;Q>DAuhuJl%0skM;tyZw;ZF(g*8l7Tzx=ay#x+0vYI}t9 zdOrH&ONpNG_$mF#Njv(5b9O+l_T$KYGT+7d>UZZGXi*r;vko zso&Ud#RvMyU-YwI>b?2HopF;+#3^`%Gq^TSCAp!;e&8IQ$c5dQhx*z+*;Uey{LlR0 z7=QU}of7-J=)kv<9PmL7#>1C+c_v5C^vS=XM_AAkvNF_$#c>V{>d4?SyuaBTOE9D>l3_kTkhrQqnKkU`~*aO_c ziT?QTUj1Br*cYqL*R<1*eagSck3Erx_W04Bov7EJ9MDHcJA9aL`dRvjyqMo6CvqhR z^z>&3X@4vKXI~2r`APk?r&s#YpPupQSv${J2JRE5)%q7*IKU5j{xSRrKTE#Zu3Fc_v2nD|GW17Z#)m8O@@yXDi9Y?| ziypKicX)<>_=11!%pcCtsmIHC$2>eMjfcOQf8*!4fhYa>F>)b~EJGjg0T1R6XG(Ny zTx7qnOXJcLzSVnX2cC@szv<^A4tQ@q$#)e`w38p=S3mMXpMJ~(U2>>P`b~Mb2O0ag z>fu&@dD!FZ$Xiz3|uQ zSvPyeC%WWCU-V%d{IUn_^(QxaM&Eemm+}nxpp*6T6u65#fu8WCJ-m?v{_varDeTie z5x(U;=z0&o@Tg4r;4Anvp7(G;4*K)|?2jE7Pk;P)&iPjRWLY2ln6L54pPbAaKlJ69 zK0J#%=ApfK2Oq|ff6zC)>W^RiDd`t~@C;XJcfN~`^?&w{`HLg?&v;SgTjZo2I>w{V zI=;}0e)vu~)bS+n44>p{JpIu#f2Hxr7hY1I&JnA1x&HJ(?;whMouUL4Rb`Q`7~ zXW$V|@Kukazmj~lXAjrKeK_?s9qrWD{NdaCTE66u-?XpL6TZ-eclhjB>8Dh$A35lU4jgEwA3ECC_#kI?h7P`esuV}GH;xh?`7F+b-|}93(2gGHE#+FR zXY@C(#JBU!%5Kdw$BBK^q#u0o2k4R)e#i-(d>{La_^SO3+@ecADc>01ICYtH-FJ=> z|Heyua39QmgnfZnkEcAspPC2!*LhXaw+;$=`~p9q9sco0Kjv@V>@e+Ieqp^3dPEmq z^v8ET>mTzrE`H39JmpVyyBsIR_a6WJo$<{>yL=z~z!iCrKYHr3Uo|c}CkOJjzn0@x z`?dJX@gfg}WB!DEw1e*~L!RbI-ul%g`6}_}nVqBKJvp$Cnjd_VC*0r*zv}6mo{g9M z2F{G5#2MS7k&#aGQYxfsuT`Zq81(80H7{pim;JkwXo-@U0Qy*D2H zpku!DpdEituKIa4z7n5st0Z6j$i=wHA3KP+!hZE5H+-U-e9CL8{3z*%KFyC_@L?SO zJm(wzjHli+9HHwyK0T*iV?CMmfmise?M6TE$ytBv5&Bf(C+8=>kCR~<7kf`{2G_O%!mD_{KS_kpVH4f*e5=+p1%+K z;_u-H9s0ws`P$Dy7oXnGRbTOG-q}w+7=Flj`9AuSlm6yoJmb>?Jl6Eo(~J40d_&Ld zOMi5=H!ggW_grPrW2ZSj|6Gq3>$kba7dJwG=5HS4BER>X?~NDz)#E?sDGu&yG>C!doy$9&+9eQ2L@tmdsf`OpWyN{{e@FYjxbOmg7`SG9kUQ8}jNL-eDUjK6Rc?a{-3%@4icpB&+yKJbzC zq4#Vben)%bkSo2Jzjo${5Anh*L5D6dz34aB{P`*SK==&1 z!Uy?fDK7GNQKIjD4ZAkK`k6oQ9?to5_JR*3`6=~-UveQo_e$`e{w?xv_)WQm{<6J! zhI~?Qpq5NB??!@-TmNYJX;&lq3Gwckrzr+}3o-4;}JTuQU(*d1lW_ z{K6golJX5)>W@FT_RM~h?x92f{6pJ?it*{u{Uh|~C-o8XK(Er*JU{>YSpxq*O5lG1 D+owlr literal 5760 zcmeHIdr(wm6sHVBC&qN;WtKZ6Es!h^A5fr2Fo8ALMIc#0fh`^~1WOC};w20ZM-bz+ zJa@SZEGq;k$jstKOw58a2!=dHBYY+=vw{Ytgo+66xj3}WNaLS-|GD@3?)iS_oZtDK z^Lemb-CXw4={xBv9#6VC{WE@e5MK}$Nsr>w9a)|p^vEb)cmPijpu&D`%Yo(T#^S2P zGuyFk4v!}g@J`SJcu_oh)QN~d)hLU*YvU&7~8*Q3wYd@ zhR1PdIkDK!iO*vB@h#vzC%*rIx0n6VGT2d@su@r-pk_eLfSLg{18N4A&j92U*F4;( z;|$pi>%=3IwUE80!ywX(0eNmy!I?g7ki*fsJIl>ov?KrLp1ki=w>}24nvaPm`o|#G zGDB*fD}e%$@h9hxyh+*P`K5MfkNyxZ=N0qsG^AYWyl`0F6XbQD@{=!p-1q3p87jWQ zzm-MPbx}~@@0j4SWjmlf;^)5?cY3Vr5oMn@Xej+q^XX!qt25RR5ae>%oUYH^1Z?y4p~5{*se!&L${M9!onJx|Q0`4?6eMTdx}m3_JB^ znN5IxmhzY)!{S|VRR=@QnjsO+^K=g`kZVeD+5%Jsre<-qkVBKEn1cYzZ!C8}4 zKsfb|C$Khv#O_RDM6(fy&Iz(xAqu z(st{*UZ^#>?R?YhC`ie=;W^QvKH=G_7${~ZHjcf;q4>o7PzT6wsk5ZY;95OY{Lz)3 zm@0;^j&Xk;@?=5%h}(^u{U3nbD)aczHFs#-q^-FA8UYs_^!GRDCqUKaRS6o**P-l} zw4}oBcPZXcZ>TSeb3=aE1n(AOBRG-2m0Ssnf=Bu^`=OR};j3 z8_LK&FCpuIJfe=!FY2h2#4i);cJ_HYP&`+)8~NQD7=yaB($?Eolb|jpcug~HHB@Ws z?LNet0BH(y%KA98^p^Z^@Z~v=O+5cd4+w9=c#?{dvZ>rPH^tv8*!J;dt6ig z$a@Fl!Fu3%yhoPeQfu;F`m(9h_^*t|)P0ag))W0uJTJ^M4fUgDK+S-ffn_rA4`l!G AjQ{`u diff --git a/dataset/imaging/simple/psf.fits b/dataset/imaging/simple/psf.fits index 3d342eff1fcd16b4b65a921b9dfe0059f8d5cb75..940efcd879ec6e65b85b23feb4bf927982825a00 100644 GIT binary patch literal 5760 zcmWIc^bPQFRZy^1zyd-P^c9lx^V0H*a*Gv8@)g`XLxL5GOA_-^5{pu>s0Ru;d4>dd zMqo8h*USWHUSd&EVx>Y#Vo9PxNo7GQc2f|devT2I!9Xrn_Zb1r%PY-IN-a{zPlK79 znVXtdoSC0jj4%X=1vbwRyLpB{aZ&-;JR|Cv=j!J|suM=lkA}c#2#kinXb6mkz-S1J zhQMeD&@lvTn+vTYDyNCrzK)pC>Jq?ZXX7*No?&R4U1^WdvYWkmb~98$HaoK0+bw=} z$wy;)7n~21hp4{|(~n_Zu(fHWO~!1v{(F8U7Q%cA_6)z`{bg4@w&(u)<<`M}b6|Ws zh&&@)Jybu;e5m{O!{xus$ZFd@QQKZW=IJ_%BMSC`6CbRy>y@yN`MzJGSG*d^SA@zN z*g@69^+VjhAMRhM{V?_R=@0GidX!Y!cOA^$ec0l?{WQLmsYlH_p!^qbd5HZ6Q1>gq z&4Yv=Gt7T*_dv{hVZZj5T^!q5QTrX+BGgUZJ%RIK_Cfrk2-6P_e`tKb!w=?PnEPPn z9f>z6Dpp~In-7zRsfUI?JUnRSURZjecKE{L50*Zulpf&ejmqg4BY#jSzcGx0(GVC7 Sfzc2c4S~@R7!3isg#ZA?s0Ru;d4>dd zMqo8h*USWHUSd&EVx>Y#Vo9PxNo7GQc2f|devT2I!9Xrn_Zb1r%PY-IN-a{zPlK79 znVXtdoSC0jj4%X=1vbwRyLpB{aZ&-;JR|Cv=j!J|suM=lkA}c#2#kinXb6mkz-S1J zhQMeD&@luUPykaA$APmQ=5rq}7PzN4vwiXsPSHIU8`tC?xi20)XPe}Q@AWV~M1D>q zOg+RHrXm#M7?_H97(`D`oF0@4(=T6f@Zrm|j?Y!X4Ro9jr(RI~wt=Z)bHyzfA0`h| z57Q4bA7&p+UTOPgEq*_5Q?0A^k1A}|Jka0#Tzu{-)vt!Zm12F?x6eTM>Q`X$F!eC~ zF!Nykf!PmJZ}geBw0>1nqWQ%e0kzqzO_qhF>we4b>w@uN@-X!<{pjHb^B>GTF#R?g zwjG{#V&g))fT*PEO?UWUe3(2;Jxo6={9*2gg&)knF!#aCvw#2p|9x+31{fbE4^t1* z4>J#DKFmIt{V?~y+y^rc#)rwn)Wh_{%%fHK!r~9+epvj%;sfFy($fPhy}{BS%zRk< z!u$sdKNue-4^t1*4>J!vf56<2o?c-3(Zdg>9=my{rjPPQLtr!nMnhmU1V%$(Gz4f9 F0sxrPID-HH diff --git a/scripts/imaging/likelihood_function.py b/scripts/imaging/likelihood_function.py index 56a7e24b..483997d7 100644 --- a/scripts/imaging/likelihood_function.py +++ b/scripts/imaging/likelihood_function.py @@ -82,6 +82,34 @@ """ aplt.subplot_imaging_dataset(dataset=dataset) +""" +__Extra Galaxies Noise Scaling__ + +Before masking, we must deal with any extra galaxies in the data: nearby galaxies (or foreground stars, or +data-reduction artefacts) whose emission is not associated with the galaxy we are studying but blends into the +field. If their light is left in the data it will contaminate the likelihood evaluation and bias the inferred +model. It is too easy to skip straight to modeling without checking for these, so we make this step explicit. + +To prevent extra galaxies from impacting the fit, we do not mask them entirely from the fit. Instead, the pixels +are kept in the fit but their data values are scaled to zero and their noise-map values increased to very large +values, so they contribute negligibly to the likelihood. This is preferable to removing the pixels entirely +(e.g. for a pixelized source reconstruction, removing pixels can produce discontinuities in the pixelization). + +The `simple` dataset includes a faint extra galaxy, and a `mask_extra_galaxies.fits` covering it is shipped with +the dataset (created by the simulator). If you are modeling your own data with an extra galaxy, you must either +create such a mask using the data-preparation tools, or shrink the circular mask below so the extra galaxy lies +outside it and is removed from the fit entirely. +""" +mask_extra_galaxies = ag.Mask2D.from_fits( + file_path=dataset_path / "mask_extra_galaxies.fits", + pixel_scales=dataset.pixel_scales, + invert=True, # `True` means a pixel is scaled. +) + +dataset = dataset.apply_noise_scaling(mask=mask_extra_galaxies) + +aplt.subplot_imaging_dataset(dataset=dataset) + """ __Mask__ diff --git a/scripts/imaging/modeling.py b/scripts/imaging/modeling.py index 278ce865..97254637 100644 --- a/scripts/imaging/modeling.py +++ b/scripts/imaging/modeling.py @@ -105,10 +105,40 @@ """ aplt.subplot_imaging_dataset(dataset=dataset) +""" +__Extra Galaxies Noise Scaling__ + +Before masking, we must deal with any extra galaxies in the data: nearby galaxies (or foreground stars, or +data-reduction artefacts) whose emission is not associated with the galaxy we are studying but blends into the +field. If their light is left in the data it will contaminate the model-fit and bias the inferred model. It is +too easy to skip straight to modeling without checking for these, so we make this step explicit. + +To prevent extra galaxies from impacting the fit, we do not mask them entirely from the fit. Instead, the pixels +are kept in the fit but their data values are scaled to zero and their noise-map values increased to very large +values, so they contribute negligibly to the likelihood. This is preferable to removing the pixels entirely +(e.g. for a pixelized source reconstruction, removing pixels can produce discontinuities in the pixelization). + +The `simple` dataset includes a faint extra galaxy, and a `mask_extra_galaxies.fits` covering it is shipped with +the dataset (created by the simulator). If you are modeling your own data with an extra galaxy, you must either +create such a mask using the data-preparation tools +(`autogalaxy_workspace/*/imaging/data_preparation/gui/mask_extra_galaxies.py`, or the manual +`data_preparation/examples/optional/mask_extra_galaxies.py`), or shrink the circular mask below so the extra +galaxy lies outside it and is removed from the fit entirely. +""" +mask_extra_galaxies = ag.Mask2D.from_fits( + file_path=dataset_path / "mask_extra_galaxies.fits", + pixel_scales=dataset.pixel_scales, + invert=True, # `True` means a pixel is scaled. +) + +dataset = dataset.apply_noise_scaling(mask=mask_extra_galaxies) + +aplt.subplot_imaging_dataset(dataset=dataset) + """ __Mask__ -The model-fit requires a 2D mask defining the regions of the image we fit the model to the data. +The model-fit requires a 2D mask defining the regions of the image we fit the model to the data. We create a 3.0 arcsecond circular mask and apply it to the `Imaging` object that the model fits. """ diff --git a/scripts/imaging/simulator.py b/scripts/imaging/simulator.py index 45081114..34cf66c2 100644 --- a/scripts/imaging/simulator.py +++ b/scripts/imaging/simulator.py @@ -1,225 +1,274 @@ -""" -Simulator: Sersic + Exp -======================= - -This script simulates `Imaging` of a galaxy using light profiles where: - - - The galaxy's bulge is an `Sersic`. - - The galaxy's disk is an `Exponential`. - -__Contents__ - -- **Dataset Paths:** Defining the output path for the simulated dataset. -- **Grid:** Setting up the 2D grid of coordinates for image evaluation. -- **Over Sampling:** Applying adaptive over-sampling for accurate image simulation. -- **Galaxies:** Defining the galaxy with Sersic bulge and Exponential disk light profiles. -- **Output:** Saving the simulated dataset to FITS files. -- **Visualize:** Outputting subplot and image visualizations as PNG files. -- **Plane Output:** Saving the Galaxies object as a JSON file for future reference. -""" - -# from autoconf import setup_notebook; setup_notebook() - -from pathlib import Path -import autogalaxy as ag -import autogalaxy.plot as aplt - -""" -__Dataset Paths__ - -The `dataset_type` describes the type of data being simulated and `dataset_name` gives it a descriptive name. They define the folder the dataset is output to on your hard-disk: - - - The image will be output to `/autogalaxy_workspace/dataset/dataset_type/dataset_name/image.fits`. - - The noise-map will be output to `/autogalaxy_workspace/dataset/dataset_type/dataset_name/noise_map.fits`. - - The psf will be output to `/autogalaxy_workspace/dataset/dataset_type/dataset_name/psf.fits`. -""" -dataset_type = "imaging" -dataset_name = "simple" - -""" -The path where the dataset will be output. - -In this example, this is: `/autogalaxy_workspace/dataset/imaging/simple` -""" -dataset_path = Path("dataset", dataset_type, dataset_name) - -""" -__Grid__ - -Define the 2d grid of (y,x) coordinates that the galaxy images are evaluated and therefore simulated on, via -the inputs: - - - `shape_native`: The (y_pixels, x_pixels) 2D shape of the grid defining the shape of the data that is simulated. - - `pixel_scales`: The arc-second to pixel conversion factor of the grid and data. -""" -grid = ag.Grid2D.uniform( - shape_native=(100, 100), - pixel_scales=0.1, -) - -""" -__Over Sampling__ - -Over sampling is a numerical technique where the images of light profiles and galaxies are evaluated -on a higher resolution grid than the image data to ensure the calculation is accurate. - -An adaptive oversampling scheme is implemented, evaluating the central regions at (0.0", 0.0") of the light profile at a -resolution of 32x32, transitioning to 8x8 in intermediate areas, and 2x2 in the outskirts. This ensures precise and -accurate image simulation while focusing computational resources on the bright regions that demand higher oversampling. - -Once you are more experienced, you should read up on over-sampling in more detail via -the `autogalaxy_workspace/*/guides/over_sampling.ipynb` notebook. -""" -over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( - grid=grid, - sub_size_list=[32, 8, 2], - radial_list=[0.3, 0.6], - centre_list=[(0.0, 0.0)], -) - -grid = grid.apply_over_sampling(over_sample_size=over_sample_size) - -""" -Simulate a simple Gaussian PSF for the image. -""" -psf = ag.Convolver.from_gaussian( - shape_native=(11, 11), sigma=0.1, pixel_scales=grid.pixel_scales -) - -""" -To simulate the `Imaging` dataset we first create a simulator, which defines the exposure time, background sky, -noise levels and psf of the dataset that is simulated. -""" -simulator = ag.SimulatorImaging( - exposure_time=300.0, - psf=psf, - background_sky_level=0.1, - add_poisson_noise_to_data=True, -) - -""" -__Galaxies__ - -Setup the galaxy with a bulge (elliptical Sersic) and disk (elliptical exponential) for this simulation. - -For modeling, defining ellipticity in terms of the `ell_comps` improves the model-fitting procedure. - -However, for simulating a galaxy you may find it more intuitive to define the elliptical geometry using the -axis-ratio of the profile (axis_ratio = semi-major axis / semi-minor axis = b/a) and position angle, where angle is -in degrees and defined counter clockwise from the positive x-axis. - -We can use the `convert` module to determine the elliptical components from the axis-ratio and angle. -""" -galaxy = ag.Galaxy( - redshift=0.5, - bulge=ag.lp.Sersic( - centre=(0.0, 0.0), - ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), - intensity=1.0, - effective_radius=0.6, - sersic_index=3.0, - ), - disk=ag.lp.Exponential( - centre=(0.0, 0.0), - ell_comps=ag.convert.ell_comps_from(axis_ratio=0.7, angle=30.0), - intensity=0.5, - effective_radius=1.6, - ), -) - -""" -Use these galaxies to generate the image for the simulated `Imaging` dataset. -""" -galaxies = ag.Galaxies(galaxies=[galaxy]) -aplt.plot_array(array=galaxies.image_2d_from(grid=grid), title="Image") - -""" -Pass the simulator galaxies, which creates the image which is simulated as an imaging dataset. -""" -dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid) - -""" -Plot the simulated `Imaging` dataset before outputting it to fits. -""" -aplt.subplot_imaging_dataset(dataset=dataset) - -""" -__Output__ - -Output the simulated dataset to the dataset path as .fits files. -""" -aplt.fits_imaging( - dataset=dataset, - data_path=dataset_path / "data.fits", - psf_path=dataset_path / "psf.fits", - noise_map_path=dataset_path / "noise_map.fits", - overwrite=True, -) - -""" -__Visualize__ - -Output a subplot of the simulated dataset, the image and the galaxies quantities to the dataset path as .png files. -""" -aplt.subplot_imaging_dataset( - dataset=dataset, output_path=dataset_path, output_format="png" -) -aplt.plot_array( - array=dataset.data, title="Data", output_path=dataset_path, output_format="png" -) -aplt.subplot_galaxies( - galaxies=galaxies, grid=grid, output_path=dataset_path, output_format="png" -) - -""" -__Plane Output__ - -Save the `Galaxies` in the dataset folder as a .json file, ensuring the true light profiles and galaxies -are safely stored and available to check how the dataset was simulated in the future. - -This can be loaded via the method `galaxies = ag.from_json()`. -""" -ag.output_to_json( - obj=galaxies, - file_path=Path(dataset_path, "galaxies.json"), -) - -""" -The dataset can be viewed in the folder `autogalaxy_workspace/imaging/simple`. - -__JAX Variant__ - -For an order-of-magnitude speedup on large or repeated simulations -(parameter sweeps, mock-data studies, batch figure generation), construct -the simulator with `use_jax=True` and wrap your call in `@jax.jit`. The -simulator handles pytree registration internally. - -```python -import jax - -simulator_jax = ag.SimulatorImaging( - exposure_time=300.0, - psf=psf, - background_sky_level=0.1, - add_poisson_noise_to_data=True, - use_jax=True, -) - -@jax.jit -def simulate(galaxies): - return simulator_jax.via_galaxies_from(galaxies=galaxies, grid=grid) - -dataset_jax = simulate(galaxies) # Imaging with jax.Array data -``` - -The `dataset_jax.data.array` is a `jax.Array`; `aplt.fits_imaging` and the -plotters call `numpy.asarray()` internally, so saving / plotting works -without manual conversion. - -Note: eager `simulator_jax.via_galaxies_from(galaxies, grid)` (no `@jax.jit`) -already runs on JAX and is sufficient for one-off simulations. The -`@jax.jit` wrap is only beneficial when you call the function many times. - -See `scripts/guides/api/data_structures.py` for the broader "JIT-it- -yourself" pattern. -""" +""" +Simulator: Sersic + Exp +======================= + +This script simulates `Imaging` of a galaxy using light profiles where: + + - The galaxy's bulge is an `Sersic`. + - The galaxy's disk is an `Exponential`. + - A faint extra galaxy is included offset from the main galaxy, whose emission must be removed via noise scaling + (a `mask_extra_galaxies.fits` covering it is written below). + +__Contents__ + +- **Dataset Paths:** Defining the output path for the simulated dataset. +- **Grid:** Setting up the 2D grid of coordinates for image evaluation. +- **Over Sampling:** Applying adaptive over-sampling for accurate image simulation. +- **Galaxies:** Defining the galaxy with Sersic bulge and Exponential disk light profiles. +- **Output:** Saving the simulated dataset to FITS files. +- **Visualize:** Outputting subplot and image visualizations as PNG files. +- **Plane Output:** Saving the Galaxies object as a JSON file for future reference. +""" + +# from autoconf import setup_notebook; setup_notebook() + +from pathlib import Path +import autogalaxy as ag +import autogalaxy.plot as aplt + +""" +__Dataset Paths__ + +The `dataset_type` describes the type of data being simulated and `dataset_name` gives it a descriptive name. They define the folder the dataset is output to on your hard-disk: + + - The image will be output to `/autogalaxy_workspace/dataset/dataset_type/dataset_name/image.fits`. + - The noise-map will be output to `/autogalaxy_workspace/dataset/dataset_type/dataset_name/noise_map.fits`. + - The psf will be output to `/autogalaxy_workspace/dataset/dataset_type/dataset_name/psf.fits`. +""" +dataset_type = "imaging" +dataset_name = "simple" + +""" +The path where the dataset will be output. + +In this example, this is: `/autogalaxy_workspace/dataset/imaging/simple` +""" +dataset_path = Path("dataset", dataset_type, dataset_name) + +""" +__Grid__ + +Define the 2d grid of (y,x) coordinates that the galaxy images are evaluated and therefore simulated on, via +the inputs: + + - `shape_native`: The (y_pixels, x_pixels) 2D shape of the grid defining the shape of the data that is simulated. + - `pixel_scales`: The arc-second to pixel conversion factor of the grid and data. +""" +grid = ag.Grid2D.uniform( + shape_native=(100, 100), + pixel_scales=0.1, +) + +""" +__Extra Galaxy Centre__ + +This `simple` dataset deliberately includes a faint extra galaxy offset from the main galaxy, so that the modeling +examples can demonstrate the `__Extra Galaxies Noise Scaling__` step end-to-end. Its centre is defined here so it +can be reused for over-sampling, the galaxy itself and the `mask_extra_galaxies.fits` written further down. + +It is placed inside the 3.0" modeling mask, towards the edge where the main galaxy's emission is faint. +""" +extra_galaxy_centre = (2.2, 1.6) + +""" +__Over Sampling__ + +Over sampling is a numerical technique where the images of light profiles and galaxies are evaluated +on a higher resolution grid than the image data to ensure the calculation is accurate. + +An adaptive oversampling scheme is implemented, evaluating the central regions at (0.0", 0.0") of the light profile at a +resolution of 32x32, transitioning to 8x8 in intermediate areas, and 2x2 in the outskirts. This ensures precise and +accurate image simulation while focusing computational resources on the bright regions that demand higher oversampling. + +Once you are more experienced, you should read up on over-sampling in more detail via +the `autogalaxy_workspace/*/guides/over_sampling.ipynb` notebook. +""" +over_sample_size = ag.util.over_sample.over_sample_size_via_radial_bins_from( + grid=grid, + sub_size_list=[32, 8, 2], + radial_list=[0.3, 0.6], + centre_list=[(0.0, 0.0), extra_galaxy_centre], +) + +grid = grid.apply_over_sampling(over_sample_size=over_sample_size) + +""" +Simulate a simple Gaussian PSF for the image. +""" +psf = ag.Convolver.from_gaussian( + shape_native=(11, 11), sigma=0.1, pixel_scales=grid.pixel_scales +) + +""" +To simulate the `Imaging` dataset we first create a simulator, which defines the exposure time, background sky, +noise levels and psf of the dataset that is simulated. +""" +simulator = ag.SimulatorImaging( + exposure_time=300.0, + psf=psf, + background_sky_level=0.1, + add_poisson_noise_to_data=True, +) + +""" +__Galaxies__ + +Setup the galaxy with a bulge (elliptical Sersic) and disk (elliptical exponential) for this simulation. + +For modeling, defining ellipticity in terms of the `ell_comps` improves the model-fitting procedure. + +However, for simulating a galaxy you may find it more intuitive to define the elliptical geometry using the +axis-ratio of the profile (axis_ratio = semi-major axis / semi-minor axis = b/a) and position angle, where angle is +in degrees and defined counter clockwise from the positive x-axis. + +We can use the `convert` module to determine the elliptical components from the axis-ratio and angle. +""" +galaxy = ag.Galaxy( + redshift=0.5, + bulge=ag.lp.Sersic( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.9, angle=45.0), + intensity=1.0, + effective_radius=0.6, + sersic_index=3.0, + ), + disk=ag.lp.Exponential( + centre=(0.0, 0.0), + ell_comps=ag.convert.ell_comps_from(axis_ratio=0.7, angle=30.0), + intensity=0.5, + effective_radius=1.6, + ), +) + +""" +A single faint extra galaxy offset from the main galaxy, representing a nearby contaminating object whose emission +is not associated with the target. Its emission is removed in the modeling examples via the +`__Extra Galaxies Noise Scaling__` step using the `mask_extra_galaxies.fits` written below. +""" +extra_galaxy = ag.Galaxy( + redshift=0.5, + light=ag.lp.ExponentialSph( + centre=extra_galaxy_centre, intensity=1.0, effective_radius=0.3 + ), +) + +""" +Use these galaxies to generate the image for the simulated `Imaging` dataset. +""" +galaxies = ag.Galaxies(galaxies=[galaxy, extra_galaxy]) +aplt.plot_array(array=galaxies.image_2d_from(grid=grid), title="Image") + +""" +Pass the simulator galaxies, which creates the image which is simulated as an imaging dataset. +""" +dataset = simulator.via_galaxies_from(galaxies=galaxies, grid=grid) + +""" +Plot the simulated `Imaging` dataset before outputting it to fits. +""" +aplt.subplot_imaging_dataset(dataset=dataset) + +""" +__Output__ + +Output the simulated dataset to the dataset path as .fits files. +""" +aplt.fits_imaging( + dataset=dataset, + data_path=dataset_path / "data.fits", + psf_path=dataset_path / "psf.fits", + noise_map_path=dataset_path / "noise_map.fits", + overwrite=True, +) + +""" +__Mask Extra Galaxies__ + +Build and output a `mask_extra_galaxies.fits` covering the extra galaxy, so the modeling examples +(`imaging/modeling.py`, `imaging/fit.py`, `imaging/likelihood_function.py`) can load it directly and apply +noise scaling without a separate data-preparation step. + +The circle is sized to ~3x the galaxy's `effective_radius`, derived from the same `extra_galaxy_centre` defined +above so it stays in sync with any future tweak. +""" +mask_extra_galaxies = ag.Mask2D.circular( + shape_native=dataset.shape_native, + pixel_scales=dataset.pixel_scales, + centre=extra_galaxy_centre, + radius=3.0 * 0.3, + invert=True, # `True` inside the circle, i.e. the region whose noise is scaled. +) + +aplt.fits_array( + array=mask_extra_galaxies, + file_path=dataset_path / "mask_extra_galaxies.fits", + overwrite=True, +) + +""" +__Visualize__ + +Output a subplot of the simulated dataset, the image and the galaxies quantities to the dataset path as .png files. +""" +aplt.subplot_imaging_dataset( + dataset=dataset, output_path=dataset_path, output_format="png" +) +aplt.plot_array( + array=dataset.data, title="Data", output_path=dataset_path, output_format="png" +) +aplt.subplot_galaxies( + galaxies=galaxies, grid=grid, output_path=dataset_path, output_format="png" +) + +""" +__Plane Output__ + +Save the `Galaxies` in the dataset folder as a .json file, ensuring the true light profiles and galaxies +are safely stored and available to check how the dataset was simulated in the future. + +This can be loaded via the method `galaxies = ag.from_json()`. +""" +ag.output_to_json( + obj=galaxies, + file_path=Path(dataset_path, "galaxies.json"), +) + +""" +The dataset can be viewed in the folder `autogalaxy_workspace/imaging/simple`. + +__JAX Variant__ + +For an order-of-magnitude speedup on large or repeated simulations +(parameter sweeps, mock-data studies, batch figure generation), construct +the simulator with `use_jax=True` and wrap your call in `@jax.jit`. The +simulator handles pytree registration internally. + +```python +import jax + +simulator_jax = ag.SimulatorImaging( + exposure_time=300.0, + psf=psf, + background_sky_level=0.1, + add_poisson_noise_to_data=True, + use_jax=True, +) + +@jax.jit +def simulate(galaxies): + return simulator_jax.via_galaxies_from(galaxies=galaxies, grid=grid) + +dataset_jax = simulate(galaxies) # Imaging with jax.Array data +``` + +The `dataset_jax.data.array` is a `jax.Array`; `aplt.fits_imaging` and the +plotters call `numpy.asarray()` internally, so saving / plotting works +without manual conversion. + +Note: eager `simulator_jax.via_galaxies_from(galaxies, grid)` (no `@jax.jit`) +already runs on JAX and is sufficient for one-off simulations. The +`@jax.jit` wrap is only beneficial when you call the function many times. + +See `scripts/guides/api/data_structures.py` for the broader "JIT-it- +yourself" pattern. +"""