From cc5a287324a9c90f6a9a49b492ab2c10c9b16e4b Mon Sep 17 00:00:00 2001 From: Jeremy Manning Date: Mon, 26 Oct 2026 21:54:05 -0400 Subject: [PATCH] Task 17: platform-independent whole-fit check; derived bound for a vanished weight CI on bb41822 failed off macOS-Accelerate: 1. test_fits_match_3298278 compared against fits saved on one Mac. On OpenBLAS/MKL the rounding differs, and case 23 (10^3 grid, K=6, noise 0.5) lands on another optimum: on the Mac both versions end at a local optimum (centers 1.75 lengths from the truth); on Linux both recover the truth. The test now fits each case with the current kernels and with the same TFA on the frozen 3298278 kernels (tests/reference_kernels_3298278.py; tfa.py is unchanged since 3298278), on the machine that runs it, and compares only cases that are well-conditioned by a computed criterion: the oracle fit recovers the truth to 0.2 lengths and its fit on 1e-12-perturbed data agrees within the tolerances. At least 30 of 36 must qualify (macOS 35, Linux/OpenBLAS 36). Tolerances unchanged. The saved fits and their script are removed. 2. test_ols_fit_with_a_factor_that_vanishes_finishes required the vanished factor's weight to be <= 1e-8 of the largest; macOS runners with Python 3.10/3.11 gave 1.2e-8. Not a bug: the eigensolver mixes the zero column's direction into the eigenvector of the smallest kept eigenvalue by about eps x the kept condition number (6.3e7 here), so the bound is (K_kept - 1) eps cond (2.8e-8), computed in the test from the gram. Reproduced in docker python:3.10 (scipy-openblas): the old test fails on case 23; both tests pass after the change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Ue9wuWfgA8yRXEuEiydTVh --- docs/changelog.rst | 4 +- docs/validation.rst | 26 +++-- scripts/make_fit_references.py | 35 ------- tests/data/fits_3298278.npz | Bin 79640 -> 0 bytes tests/test_fit_equivalence.py | 174 +++++++++++++++++++-------------- tests/test_tfa.py | 22 ++++- 6 files changed, 137 insertions(+), 124 deletions(-) delete mode 100644 scripts/make_fit_references.py delete mode 100644 tests/data/fits_3298278.npz diff --git a/docs/changelog.rst b/docs/changelog.rst index 6ee6938..d41cea0 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -36,8 +36,8 @@ Added: ``||X||**2 - ``. A cost evaluation is about 3 times and an accepted solver step 1.4 times faster on development_fmri (K = 50 and 100); whole fits 1.1 to 1.5 times (see :doc:`performance`). - Checked against the earlier formulas over 492 generated scenarios and 36 - saved synthetic fits (see :doc:`validation`). + Checked against the earlier formulas over 492 generated scenarios and 35 + well-conditioned synthetic fits (see :doc:`validation`). - An ``'ols'`` fit whose width shrank (with a small ``lower_ratio``) until a factor was 0 at every voxel raised ``LinAlgError`` in its final weights; ``ridge_weights`` now falls back to the minimum-norm weights the varpro diff --git a/docs/validation.rst b/docs/validation.rst index 0513dd5..7c6e6b2 100644 --- a/docs/validation.rst +++ b/docs/validation.rst @@ -185,16 +185,22 @@ keyed on the widths only failed 32 of its 40 cases, one keyed on the centers only 32, and a cost cached on the point alone (ignoring prior changes) 27. **Whole fits, well-conditioned synthetic problems** -(``tests/test_fit_equivalence.py``). 36 fits with the defaults (2D grids 20 x -20, 30 x 30, 25 x 15 and 3D grids 10³, 12³, 14 x 10 x 8; K = 2, 4, 6; noise 0.1 -and 0.5) against the saved fits of ``3298278`` -(``tests/data/fits_3298278.npz``, written by ``scripts/make_fit_references.py``). -35 reference fits recover the true centers to 0.2 lengths; one ends at a local -optimum in both versions. Tolerances are those of the stopping rule -(``threshold`` = 1e-3: centers in lengths ``sqrt(max_sigma / 2)``, widths and -weights relative) and 1e-8 for the cost. - -.. list-table:: Whole-fit differences, 36 fits +(``tests/test_fit_equivalence.py``). 36 candidate problems (2D grids 20 x 20, +30 x 30, 25 x 15 and 3D grids 10³, 12³, 14 x 10 x 8; K = 2, 4, 6; noise 0.1 +and 0.5), fitted with the defaults by the current code and by the same +``TFA`` running on the 3298278 kernels, both on the machine that runs the test +(saved outputs would carry one platform's BLAS rounding). A problem counts as +well-conditioned when the 3298278 fit recovers the true centers to 0.2 +lengths ``sqrt(max_sigma / 2)`` and is stable: its fit on the data changed by +a relative 1e-12 agrees with it within the tolerances below. Only those are +compared; at least 30 must qualify. Tolerances are those of the stopping rule +(``threshold`` = 1e-3: centers in lengths, widths and weights relative) and +1e-8 for the cost. On macOS (Accelerate) 35 qualified (the one excluded +ends at a local optimum, centers 1.75 lengths from the truth); on Linux +(OpenBLAS, Python 3.10) all 36, with differences of at most 1.5e-9 (centers), +6.2e-8 (widths), 2.8e-8 (weights) and 6.3e-11 (cost). + +.. list-table:: Whole-fit differences, 35 well-conditioned fits (macOS) :header-rows: 1 * - quantity diff --git a/scripts/make_fit_references.py b/scripts/make_fit_references.py deleted file mode 100644 index 8bfdf0c..0000000 --- a/scripts/make_fit_references.py +++ /dev/null @@ -1,35 +0,0 @@ -"""Write tests/data/fits_3298278.npz: the fits of tests/test_fit_equivalence.py's -FIT_CASES by the code that is importable as ``htfa`` (run it with the code of -commit 3298278 first on the path). - - git archive 3298278 htfa | tar -x -C /tmp/htfa_3298278 - PYTHONPATH=/tmp/htfa_3298278 python scripts/make_fit_references.py 3298278 -""" - -import datetime -import os -import platform -import sys - -import numpy as np - -here = os.path.dirname(os.path.abspath(__file__)) -sys.path.insert(0, os.path.join(here, "..", "tests")) -import test_fit_equivalence as fe # noqa: E402 - -commit = sys.argv[1] -arrays = {} -for i, case in enumerate(fe.FIT_CASES): - for k, v in fe.fit_case(case).items(): - if k in ("centers", "widths", "weights", "cost"): - arrays[f"case{i}_{k}"] = v -arrays["provenance"] = np.array( - f"commit {commit}; TFA defaults, FIT_CASES of tests/test_fit_equivalence.py; numpy " - f"{np.__version__}; {platform.platform()}; " - f"{datetime.datetime.now(datetime.timezone.utc).isoformat(timespec='seconds')}" -) -os.makedirs(os.path.join(here, "..", "tests", "data"), exist_ok=True) -np.savez_compressed( - os.path.join(here, "..", "tests", "data", f"fits_{commit}.npz"), **arrays -) -print("wrote", len(fe.FIT_CASES), "fits;", arrays["provenance"]) diff --git a/tests/data/fits_3298278.npz b/tests/data/fits_3298278.npz deleted file mode 100644 index df8c89848b011763410a732b2960ebedd21eb4a1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 79640 zcmce;Rajn28Z?>!!QGwU5ZoPtyIXJw?(pFr+}+(FxVyW%ySux?$^Xx0&+N?EH)m$f zdf@3BF6!;YTUD#Ox@E+HL8v}__yGRv7yrWoi>t}o9}mI@g>zc8CQ(&@a;~Q<4^?nxdX#k(*%=ZID@-oskx&r63oR zfR?Tm69tLrAU;ek^(hF=kO(yg=mArR9a1QX9kK>e54bed^(?;dCDw|}zID|_=10%Q zlL43g0x&4-gX_VYba?^WhYvD;ftuRrMf?sG{!ge5#(K6!e&iJFJa2;lOuT|C~VwePcr-+dpxA0R9J{U7&xAsr@euQy;P0 zFWdrqdQRAtSzyF1j&vHLE@e)xS_lESe7OL+Arj$jX4$yZCR+`=(o3rnLiGmRxY$xP zTjBJ;06Ak;`GlAeg=7`Al;TG;=eJvT#crICP#!)RNHWCQuVhka75s_E8W}#p4!L(;>wkQq3^#<4x#b>Y=NZbmuM0-=+H$4rkb`J%! z#%kE6Yx4NXUq3#5ghuj%nl3iQwRdFSCio&cr?Z|#M6nO#n^Ii_E}wwD&AJs(YIqxb98jYjoU8Tbu@z^EB z&a^sk+Aui2cI|(YG4}2_tAk;az>zw#EH@>YUZ2gJH-S(1O41wp^i)0HYM{;XIQsp< z0c}nLRx6Y4=K3}*$6-Hho7~J`U#9r+n|5)M7qSC*HtnPJTC}g>Cw7=-Iz_6R9F_9w zYX^D3){moZ1y`{a3|lPDGFuVffecI1!$sX?6fP9rQYoWR7$dhF($X3}nIuWC*$uGL zs_VA|zWN-({k4B7!v)3O8wK{y{?)a#v3=*6fIqrc{6Fp9fAsFhQb`U>lGnx7<&VEn zNvoHx%R7|!o~8Rc)-(K1(fkYRQKLZ|RC_Hg!defQP>u#G;y|a@s_Umx9LU^D=4@pB z@C5!1shM+c<2%xK(xcP*e8VMTy7{Jo!djEo98RLgV*iFJVe<#eJF0)q z-`{ZF=I@**TEi)n+EksgBD#jXF)eE0$|=C6g!|=EqII9uCn8^s{JI^2sphbedasz2 zl2qS~lhQXuESOWDp@&3u9yR(h0?DXl*=juz!`amL;dFs#BSp|DF zDxukWqoX2V6SHPMeX^!XnWfDteQ5rV8<{wKE1fzyuW=X=3F8GT8VU|h?E@k1+ z9IHmh%qEH1kHu=f%jHo`^C?w1h4xB1{o-2Tm)*ScxZ1re#)%-J0yw086u`n=qREBJ zg6OraL>cB--kUU^QcPy)h57&`HxR^_=`lt-Sfd+)-sdkaJ>e7W*Xog`-XMYf$p=hL zd+Uw8n2<+KIOmwIB6#6CM(7iv*ETU9(%sVv^AG(sZc!d;a=%;~OokPkJcC}QhJqj> zYrnl_r0-M%xwbWE*_FoVP+x3!O6TiS>NudnBo2RtVzcU3r=-mEs*-_fKQp zI{PdKmPWPG46eGX0fZ?GQ?@01qG7BnqYHn`JcQttVF}aqq=NL=iqLK1U4#NZrBesx@t7@x&mz%$jD4^^r0pg9A!mw6b+cSBx+ z>%|F)JuUj((Yj#5H#xj&SOnUgtx=-Qr1Y-MC4Vh`i?^aTv2I+b}#blY^-%+h@6`Jo3 z`P|KsSpC7|2=A|acT9$k^xh=zeD{lF{+E1*M3G3X>ZX43yGL>5XBl{h^4_oXzw=%D zAN29g!vCLL9?8GByt%e;(}TT2j((Rw5Wi!2+n+-!BhF; ze-L6R`zHv{j=)P;V)edOyCh?_;T5WSH_Jrdt{Yf#xH*=_-(kKahpW89{O64S#*jZg z{R5!*ALB1q$b$rq-EN&c$S;T)o8L23F&I?}1}C|uLnbX~EryX)pwo~^Dx5laM9r10 zo`TEPc*L1g8gNEUjzg8`-z&$rsD!SIQFZj{`}BP4WC-R5Vz$CN&Ylq~rqe@ovLUNp zLdICj;A-wT^TA*Z#6E9yfS^R4GxTJh| z6I$!A`Pq&(R^}QND$VebDzvD_#}dnbI_vtbNXx6&mhCn|dVb8|G8;+KBJtlOvdyP;@fhgDRzVK=AAV;xC#GycT$DA|`nw za}wBGN^7eiK?e@{WklxiXI}jH^-Ab#HXMVL>@ocZQLeLwh9Q5{FMCUDa%Z1k? z7@kI?zTQpJENLQ6c89~%pQuw=PB)zI$OU%7-g@nEJ6&Fxp|YEwF3eRqO7`|#6;fD+ znlWaAt=S-i?4)4GDobVO=<+`o`R-Z3&x0wLUar3cn{1QSLqIT(I3Bue#hjbKARlqn zG-OL>^v(Ioior^60%eI%KjB@yi9@}{{~cWFD-OsEQKZ${zEIO-L(Yp7`R z9+p?p+<$>QA?w^1JA10{k_OXF5%@Z&z_1xivMwg|2_q91FNqn%2DICe zX?;I4JBd3X%S;Vhk0E)jvYEeq=;_1K*UD|7Jj{iZ`9zpk8d~>8hE@AkbQ2r0IsfA1 ziHG7a5j|!vgn~3~?vH%aRWfde5}2O?LpX=1^+72_CD~bOgse5Un-7SNSgg8l2<-+S z*vt-OC=$U;by!D^G;;ajvhGv-w()_xl8l2lZoYHkWL9PBerI+w+skU>E|cRG!$AqM ztvMJ>n*l`=+i?um*X+DG>GCd$zqWezMPKX6DqC4u| zD#)I%3m;!s!xtFR<;PxvhWQ3yKC?6Sfa9t3eWRU01ZQ4#{#?uq#WI+UHF^#jiShbl zNc^a?fYN_OC`#%lyV~+m9Mg~L!){C8G)SVxVStLjRtLn9Sk2RcoI>mJ^y9$+a&PV#cYU8J1KVWJwD@vt0By*a*Ms;@ zTjFmIkJ24Z^)X_45*Y&J1kKCOP{##ERCT)7z)Ox*dp22!QD}~Njf_y|#ef*`NtBKH zqR5p~Z~_l~xsHa>K;Ck=ZwU5JU?`(xr-}5iq7mbc%@yl5U_^cDA_G>Y-AwUh6UPP0 z+AgCL6KAE!r_NfE97=f|Jh}O3`U)+ciaQnj(ZPnFV`8(=4hwy(7e>Pu6t@yu{frk; zXhA>8MLO~7gAJo`22rO)5Y}x6Dv^u(>ZPl&SfsaweIH$(vuHB-oYfz5UfF>6{R~Pj zZ2?Zlp)p>?cl*;b)L*5zV(01dosHk6_?Q0uFQqt!f&MvR8ZtCevIuqZtl&(of}D+Ptn4zn7?aW)j62idnxZy%j^1vqT9RZt&n%oT ziJ>x0P;gE!Jg0cypBcJz1l6}*3V+pAse zKz=)uV%xwueS5K!ysnY3H=CLZns{tOsCRP5T{Z31{^??GQFH%kzRGI|OK6Q(-#u+; zaOR}Ew%|OIwE>MMuqZBwun1B-V5;4aN0S2W|Htw?IiMNf?6}0o@_Q)-*1w2f7~FOp5LVN3&1IJ zv<^a>yaR(_fl_8t3Y!A_(DTY;nB=uO6OL=!Gt5S9(VED{RLms-1RCybc2N*PSAEjm zBzy;*eOC2xjgwjkpf4^zoG2>uBq66A?YvP@B)Pr>LzplrNmEP!9_cZ6BFO#5lyqaS z4x={KJq>Mh7nnsiJ->f<5Mhn?#Y`g9S~=bG=K7+ebflzz5l7U?aM|-x`haRzdP%!v z1TMIZsA9rj-bEO2%a9gSSfI#+(Hr!w)mcx4*e7o(y(%=Gk7?_Ut^*Ty{<Aq_4hH= zkFMwuoKCkj@MjkF^wal_&~?40xN(4n`Gg(Uqj(f=qtB#M8q=#lO?d+hE*G?B(~B{M zsuRA5@iA(s7ku-vQChB*zy8g(uAf=2BqE5e(NZSOe0dU|@oZM~<}>h)^ag7gM%(fBkCIYtaDxI--^>6jfzumAh&-aqsShWdx3u~JMba?L z)#-+)H-A#BOTMSwP*93yZvSCV4u?Ro4U$gHT=DZt1KxERgurCk*~O-IUaa*6=c$ea zh*@Ri`kMY0{pQI>vM-t=TTyM9tCE1You0XG18R@NTwlS7{elyHkX&QPV}c)C@N}Bx z0b?7&u8d-AW9W@M5Lw(sKn}$CB+RQ-E(MIf%oUGz3!hrb1)aob4sXezqXQV@CqE`| zi*sekk^m&#X-q*dSTwkL2CfGQkY_ACnhY9F1hB}evg zYzHF>PLF)t0aRHS^rlnkrG}q}MVxosafzI?4Ge4nFB@I?v4gpDW5^;~XF0a)v$lXo zyRjTvkCg@FRgk=*B7Hu7xBY`|Xq2>weLcEQ{B3SInNBj`ED`UIls3OzmzB} zbnHtyEAt~I^S2O-AL}}x_}}4SJ8QblS~=F0X$?A%M0H)XjHzJ{Jjmd9eTh2Q19_TJCFlfS#71LW zNh?aRQ8M@Gy>rwNsFG7Bd22+CaS9|=iJ}xY2hbvm7||VQV5%5FSg&%1X;R#QS!o_9C_bTWap*nW1?1e9@QF?E?VZSLXHo~ zLl$F}Ph>fu_=EP!&{GUHpEMpu7oI;Ft?buv*YO2FBt@fQptkHE0O-9i+5)gQcaRdL z!7_((7@F;akp-vMEFR)-Sbow}Yo;=-fS;&BGSouw!G}Gdnwz0~?bYA930m1g|CaFu zRqrLZc)gC8ya%uR`_e8(lQhFRgME!}la$ZI)y6PFEPnIe`pKzd38rs~tnN?{eX+(b2P092`c;F*?h1J6k9cKE`W&*C-oSJrWrp@IffCg)&#*8P z9@uhnXq%#LhQ45dC*#zaC6q06!5_^z;zQpewXIoY%lF$q9eR_DjSE|Dj2F_Zq$(ho ztZJOhce_(?$VF*O!=E_gm|_RpZX1VGt=xCmLv|@Sv_~mw3gkKh zefa{63=CCIR@Z;_%}ljh7?^_8Ny*1#Vp``eX4oFt7aNm@Y?qLWStG2@AGrEUr@eR^ z@%_s2zKmS-hZh$2^HFm6PJf+H%sYRA^!9?G5#uS#t*vkKqw^@y6d1k_4uO-qmY3scG(L%8~9|O$RQQXJKsyw{ti8*S& zi;)WS%^Rj&s`}K>_a>!>jqqmlv&ZVUg?MUsX}`K~LmOGHBNm&iXJ5y1gpY7!dTktO z%CSP^HZ^!{oqz-#)^+KjsB!-kWs!`Lyw3Q!@>Pd0@RNiSTvp}?0sD};c$ovvkGx}) zSC_-VbvZzpTZBHA<5*$Pn7;Ld^J9ZX#t5s$*~Bq6V2hAxVq!gc3#e4w=_V$c0|IjK z>moBp z`^Gy8Pj0VxWLJA&la6ttqEfmH%dZ zBT>+GD}aTRBr=MsrVLOgM5o%PYRMrPuO0K*>mjDM2T^_n;%*;yT7}#Fp12B9vDnp( z$|kS8R6_p3WF8upm}qYGFx_kVqI|U0&)v#o#~1|IjD{qOJ5poOum6L0|jqV$8@0yRvLmb&&b-3g> zyh`-#{3%u~*Fe)L30D}eF?QoEEOUg?;x5mb-lk8`+g{%nD6+S)M2``%h6InM<#k^n zgvx*tY!?S^YU7#KmnhnBy-%LGLoh00SKci-uD8(B2XxG4ifNv7o-kG2H85PmfbfFe z0R_v#P!s&b=Y6xpHR5qYFU*Han%(62{aIbE{Rt+YT;uMDtrjY)ah5s}r9c6d%=1-1 z_Z0U&2C`s!zT`?PxCwX)J56p;T|R~(o@wB^>(TDA3B_oO&OP6I2)#D zzaJuu7rww=V&*y_2p2dG$c+}qv4=wM?eC6f882rVC9oW~1lSAr=8yFc^ZF^7=uVh1 z%LG__bjl$w$7KiX0LAg1S;l>Lu7y8QI6Gv4N5E+o8guCBSKgc2nylO;<>i^4k+)bf zwqhPYV1&Ky&<0@NHFT(nuxK8|-y0!sYoZtO6rBm#sO`Q2)uMM5VezpaqRBAfOZ=4D z+?Y(slfTgTc&K&2a^H`EXL0$*Wd!u}0-|_`WghE^jHW$XT6p1Ni zJRt;|0Vg*`8LFIc$+a}=f@9p=lsYUm*=v5*6jsjHay>4o7pX_L68TKJl{^ zt6#G*ZuIdSs>a!t?wFt#fAyJi_0=+fhse{DI{k|xe~^mR?6~?COAYrinKsFAqgBVg zW?4rx6Jm@4#XzA4!seNZu%_d~4Zb&Z1(iASW;`&Pg60X4s!a*SrjIw-F3rQ^XlRDH zm04!x#)lX^oGE*VWUujg71@bE;(MCE8bkaG~T-pE9ATlYX`ZLc`S^d6j z(Chnn{to5cdocZ-K`i=*_4xlDjrmuDxI%+lDB3n_lR#Eo*TSsPlc~Y(IMk{81(tqh z$f6B?1~6}So_BEutpk=HVhk{5+k{n755Fx;rso*`M&9tGk)771coKS3QC9X{O85KJH!WJiqaxnSqH1< z8uu+113Le@HTT&;HzmL~t*eXeNCPMTXxC3PElOWlITaF-5*{_ze&5Sij6ffL2ma63 z^fwl<;N2odAMhXtVg5J?#E#wZS<@H8;2V_MnIXXuAQ!ZYs7M`TPN39>n91}(#jg}# zpP4@#&vA&*8{>{F(Lwh%ICiOfqFLTsx^nf13QMo8tk~Ow+J?f9QXimp7tjT=PEW>{ z_5Nymb$fMZa|^h>C2V*tOJ>s|X?i`>gEPdN(<&`2;~ifq$>XgTy7xd}pa}ZH6u!x_ zlRB1=;~TCcXs;3-&GdHR{lZJmG} zC`M$&O|fk%1;?=iNbOH?Jk1Dn1eI+;ea74y5uusRr})HZXI3KwnC_xESxps63!$BdV;MXh+X?iz*T=6D?zoFoofOpzS@q@$5}I&mWDc(~m~Bv;*Ym6dbUGTwGWG_ElObW|z|e090%!#KNF$T@=2 zsYj6@vT){~v{lKvFDUZp^D#+`W>^`>yBj_~e+pd|wv#e;04kazw{x_@M2>onk_;xX z!e7T;*fK*aj5)jEZ0PS@_E#kLNBJbdf5~yoe-fWDALgN(dJWDQ;JBQ;viMs$tUP;Wn zRuoLbw$1`ea}Bt7;g~GaZEJqT6cr2S4)^u&*}6D~1gt+^%gWA|gw~XIJ~fH!3t{em_5fYxXB(u|COJ*?G8z2LHR-46M?%)`=2gL^%U z*1UZ5kEFT-&PhN}jt_W1E0QjsEYRs!Oh6GzzMo-UgH)EM-I6W#8-4hGdhnW{Z~NKT zRY5Mh5rRysRo@>Gbo?Y$&SLcv)9B=jb9ICwQ1Q_eDW}6V?+YFSR)MECYiqT-9gPKe zr%M8p#DuZY{YJT9FzpB39b#ssu+(~k6~V(+X2_l|2GxYlJ&s25qZGGC5VUmmFA~|z zUH6eeMWV@&p7C;Z4`S!S5oz@|#|6&7VP$L0yzzoPdB^1@saeMJ-(Hx8hZ`YrH0hFj zdfBWd)Fjd#H#yw*+GwsQi(zSCkNb@+cfJpLyA(q{z!py@ghnG_3j&Kvk;tcpOWo4L zFlrI_T0m31mfK774YF|n}Qv|h`Ti=SO&xTG&@bI>vb?+-& z5&2Ld8tu&R80Xut*C%$BsMByDzw>PWk)g@#hjsR-lXYB{5a^<9RyldRyfF>2HV*?x zDu_Aa&=;2r)pCGe+E;>Ox4N=xcMX^vU(AAe6w0n(46b~SD@;#aLj{!UbZ6GSYo07K zE?=ix&?Or}VP2tITMbDDuV9gwJK_q#dCaz8DafzQ7pe*iuhX2%fW!nU4i>-A3PddS zEf+fa?5s-W0KYd(6~p?r;yyC#@FS86ruHYJQbjE-2oLH6JPA7W4zaP=!>py?guWHf z%o1!{=kI!m)rKEqfryfjWAY;r)5ltBNBiEECBt)-Yo4ziu-VJbt_ z-MFbchnQZlP+3$-0(c6;1oI9-seXMBHNMfn-#KXrg}t zth4tw?fEG2MZiTDoK*uULbeL^0#(5$ySwD|#re-blu#>{$gi(&>kB`y1)`=UiYJmL zcL@Ek71XrCV0z=;nZGkeTqRmV=J;6^6l?feyMLehF{Wk3;A(=sgT$0~DDSp|`FD4i@rP0RZ$=19bMz^B zYYXCF!Dh|Q8xZg}WMtwo;M1vg;U9EeR`prT!G253WJe94z9W5icFeyyyZ7Ame=kB< zdR6R7@yRI^JEj_Fg`?*A_Xr`&5~JEXs(;Sk@3!ra*!0%FLY#aZzEw{INknzDF00KX>JzJnMny?;On!bl(FqMx`J}au8TwQ5%hsc8HDgCJe zF=_!T4ZCZzoMp8mR)~AMIbCQJ|IJd(mf=QU0u0!BfkaGMNwFNzNgeIoV)atV#qeDl z^Dc0O_@O=ko6TEVeryx4p;DVh-O~``&+Yv&x;(9&yzu9^=Vce+Se>HBkKUJ;*3VM@ zW!)+>XPuJ!9NV=WR1S<^sb=al!w>oNo%b&H4ltLoj+({l>7YolM2o!@hz;zpO5v~k z>8uD=qa3zz0t(RhM1qH`WTV}`^?u_$(uJH7a^ONPDjZZaP1Y;jy!|mJu9{Mvtq{xV z&;{IEpp795HCFsAsNqFCY+pipJ#xrf$^E@JEg5Vz4Z>z{20~rWq>qyV1)WaZ*tpy= zytM`)^9W5s-%00{x0uPs`NZ)>eq$bFdBzyER*7nzH~%7)$z1HG2t6F$Adc^|1}QgV zVN*m!l)eo`9#{x*CN+SA&=%Um`hb%^VgYvo2Q4D5V{SvC?&l5GWr@;3T{KT-i?EX# zX;%OY;+yA<6Vd{bK@CY7T@c0`NMd|a^?KoCZR!&jP|u`SaXY_|BjK|(XN`K;H2bDF z_W|kj8bp}V{7|{vfktwiLHo28O>vQ+bj002_^ORI@+WS#WXsY04wrK3{?<0Il@I}8 zMQbESyKgv(;!L!$;qqI#py*L7eUR$7{9c{L<6h*w>shb^2ASAmIRYlmPe99JG;~`6 zhZEVO1Yv^An;Yi7UIfC3srIX)c3Q94UC$-S0Kt{<;iVRHmxM%#H$XkYYM8c%1is$# zOMq*EbZxnE;KT-lyQvf5&+-ta8-f~b6>y#dW|jqbL;5+14( zl*D_Nyc6CpmiZqN9(5*kq&^CRJ;(2aNBa4Dz&n)phGqFX!uxNY6rslL+rFM;V~0g6 zJ8j`~$jb0$x{E184Qhn-(B;9>^PvB=Ve!{lgPi>FCpW({9?Rb_-hc6=$P^3C=~8|I z$FSKHCpnh+>F=twc#TzJHlG7qMIafd5aU-aYkRZl5n zAC}q8m>J$FS^(mRQm+Sy9MmQ867m}$#t|^s)3s&PO3ZDxxd{h7@}fz?L+9tKb~1@& z7w>}p;A1k70C}#bSdGrn;&)qH3W5|Y8e+|JBXUabBPL(jn(-4X=q9O?V*`d;a|w5K zec7|(hBSXi52b&T*>0}s4Wuguo)=ahz8_Tg8Nlc3k4grI!3i+-FYRi&zlnl2#1-gdk)44EUI+j#(P_D_vyDBZt zcfh1VYgqIoC^bRogEk8?D|YD!xlm6*;rdV)9e%Q#Qy6(0wcH-A%; zkw++J|Nf9yQkdA&#T`r`J_a#DRx#6clHSMl>JSB?0kS{IBkN%{Jja@)nJKQc5FO?X4;HZUV#nz#4}Xc-qH_#<4T&YGi0? zJoN+YVxdQ4b5Z*ZWJWdf@p3uJ_7Z5}yli~Kll4bvP2XVu&FGP?mD9aLdFMUW-+8b752F|O$M_#4j{kW=uBExCCVo`L%IK}W zNn}Q(-_s6iX+a7U$Df_B#?UU|& zqgg(xo4d+-;${hdQO!ZhSu=V~5orTo1pQ*J+s=-eDh9=hXie#bv+F|CnSPaSsqZfY zXm!$q=SaCap)xyASW@JfD5#izOV-(^Nxp#sCi!Yb*W9vbjNP}*`EuH-Q>}Uaxx;n~ z#DGr5>O~N&6#gi6c%C7w+IMa#I4u%GNO;+utAG^f&XB8VbVmjqZq9Tc8q5cY!1MU0 z6k*CPYVEE#ZkG5=sSX`Kv^*&E>1)kPJFEij<)t%z)uy~ULO#}+zu}pb-S~Qg#P^EA zMi?cF$9eS3Hw?Ke6l8-gVOFDjNIlSk3U-SDjmQDRdD1IP+N$GDstUXumQjgp0}=Gp zY?7~z=8+4*{96h3&#-GOF4(6S#PTN0jC=>h_3VhjPVIR)5OGT;+dt5)O)85NjjNCKyk2F*p=brC8OpZ_*|c zg`g13e zlh3WAjTnVNvKcz`T;d$2*wGrAt#n!?9~)bG!pOGjv==ec$`NvC<7DE|SLgh^Q&o9f zydYDX6=h=s-6C%H%ceM98rMhH!{X>`K=Ro+U^YjX?qit7^rehQa}FWCaalY#PRTBF znwNqc=EqY|&+pUo#TYU%aGA=~RO6YqDSvj@({(7rh#sG1{i#bk4d-LsFR3_%VMbGQ z$m(JJ;U}jl^#VW4_pH@?y<|u}NNV(!rAcU;AT86YDgE`v34@!FAyA(pq7^FsS@6zf zmnuzZ7XGq1i!?U_3Qpi?uUdI7*4YOpRb4&;`qWv) z;QvsMsdKmMlya-;*nfA?AE54U-=VyV1KZ!}anOHJkEs_HRlY1NQR0L-aTh4i&DqR5 zoy`r^%16P9-+C}Mv+=GB>MUd<<%ym_W7~sKv9>A5jQ%dj z>KgHXE@-_AGTYw>@_#QvP3?zP`yy-ugEL#bK*>{}w>RFRA#1d!Vd`(ZHa3%nNNe}| zq{}l}{b?NjIq83+$o%h$9DeXCLXG`HL$=F40(D52fFp`s5~7x8$P3pGbDFcnYMN zB(GZFB&5y7*9@9E^J~jd(fCjqZ3j$z%@sdJfxX!mfhou~NvRg_JKOimF=9P@NWbwSdg}J|wWkdIOefpokN0#h64d3mF8qQPfI-EdtMx%+!4@IlOYzq(eQJ>V zdA+k5^qwQil4)^R8&>$fk;4~SsNh6WA|zkO_Q#7|S-|aQT!vP}vU7cqAVZLqP(K4} zNrx?ylP2bn8Y_j%W0f9>3p>7#`q)~uP=!+{DhW6FZOQdqV1Ukap9o_+kyR@%TMI>^MQU4LY;x@TT-LDsEN=W-x9 zpIGV6+7VnqJq>W*i`uUOH!GkG`|d$JZ&spJrt};Xdt%2nfFaRL5bw0mZik(8V=D(q z;Aj9thvO?znaI>UF!P=iosj?A^~P8F|}vBMB5Z2e!eu-5!`ui{0#k zs|Voq4Iy(DC+<*zN=Y-G1DNF%dR#zMdY4x8>9Ii^CI8iZqf+qX7n9-{MOq>O34vJA zQCVGy<1d^VV)$hl#%+AMs-N&a9W|~+<@DkA1`k4AWoEiq7O00|VBLWhDAZdHXNPM# zNUA=55kl2M3SMP~W%37F;iU5S6i`J;y~^!~bbb5qN&?mpW$8D~bLZJ|)3hmkO^zKG zP=3@T1A&VY`D-~!UPQ#DWLl5oBl#xiSnVF*>&cG8UQhv*W4mN!o|hmPkXisYiC6z~ z;(c0Wzko}q>{P>sbUS&%F@$AuN5!H9%+kK*`H3yUH}G_y-Z;!;kX$>ml%{!n8*}7O z=slMQ7~&a<)N{%3%@F1W$S_0C*5AaHB1##x`3%*kFOXL$1)>x}38{xY0Zm?@yvH;V z?%W#uX!^ic+JqLN;I*gh%V_lJOu^oXWHh3peVHFQbSf<0^0&?m2k^dl%f>;F4!EO2 zEI`XM=y0dm$7!H8tpqC_g~G`tT%D^f92UyJ!$9Xsv1y1N{aCZu=P)BcqZG%R$fB0? ziqnj!ymh4TK?nT;TRLU+b8%T{lMqefNV&zEgGWbabLe=nlCGH~mqupHpgr;2bk%d> zN=75SXp^+rr+`T~@7p%{&l$Ge*rQ8m>DSlUN|7owN4eXeqyq2ytb}A#!DC+O6@{Rt z6=iTq2^lQzi%&maWS>r5qzP(zLS5py15lrc2XUFHbO0urnK)|YGktcdZv{8dtpqN# zkd~4isVYn^Q(0pXmHxVVuW^BF!i9+)qMal#!Vn9BWvHI%dQa9aSV&XRcW^J=`QMMt zEZk)-7y=eh*vXYR{N8~6s=}0{-dOL__pZXfr0;*I!lmL3YFc1DXa6P$ul&+F?@$o` zxXgj>-!6_>|3iNNtMTdOhaxrQW=02AOTa_pn+CCeE%X%L!KfE3=1M2mqaO*Z^feE+ zsb;}exlYe@>X%v}SQDhhRb|OJ>-#bEOfuWA9r&ZK;X?KeiNJ3gwVlg>2C%8ah8DKN zsMc0TCbbN!W&gHv*A^j#jPEy z!*kY1aT7DFOU>GL%JcfR)#GAUQ>$|vEELVNSyi&p(eJ`7zz_Q8$+v%bProVm|Gf<8 z*So5IkHU?Mj%;Vs#pCb;);A^reF7))e9HbYhwS6!A@q$RkOWv*W3| z554GoCi3y++v6>E;AI6EYuW=`q?W41|hhR>LyT##C!Y0sRc^tpvmYj4akXQGIS5z zNgiR|DbAohd)MHHRq+_LLf&Kw!`ry9R9@FkpjW=(==_V;mVnrK&Q>LkwxYHdio?HY{E(M)r1EuD-@N2 z`yn*878V*g4p;5`ushfWB?9>zVP=SM*+@Z_e&NzPD4zwVaPrEDXph`1O^TAdM6H7r z@3k^}416%jm`9=P#phNBdo>Fi%veC#dL)y8OvKh8?;NarLlRJJd|`?VZGh7s7Oma9 zpdry#!u#1WdDZGmuS~!sM!0I1B4#;HS3#m+akm|Q5Jumfl`P{z@ssr&Bq|MCTEh2{ zw|$wbww&L>I?_48+lnWmLmfisu9hi5$On zGL!q#aD?(l2uwFL-wxVXFZrs;kAmKoeqvR=6mPLPv}Si+wnQBh0)Zvdz;1x zjm@Rn+y`JSUi&bIC07cCfWQf$hk8bH&UE9mUdQZi%uV~Q%7Vpo-d*6Un7KJ+{Y2L_ znDa_5q{CuQf*IRiD7Pw#Vse2I^(o%2b$vN4!4wJpnJ78D)fu{tc9a(BGoC?TWa?xN z6tx8k8!jxpVT3BaQ@F_TnMfdM`~>tfOG1~xATE1{|1?+xtEFf&>^dV z9!`^ON>BWoFxwrq9omRV^08P;NbT3Sd=W@H4Ck=MA#Bm)&ZkQdo-upwSxx2JyV=%& zY5vq(OavSBil6Kj+e3nC`DQm6&L2xP5OyOayD!!ES*c@6brRVadX-F_4`t3R|%H3l?7Nj~2ckeb<>wnKR*tJg_FSZCQZ zJ1FMC%Z|=EGPZ%CjwpqkonwHgJ_Nf}(+!>8ClrcXhI<0pvUV(N%uzJ*O@WlP2gzWFD~ z7n~OdLzu57Llb#9Oo8=kd9rfcnV=|%udo%*Y-DtG<0?sIf}vO9VRD&r+j3{jExq2H z@N5@g@>FggR*8pHcbiZ%ZnU0<-&pOvqpfYloKA&{Uo~B)By6y2YcRhB7eSrRm@hec z$$S_~_GJ%J5C9y5bHl0x<9q#NBOvwkF~P43p`=)VR#LWz>aAWz!7z5`?yvg52)jN_|? zyZ&e3y-s%1-s*t!1>oQ64ao}GY$?`gO5FCvh7jH{kfsKVt0-!UJ|+M*=yGG zza6wY0cAXScPw!K;?MqW0_sA@_aAQv;cgf?^z$66z^)MknPjf?o^7OcD? z_GNZjHhG1?c^&J?J0}|(37pe2cIFMOq+I^P5{6(qC=?tM{!t>D7s{>QE!Qz3Q1pAb zf11VL@9cUvi~q&LcB~6!)W zwYoC94Gl-azyDiMhy=T&=pFdKyr;kMiuvzeG3Hm48WPJ-5L?Y3Q$SW~a9)leHdoje z9#OPb3|cUlgvXs2*7BzjuNZhLi*O+Otk^bwFBzon;OMb^D25NO<(`&fndRNh$tz&R z!QKIIp^lGAjVcQLS>u)zX@oHl)w3xOG^x8ic$>Cz%{Llt;A4(QOUYWpP>$B?B}SXf zWLb%>NrrKy#@ts1v8y-1S@n2i^8lEy^7wON>Iq#Ocz2}}3;jJBq+?|jG9{nHT53vq z@9)He0$K9!&-0Q1$>S6%MX0gta;;7un2{|a zlts5PM=pz!>Hi0u-QC?GxVyW%1$TFMcXxMpC%Eg$Uh5^- zT+F?`oIU6G1^uZS)myJ#t?tTawq)XY?t?C6a5Y&oO%>GBp^T4Fsho zYb%4RR&$Erpd4+D7HVKWuD4m;<*+#xpMj! zjmf{WX6)l4diQbG1w(bxI`++3X2^#7fjKO@KBUeCxLifd%lk{>Doe#HV|8>FBWPSP zMCDpi5ixafw!jiX>)+}sT>YTVv31VM0;jC0HEb!p{|rJCll6@d_vjY}Imc}V4us4a z2un75!&=rXjUxXimo3I(%GM`uR$sL76nQ07WFu0^Im}M%>N4lo09+>Y%p#mtwm@_* zM19`6PiA{_vyM~{6-jo4&1>e6CIa~B*jb-DS0f-HW4x75%pF&e5$PbRV$Y^l^N_ui zJyc#+968Mi?TVjqnE@aR@v?^2@$Qi)77z>>0+5b~5YH6r-CtUi@$m8WdY6STA}brf zs;MZ@f>c>Mw^xSdh|57TW~zwFSfT@d(PwX+`ip;S@B*{Ll9hfn5?f{V&Y7)b5*YEu zK`2R&e<)j~3CUSN2=COGrgoi7GnX*xYQ9V@a){;hRCk<6lu3lOG%`_hx5PzQ-U(Z@ zdvMK=QJLNE)sMu(c9yeoW8{mvI?R+Xtw?^;hT}Q}a<9xEy(j`*h>ufB+hTxV(AS&| zL%6q;GsaDp`b< zpC~gcp|95<{jz3$!wqjxb^gn^iCSQYAps>0TEH4B;m9LccY-0Ig5rD+AQ9uVoTwAp z&4KuejvPz^pNW9b%5lF;rd6y4ESA2Q4gN7h1vd?Gxq*o2h%rioJ(B5$*15xeALkR( z#M^?iuj;H9ZnCPRLqxsLR9|Ym{}unntjpja6Q3ukOumem9!Np{`On6kc`xNQ{H3sL z3d0y3O^2H`_>M+EzcD28>-O2TvQ>6ic7vEuQPsDuyU}p`;lAxJu2u20b^c))7nZfd zVa%p!?tI`60sbn>F*b-V_n(4QaI&!iD#*L-r@B)js)S*{1}n3FGFwGZro@#{h&nT) zGOl<5IfwevQGw1(pER+c=balTbQP&;^BQjD5fNK(0_RU<*Q2JF^zeQQdZp@{sX*&}>$CkwinI<9-rxhwqxbG0tc^c8w57>_B_C$@42Q5_)w+|JEfPuI zz-9#8Vt2qm@=*?XBbd3~2KhlDM(#_Z`FYk7|BT*aM+%ZwXlMGKd8GMJfNm+CNPK+H zucKo~UNVT!Z z`=vP@83UxZx+6@q_+FX@NcDC5fq6}#8BxZjRJDd4xfKoz~ zGldl-85tBP)m@)!8_NjLC^%p{MhOhJ7*f8}XfWrLMlCnXdBG+PiWF#_9e#ybsls}0 zlP@zte$F#!CQq1BsF1iJ3HCe9_c`CGiyEBn{GdXKQ&%~%8S`71V{b%2s-i_#)IC=g zoYN$=p7j`H$9TIc%YBPVjF3nmDssCBmw8wKb1(hRs6Up-M8ynXW>cxj6p!o4^p3F z2=S*PL6GxElW$_c3%<1OajR<7HO%7bF@u^M)8#E!1+%#U|4gi2whUzt6*(53_4hM1 zUD3hMC++P`KQ)MO6n5{FwJ0tO)1p(_o0vC$TL7t9M|+iT8?L~!Xywcd316qZY>Hp# zmim@WF2Fevp|~^IJ75K0`G^yHxQ$a&@pT~)0?s3by%|F6{sWtE1!>{&UF!dohyN~|kn-#*>KO9o?r#Pg z0IIv%_gMav9e<+}X8fm>_)j_u%baHGnR8*Tn8cLkQ)?Ji`m}DW!gKc}&{$Gd@5SQI zdSP8YI8(-*W=!5={Zre0#1H=SW@<_>D~ZnHPIH*4pg3_;i<^&TYDgR<;PDzKCQ;giqDOjp)eR-rfxG8P85FQ8Mhd}nwE`CNsCG1CPIn|URa!Z?iR5p z2*DzRGPqJQP2^F|h}prG`Bt6#=3KJl*k;?fj}y5g9Piexk4_)6w?8@}epGMFStl4!xx&=96 z4-!C!+6)NKL~PI9rWc=VXwD+{u;XvJxaWw3p1lqCiK$of9fqE6ELzNt{YjyBc}|P_UMV1^+Pe3x{-s-}nn9o^foO6d5>96` z^&9l2t|-;Gj@jOgC_c9S`Q6qU3A=&dB8|Xl#sTONCjI1DQXg)E`+mB{3$+q32=sg+ zBJS+VFlpW*_-vbM3gvm}Ppir?(-O=skZ4QuNfG$GV96sATd-gD^Rvk|pK&^vDTE!i z)_~ifsVE=T=pjhpQ%U{5T1p8k-qtf_*FTq?O~nADtPZg9Az7#{j>~Xf?tGuFpwJZD;q#t^+{Alen38)`-kB%Q# zvWiskR5nim2yOApP8LM{9OZu=p-e4PbiJ>XKZ)<(MZrfN97{G1~bvEZ$rN#XZ_xc|IcwVftN6P zkM&RD`v`k`PkjHJr@Cz0 RIjf}PA%t`|*x$xj)T*l~rRPS;9Weoq9%(wFwQ}ut{ zBgrks)l-z2q*GsQh+B~nxRwhTTK={&8G zz5;29GJ3Moutrimf38R7nZRI#w zoz}D>24gW1atba@h}nwm%@zB0aQBGYr#M4xI6<% zClt3nfZd86v_(K&5ZE1OU<=^9r5#BtggYTMe?{ev=^Uk>K3?U9yg51sq2m2K5Fb@n z=X>B;@wD1V7qVY)US6yqj$RQns@k&S4g6icu#2u6ylk%8arOa0h!S2FCCZMA1cz)0 zEuO1vz<>8f7gCe;rjfD=vm^Ox)`^FherN76iMcr_QclMN#pDL!`RFv5x z?_V!J$yqYSMKD2=-Ex8QO?J-pGbXcm*`R!?BNE75&~YO;jB#P4x1?a(@vi*=)uM=D zS2?x#v?r{T{cZ%WDUH935Rhb6&^5TWsxw|w(L3d$A^V)8uN(eApYTd76*o)jVQJSp zdY&tm15!IJ>M)G3%BWOprOnYB058`l zmM58@SUpxzD1jS0w&>N^HPYJLL0PtbLrK6B9lQ=K0W`bJRqdmdac%J#KbSNAdU}WI zgEQen+eP3d^9oeHuUCYZiO~J)W&OVK+1vH8ZQ6Dpgo0$;nG95h?*6n^Sbl*LtK#GU zj3Tg`8rdpgqwoM2!0D|glKm`m!scnQ;up?OT+I`cqg60of?N~AUlb<) z^XipsIsG3lk}wEJX}_b$ORT`IYmPDr)G)FmAw2r!qk zPC%;IUKnBQgQ4+@&q~KK_#0ia_|*gd{cHYY#*b{R_ay!gnyRy%EReoCkz>J^(Nxct z&{#3uo|OG|g!ljjNZ!!a$+LWc2f)iNd?z2nFK= zd1d^GNaC?71-e#GJOXx=Av4`JA~WF)7tn zevy}MUmMLt!Zg`U=?<^n%#@T5kk{&*zM7ExykAfkK=We*!$OuTZ%kS>sQV9%i;T^7vC&n|Zj9BrdQl6Z);eH& zo$H6NAC=jBbu-v?f_6i~zn`xE#pQX*pw;OP3S3v9HCi|3h90o0YuNHk61c#;;5=>* zAgZ7BqRH+Ihyj5ZMbq{-@;-)5LP7h7XJ<;)1F31oVkh@SO2j=EZUV651Lk^6tU;m>C^Eq z7F>q~pf3J$u;eu)JMwxq=L(;O!^^(%m^orBL6xnvhHPvoPoq(!e7Gs|}CJ7|?$_Q^5g z3x4vZ3R#*$Eo@-Fn16G-^+XG1^SCTWH@6fYq;j`gHN7}CH4CHA0j~KX+>M9OgUMA- zLTeSgkk!|?<~FodD=YHwV&YcyoJ&xA3?*A#C4@`m{a05&kL#H(*+&Q8KRF zfdhgCGu8lXjQoA=w@0vRe_q$jYg*vvHXx0cVkZpSa}oK zYzaeGsU);lQK%x-8}nR{l1>captI<9H!8CDABn`Z!Nhpl^;bdB2@XBdVr}U{H>`7p z?r)3Vm|~`}(I_K$7le+R6NIs~&opN*n%EmO-%bgg;ZT~vLfNpesdPp|sjF=bk?7Y1 zqINKu?qzb;H5}bPMO@*3=Mmq#U{|y%3{(zbd3@2~r*U1$N5Q+DRP{bW+lpLKA?;{< za@ve)JT5>^E=HDlG?8|jI-$pH#wOBEB^Q+={i#PsM(gGSqC4L9(g>pS1%A|&f_*wj zevcJ16q0izW9m^=q0u8lP(L1V*9bWVnm>#~ILT+O9z(#}Ay_uoe_n1^otWnMOF4{Y z^-B;99ktMLS!@pF z7SKzoWbK40z9V^&(mpX)kOP{LGqBG=vdZHUHVT9gCtm9Cx2lPv=N~XfSNO7jdLDn| z_+RE6A&5n<-J*~F&B|rn$?|t$_#5cxU5*3)Bg_BD4JSK@+&(>&M9dV}8Nc+jp`EYW zu8SuIrQy77D$?SkFf;*q?3Co2>2{ddM-}mXC2}?sFYN1YE@3r3FXZmKVi>|7@49mx zw2c2vj6Y(G-o^NzGgOxiV=B_c4b8>*n5Owg1 zqH1{EY@mqkC`;Ox-Gu~WquQP;*{E$?`m&fjqrULY5$e76_a+fgnSx**H)eIn z{&foa$o?e$Zz%HRFFh#rBIq5&^-eD~wUQlVX)#4j5}Q0Ll3}eEFHET0h0}InKn5do zyJ|l_d#(Q9wOHCm&bd#PP9MU?&)o7Hg^A@USp623WuDHf4Z5kk6#S_yE`5wmjp|eP z*e54{Neqd?-kwm(xlhVMD>>#ja6FGa^9V4M+@UUEc%AaWSz<$*1j~1pV=6B^4&-1C zdWkq#2*cB+&(8y=!8kwa0{{<;!oJAXFd4;sLgTwbIgCLJp74;*+t#z^gr{^A??Hx; z5)+M%Sx@jjVJnNC7lt{QGI=q(mE+|h&(z*_7`-Py^Apt@V?)o0Z@Vxx0nCcd*6@90 zZK>3=QG+kp8187bJ)@=bePUlz`Ksh|&GDDFl%?w)|mtQh41@qJ1!A;8u zs-mBHoG^?SFm?S~QXB+A1tDF9s)^DTT_1wS^InjA&s<~h7#p1yivJG2h3^%WJs6vx!Q1bJvE4hR$)$%zFr z)u!hy6RT#`1(1ELtT|@dT)`H}NaF)@p$vaM&L}gBMb1UuJu|H4V1q81pS*PPLedax z&&S1#XFCzdE}nH#>n;g?yT^c#O7FhAs7UD#c#vlw{}m^FZOD7L(7{CXWufHz=-ohr z-r-kil&rv&d`8?fzmQ#3Lj}?FjaJkmXHw+zb5iKLoNQNA_Z_2oYI1<5fgIu$e3^!H zLG~G|z~afRnw(b(73P?l+>xsv zMN~?r2c$qUo;=xUkC<^xJ1X~0bDi4qVeSVdsO9Z2HCW>+{(G&If~QNjpZZfa9z-Fw zDBZ6r3yk5=e?gSx@wagMv^2KT##n7c)PHN&3-7|LR9H$YpOAA?$X-({guP9~S^Ci~ zo5pZmJY2gL?ZWbc(I+mwh-%$@B$8?6!=inR2CzIc}p9# zPpAUj>87!b>mzPs^dt%o902RniY$F>1%d3aFA$|@xu_nf>3KX+mfQBaNv*jWkWtu&JCalO7 z->CL;ku~KS-X4Rc6vr*0EkFKXG-3Y5eTx~c>A-y3Y{b}h)t)0His?@_S}J|2K6;Sr zJU;(*4+xMyUQLDC3l}hUoCGbzGNZ_UTsU5PH4l{$L>T437;!xaR*erwOo3&QDGa`d zJjD!eKIej>maF=n;`@SuLHi;znCau4$t5etr}dj~N;^kj2*}SS&TE`Zxw4|au=aGu z$clR~DX>u~ZrUW4NPyBbqWM1PpXVw6CdTZBkB#rK z{LzKKab-6CBf|frubR@YAX(Iop_19~95H>tHg%p`%&O`xI+qA$JfC7TJURbV>Pk;^ z4p1MakYE>4 z=O82PDND{HBJ?BP*(7H4^?43WdY5oPP%%`~fY72`7mXgu&#S!e|qYVa^yI7n<|Duvnjmbg}+et-g<&X|Sf?u=l|KGEE=k``uPu@@|z0Kg352 z#PIA7Opx$|$R>mJWo3u%`L0T+Js{fy zRYoI*jB6w(`7IqXteyeVu^^Oud5nqrt*Bo_dx+#Q{ccq(W4$qZ@hUwGR-~fyB!lGs zzn^|^wCzjeU_;YnK&(pQA`(O2N%BbaYVRKudmQF=%^}XK0+ITJ*4?0@Gw}$xz2{PYReE(GYBRsje+ZwXNG_WrJc0lF2&jX{f z^)U^-z=XIKM&#)}hgb-s7J@ix=ZE}QZ{Vzus$25TiVkYZhmB~-ziYh%XPN|L%; zHS@ZUamKV9zcyQFMEZl7*Y*eoOjV}uKQOPy!`e26em-o%Ip}&lFc-u0^=r044-Ucq zLdDIX%-G{2dr4e}1cmn1k$GTy))IHFX_Bot9s+%g1M5<&bVB&;&MkF>u4JH#@G_D- zG%eCBC?y4iv3u6cYr4 z>9A|TX0RLZ?%|U#+{$5d%9C$pv<1$aD^Lk2G;H4sv>L(dbn~(xqsRY5IDsi^$506q zai&KRyY{DP6x5sp?gvOV-KYDK>sE;5oayy_o( z;U*VRdSi2-kY(!XzD~#YEyiqRjm{rdY#ZSa6Tr1WVl$VuT558la;V%D(kx;X_SAGV zU~T8SyrLHN-hIE_nm5dWJBdy876(E!J8IuUm$GkC>e|h5;|8&A-opNfc7rc(uW)I3 z#R~bpRN*F5PhpN7Ss7BNhBjOfA=M;$znl$1kf|okk?GJ&TF{8GSn%q&l9M@Rd1-32 zG9w3!x*(ZN-z;c`i&qs+s$j>&dB%yhK0p~sj>&p{1kdPNYmuLeQ-zkz(Wpb%`Xt@y z8MInBfgM590IA>=;qT1Yng^gKlP7&pW>u#zlAtY3m%Jc?MLid{=M?( zgc0&FA-N+T^ZljVHfv2Q8Wz&+y;+zVk|$F>5pu6p=*HzzjK0nZnS7^L0361a<;n*i z1Rx#Eu6}Rk&yK%6GY_6Uz;(#FlDWl`T6VcUeKr@uY!(hz{0cmTU$7n|brmnDe+{3U z`9A)(nYS^HzioV5%PwTjtzbDH?;X)x@;2Zu=4T&BmjN38CYUo-f7K~Iro9vO@Wr+| z;!p+nGxwLElu;sEPsp$;{8$<&r-?@1>nN^KiM1Mo=}nqYdpZm$r3WDFvLe)pT0EKj zR|0buBw9acn5Ad-Ay=t2mL%tcTJUf+3#;*=UxIg^hNPX(=wf<2U?D4Xb@+RX;-ta> zLRTO??54;$ID}>D5v{mKS?T%nV!_yRvqn>wof_O{VW+-40V$xFw5$nc1DFYDN6>?R zHO<^bt-PR(@->Q_s=Y$J5w%i|-a@T;{7gENHRy8M>$1_c;!arBpqM&mY}t{XzmtbXl}-D#a={QEwb@qZD{;^6MFE9QJ+taQrDenWp%z_kT*n z|2v3F-7a<;p1aJrv1J2Evpi@;IB z;|ae4>=(chkcm&7azRCT5s^@{ZMphtE)1T7zj?mhxVdNU1G?YZbh;Y=y0a{*{*PO?E*kgiaIouP- zO#}`e)s8C`yQXv2tGpc8bTfo6xSSOc250R)s-epjIOld)JL*E}RziX!ut=SSNY;+P zYUAJ4Vv<}yEL!n$CxjGDu@{}}bNz=Z5^Ft902bGTdGQYKi6hu9vk`U7@ z1|Jg>+7#yju?0Qn>~Hd6us94X`mK$GeX#&zQ;)c&1V}5cl!UCn?Jli|)jKR-bLDCW z3By5%5pitv>{K(~fUO+Fj-aU}4fM%Q;C0-Jvf}`4&~`OXDX6;0_dSh?r130vc!@m{ zB`ugSTDzjKAQ?t7k{~nL_gg<5woe>o)eLQRDX*(PDN>({D0VIbZ<6yZOhM11NNh$W z1Im-=VP4z1y+j28&ite-2_Xo&ka><6xRAwbX{kK*QV9bNGm5l2=`s4s+{g*eG1AHI zjMFsjOjI$Et5Brfz)5GIs-=qVtx7U;#5t~Nt<;{{i z$Bf?|qeu+%#D-b>ucdIp=dYx14Z4p=C&P$QDH?Y=kW_=hGkBeT$#<&*>e&m&kq0)Z zRNe=eyn8DNpR8{gRuPBA09g4RgD{`WL?+5G&YHH(BNPoBku8_@$q`95%a3cZj+d zVm>!_A`@8*)m}*x-nU1pjGxy|y{3o-gO-TkQ}o2a@!lAgyaYU-$^cFbWl-;${2XsD}PiU63gPyO{>{kT1xv)DVny~>X* z8#QE%;I%w9566y!x<)ZJCHMH}Vzh8?)Hu6x0rBDhR5!A)6Z^nI*ZkTmc^SV`fGN+WH^5%kTb#dwUnp>qS{zPO59QU@u~Eck$ePEnA{lS_9C&eKf3rtU z?q8#voj~&Z;Yd#xFaJ;H?oWyMFEb16pyL%-&MB-uy5a!((#GD4ra$eezd;M%t0t%a z*NRnu1+wqIs#aZ~vM9#(Rc}exit;N?j;mP!bP=+6pH1&JK-6MsJ;~z&u1zMVDn>n%S&I z=|u}!OXR3J5Dmu9$R7{Ku8?T{_V=t(WBw}J8TLzE+7qQ%uMUfT^-F7K$JvF;@6Oh? zQbLXs)vJ6UVJl+vJBoG=70H&8ZEO+KkZ}f#@7I32%E+PPa2b9;iK48of6ZVG+uF1b z6d26@Z8yvhZjzJsyXs(&!4g1*-0(g7k?NY`>5#c^JU3pF8Rj3UhY*>VQ}-yka1_9t z!?~w2Q-cfSHr6^|phABi`J0{wtAa~J*zOj`1jT#NE4rwoGZE$@K8FHFaNrCvyJ860 z#^PLI=YzV+bDlM*SjZg|nN?3g4kVL&1OKvs!8PIdhG-p+`>N_>t?Q3E;2@w^nNJ!y z3gB;lRs`UU?BD_*8!-W3G{8sWrSW{5DI6FEc&dudzRB{LyQaHGIc1e`Lg@Ccp*@j@ z#Czlj9wzNI+rpn(K-Q}j!aaN~LRdqub7|}|A>LyTlQ7yexGq3Ab1qp$pwoQ~-8Ml| za2C!KU!Ok7e@LZ%5svx&uUID{$x?#hkR>fD$Rh`6>v6{=zFZo z<3h#U-tGd`!|@8?b-rY;9HV9eibt$RwMw!ESh1kxgt9KH1NaiZ1BhZn!n~on`Y}(} z#zX5#eXyz>U6eV0a&GQAm)zYNfCPjxis?As^PIOptO)P;AL&HYkLXH!S5|!Y2 z=HAW%jhMCjH6zK}Hw9e_Vw4<{A*RGZuc+<-&doIsE*h`i$ymB}Dh;R*x(*n~WpZBS zRO`dHE!Y>kBKnr6=pvY6ZLxRmn0B-|K=$rTkBN@EyJLlDjwzGMdvct>qb#YdM%%Cp zYv0D9hC2`XbFE${uHed;L9j_H7(PL`wWoev=2ED1w&U_aXNvBdh3Jyrfu~qQiEaFJ zhIcfigw5AfVdv(FS;S+03QQ?X3@SW1&4i+3GN}1K-NFE`!fx+N5QqksoutpIu$Me|>+=aaAD7bR~P ztGQUdRLyGX*pK6D6V5t|tWdMFS54nD2q4&7{GN}JcF8rJA+DZ6YAYLxYzU&sBGVEJk03W@gtrLlGa$UIflz@_J%>!4IV3Pr1i2Hp+RBDgq4nS%WeDlL7r^bFApl$*zsdzm@(I(*Z#uVsQ|=*Zsr zt@Yp<-*>ll^eU*xmb0&niNhi*Ha~~=p0B)`V--myMXE?Tp!&H;l8+apA94LqPNm(P^B;Ms1G65Q%Iu3pyqFX834b6^H_3;W z_|DhJVAJ!hmG<@inw}zb-YK^6ko3G{Jb4tQ>`g#Kfz2w1A3!lJkNmFCBtQm=nwfLL zb9}Ykvq;7V4aTZ%`9bFApkLW_49B-q$=Or)eD7b{vJ4)%FpVp?L656n>55`PF-r6U zE6HrF<7&A#YCMarn0+d8yY3Sc8G-j>5C(mq^=k=jw_|%I%qPYW^o<~!rK!gGRtS0c zU-*7xWdro1<}Tvc;T1HpzFb=-pc>PQ8MWr2{EdM!u#eJ6D9vVe+_^jK0luBX6NDIyJm#Gy}zpVGWP_frP(&fHP^Cww2kz|_Z@Rrv$MFb zT_HBc_9+cUDkhzdn9eizP&Jd6U~Mxrzp-!U8LK7%zTOAkMTG29^L*0PA)y17lBy+4 zS4fm*C*yRvM>IKNN)KH&@csnSphwaa&R~lW%>7SaO_fK!Xt_a6P@*nS;&NKvI&u-v zRsfh8OwoSpP31BcuST0J$e=RP;>|yFFx*2h0!+jnx^gz-APDKSukBEzI^@}nB+u(; z!|Cpy8>P5ziirEQDT>l7a{m1{Hpjj`X#GPX9F2H2dtWDi65+p$>~cF66)y#W{o8dD z{<0kF_gLQ7?Du~U*+ul1q7VUNZ1y{U5YxBQwMo1OrQcA6Nvbfi!!8&weTH#8N)+pg zhF5x0gDb|TytHyxsQTBsbt9Q4?y8;a)@~?-HeL@9Ol#u#Wf0dL>j7n#*#_#Q{?1Zr zg}t&sZy(#q%c+oD?T;DivJsX2{cGPd=kKrmpUAF^|5U#JoD+;PVuTnhUjo7YpppJi zLYq3hJD;MCYMr5q302iW@}cgro&7^=8HInK=snQCOv{H<`tAf1`{T=g_&-jtfYlxI zSrzNgi+u}D-&N+D+I;xgA!D^4+L#YkaJMj&ruNYLSbp?DS#(gmD+Q*+d9kPN(Q zwB%g!SX(1z@~+#(2Xy5zQdOP*DMp?1bHEq4>y1CD7WEqFR25|TIY%lL69!VzK0K)P zkO$gcZS!(_-*7yH7E%SX5B3l!&zx1^BC5@EkeD7&PiE18rV&a;j`Qi~DFv4+%i)q( zE%=T1y>sE)4bpSChI`lPE@B`&JPX%+|bo$xKJ5T4B39%4?TOvZ}$XHemrIHt0&-Pku}YXmTCPIzNk3<*qR8+uRtVx)%f} z;)#(?RuHv6(>X!2Aco@^t#OUnX<=lGO6nKkyMa(?ate)7VgWXpsYdkqQqPKuOCQ%g zkPfv3ez*Rn(l6Vr-X#ypGNDSpL@?nAvA0wgDIVao;htGqSz`>Hb0B~V&v;rvBC?${plVmDKq>BXz=_{)I(^uAm_bA~Kw?8q!)=2U4W?f(Pm>X$4I;8ZA6!L}C{g-b z@4XXV6&C1e>l~E}PzCu)DP(9BJK0Pfu^1F4e&u!dqe@Vr*X%q~;;aU^t|G#c92=Xm zR7uc^>~Ey9)W2`JjH^;=M1HE)Q?^PDTXS@JXVGl=u^Jogc z4jMWetofcX2I{-lcXeLfj-K7D1&IRcP8xW)<#~x=e~o7YZJP_rLGlu_=mMw{E^qIJ{!Ac}Ndmu*5y+P3^~gN^G*1*rv6Ifc`1n!P4al=anZem6kiFZH~WI68hbQTblA$lj(w)+|C zY_*hklj<$Xrr9T=O1%5GQAbL(5fhRZ->*)MePD^OFA-kOuguWi0YvM zRioA(@6cRJ*KFi2+VT$9tP#fpLHsB;V%sWLgIJGAd9llKhY-A~HuTXxt60H#cQ?91IaY+Os%!K#Ewpt%p|-a7KkjIsEF+)gDqQJ2(29A0&I^ z@`?IAAHPfX|B=amlbnh*PS#0cJM&gA(WShe zAJb#3zo4wS-08u$cdJ+1Kglz|!-#@}M8~UCGM^+ET{g4H&V2k_!3x;Fp}oIb6lu%+DEZ;`7zb=5kYmGp*yYlG2QlZ7u}`StKh2$`4};9L>ipUgiMl zM0cF%+`oar)h_xPw!$dA5!y41Wg&2?!xXzwGX``G&!uV4AxLeKOTX9_Vi(+IWm)N9 z1`5xJtYkC$$D7{g3+5_J&>bHr)D58rQk$j5*^Y3V$%H$#WU&8 zK^(>`e5&fgixOA*4vLGPuyq~8h%wM{VcB&MjfZ^jBp_Eq;s!P>KTk^mOC)8qvZqie zMl69O8?moQVSYbT3NI#lE@k9zD!xO$Hrj~o(oBq#O5Sthep!iYSN|(cs}k+xt47fM zz8gY;4Af1#26kUsU8nz(zT4M=rZIxD6(>xay?sK5ssD#x%^J3J;vSAPP7Ywz5$jzLy9dFCb_+;uLm=f&#bT{hjQQ{#VTB|B>puNS! zy^E5$7PWW|Pb_9T%8?wzQ)pnN>0`+qR;}1ZqYP2vE&TiF+Mt|SFc%WBX-b_up%Ki% z&Ll<3x^fS0QEZnOsc(w?Ty!W<8RWE;^y$5Cyf^o~14N`&)h(mccmmKaV7hGm>9w!!=R~4J0%2NjJrJMIAh+Wu9?pledsL^5Ix%Ld*zyXA8p_P1P zQF4}=`YMnADsCmp!hs2H`pWg{>lCXH+^wQ;dKdN+O&-V%)mE~ET%UntXI!5{WDPP&50gGG z7lGg;L$f{NP58e+Eqc(GKzHn6zmg|SuoRR}xaYF#Zr8YK&pw1wZLA!Xif*^|7PC6j z!@#v?xi&cX*R8U#(4(ZD#cy>qj3vV}kJDf};mdKSYkHEo2)>MN9Zl};U`LpzYLIh0 z`(myPb=M|vCMIUOIZjElSj*}mNI^A5T{j?C$Sv095zl_+9i&uoz@gnVR{2DN>yGLK z7e(t}e)SWcZnSlfA|M>UEM@=6=-Vk*sY|>j@M0I+8E>!XS4l!eK7Xohfav*c{jATZ zI}Hi;{PVo#$!-~^mykzd=65$lm0TnzaP<-yUxvTKAPOZQ*d{mCQs}U=UVu%`4Lu!x z7KumYpep&r|Es&1;q7GUK+26mQ9=*A_Jks1N0jV5pO~+9Cbav^iBGd_galj##$1Xv zpWBNlmacuW@UhtRsZ5t4*qLRmz;!mvsAtBnvCwHjn3Y~=i5(oojaS&{q+H4^<_Dc6 z$*20o2dSBP!OrhLNE)o{LGv}ula-lhKLMJ0Idb-c=pxWb#t`1|6a+ll38=%pet0O# zC8l)lJrn;?>;F5*<{vZW07DreMOy@XH0Ioq>L0v^@<$Z@#8k0&WugjsP#H_U3W1$jd*CL*KR7dEul#kC9P=L;?! z=8t!?ZsT+Ryy^94NB)tW^WKO3efsiWzyF8f`Og#clqhCIr;fa!0U+E_&Y4!2F#M_# zeuif7wn7q1PF(Gwi(b>S5_i5B#oby~l8XDkvOaba;OFh%r}gbG^ZG&T-;F$f5ZW)5s`%R=x~|j!s02xrowi%Mbm?5r0|z7 zBm`B&#CSSrymg6M15wLvPr??*_^-t3^pJ*m&q*+<4{d!F)5j$}!kmsI3L$amlPi@Y?(j#xRuTKO|c zmkUA+k&v9eRwnU*`-Rl%rhpW0MEgy$CQ0X|G?B6WfXg93hAGy-1WD(VR9>@i_UBlp z0+P6@R_nmB&0-katK0mu}9!y$xT@D~lp19lwXQ zgE2lohM)=?=l9qmD^nmJ_C%fr{Zv*2L_{AMNy#hnC&sEG*?jnCO=T` ztJSZ0*XhL;*UUnav=u7)?)@@o%)^OO{5u%L(`(@u6Zjlmn{^DZRR{jWrrCrLcEjk6 zJ(7Tys8IEXloq~Xal#(;exLLZ!RNIjgnxdgtJm`||3u}+$T%m1)0eQ4su0dq0kwsr zD@2Z8rNalgl3ywCdSSD^Y8oRlGP~d!2OB-(WM?vMM`yb1Vc|!@rE2wtw5x)9d)tN^ zyuf}gPx;H-l?dm%JsEkxkj;b^)x+p$6U{R_CCBpX`Y%-6Y9{5`H|Q!A+EPZQvF-vr zMM6q!!ieA%%)w3 zh8)#WH{@^Ky&{T`n+ntbPAsvZViy`hNtth@s`O#FPX8Zm@8BKQ)~cG$V(YX>;VuNw)2U5nzd4mK8v%SC51;w(AD@Rl}kOBkZ|iOiN=i*Y$7dUh)P zlBZumN_Df`Jq78)?c?g$|FU0_Mr+@q?(nR`99_;0!8(}Ovpa%;z60yOaR6gOEGyz` z=(537W-l8n1b$EH7gp;XY27dg90-`5y36;O zZ&e=|A&BHd6RfmDi{j2=6+m%p7jo*QM=~vP5lO$}_AIOqsr9)4#Ll|6=N3_RDlo|e z6b-_m7^34kqtkp?SSnG!)S0hTUK2*!lZOQe*e0bbbBF6lv(-!C3{u!&h+Rv^`;l;V zCqgH;YY92EEN~e}I-yB!PThs&7x~?D0RIsS1{LG}XqCn{IpSkIh}jFRo2-jZ0mMgX z5f`4^W5Ul5idQ=n%~+fmG6ts|lflZk*J)|d_oOsVPGR56rN^le9n2Lsv!#~-Ng(g^ zIJe1IQ1!_L8x6e{JKLj1uf;!-+xrzLw z#4v77BYw38s)H0$SAeSvLwrQq-x4X?8}}T zdO8v#C|bLuAx`>Yvkn%~0GqI$x;xZpX+xLIr_x&(&J5`JZ9d8wtg6w`F6tMtuosM;G8<(GL33MBz=jA{ zg8g+-{4G@UgnLxLhpS)f8%~Wv%Ew7N1`h#i%pc7H*hOn1*#~k#R5QeTx-e^z4vF_L zSppk_my14ydyCZaYV&?O%+giYP&hKwrz!n+rvM5iYG=K*(vIDe_*R@AuQ38*&mlxV ze1R(%cqZBAbFfpT9-r{#ob-3Uy0wBlot@B9h@^G-#P=ktp4D!|lXmofaKn2>JY?o! zR{}XQW!@%?quR)daMg^9R7SgvDP&9M^0^zB?ql)a^Guw7j&p7PVC>jX+81x2WjJ1- z!5JOl*Im83U(C_1E-11E>RFyKEV^HCkK2)k3+0{cE>iT;p#pnCP;YH9Yt z-{kZzDo{N*+(;c>c06Z#C6amx#MnPw(%RcmO3Op-fDas({TawvGL|qfK=F!0cZSBR zYSa_0yA*sTLV8TxF)q`3r6eh?+-P{PK&`Ta0ZQZ6_=-Mhi?ynBpxAuh_r@r?tT>^M z7SE~se3}(^?NnD@hJ`wlPPSWi(GFPiDM^h3^8p!hrN-JURPfV%Szvp}36R{&ZHfB+ zDc;I#_)Wpl`of%i|2fKLS1U|1&D=pR$(JAB}omS`|QlABY*e#s2F&{g(y#gB12|E$9Ay z@AS)BmUG-uW4oG0a?#-4T}Inz^)FdE zwljqAC*7~=8Uu&IBoK?WRu4+|^QP8~A9tC;gEl447Yg-M+(Jp6OOY{L2KEU&9h8{b z)Qt+U9vjgrDHRd!jho`M&*X_s^FEdJrgqEjV1195*Dxg0w-Mrjf7ZzgsQvv4A-tf! z)w@0WeeLg+|LvpxFTOJW{SINwe~)JIToDLPq$kxz77hX{ zYL6|86v^e-3dX08EfWJxhUqJin((56Fbj`^n{xQ>Z^fvO16o#W zzIR8WH}rwaDpItn=nEJl_t;sdT?8x~g%zwwo)J0C5y0TjMPCGR|C zq7nJ5)|@gDr*cKB91S-E4}i9e04h^5I#@2oed7*qCj%#rHyWWWVzaAz$UwJ6U3+sE zkqFLUrg<>`W`~@8vQ*?(5fW*^(;6_r0iYLx!$Xt|>lti3yvo}4abG%ml2yd#T3aME zrz;Y&E5h3s!J3O&VzFV4J}3Y^8J32|qt`q8zy-2)5rmHT1)^$PAox03Z?1ey3_zQ@ z(CNkwKAv_;6Q0sn%5VLX-#PaSv%^6T4i?5xb^_K(?>C1QBffShH(IM_;r2Qp z6EvVB^D_9z126~Nu~;qB>q^zi5I7$1Q9`18+&~V^1+&F5Ak~xnQqnjS+9?lN2Df5t zg5?4xPB1sPvS#81!+b!=gr^YJc6HWn{4ka;R=h!}t?Gi){@ODF4p`xpnOU3r>K!M# zmfpg+!g2_@L4M-e>XtC05A@N3W(C@QpKmSGH&BmCY6(F@DrphUYNLrPLCVNb>hnh67H;Km#UQajL?#qfavZnG8+jN)Xq~{s zV2bJisXsBVvCz}}WUv8}$W=^`zW zUUt!g*aAyz&GLX-hi*~U-$QNP!%sk)VezBDgW;G+sUszvr;SBaXY?Sih4FxP8Ow}e z9S^+Yt=qn{rs9Zf3oHwW*+nZQJb?6!V8nBD*x=G@J9oxPSl(b*4Q`&yK4TT8Fs8Zh z#4R}!x=ps~?_S?ZGU7)+*AvAHm!jrm0Ov%~_i8DQj0O1hUVYs8uh zfgih_uk8vJC(1@5CZyECsx*rnR8lWDB9cpq-)5M`lQ%fGM*0O;*7Fn6`5Bj$;k3@p z+(|$N+!y8Eg2N45q8X7*2rJp2f^DskpjEsvA6=86(VT)Q&SXhC_dy+y^LRq6vJOxN zWhXvJNfu+bIJaSQh?S<ny61q(*>%`#Es;=V^QLVQK$7CHw}4K~3V3kP{=e zgRnnk^$$q$fn&3 z^rHeJA6XPwt4Uw<+mSdklActoJm$kaCyCn{%PVCp3&?2U*jsQPqcz&Ws*b5)A^j>R zT_F5U)1&DQx8n+Zb_YRl7>~dSASmttE9YwYng+y@-)))dn0Wt7ezHg=*h#Nf+J)#jSLuc*{%(*qQWPXMSwJ)(1TCRNcU zEvwzT474X_al|U2>zvgyEos~}QXG4;Upp@0U&P6upEcV)cf{M;&oJRAY)grpl?0A^ z>Y5+AGvdk;=+^ZLLBROI^nP}l^yzqrP$9CT&77MT5{;Gx3E|h1H|KJNNTYk%&W!iS z{7(CFllP4}JyHP<*jF3rYO1U*)%QL@u~9XukFi(1F?VtX-luOA&CJ0T)2H@<@NAxr z8t2sM+(Eq^sOIjdoGU}hQ(|~FhcgG|+I_a_Hg5us+EvR^L4^Y-`&3FN)s%&jPB|kR zYu=CL`CbH<&xQ9O$(1gQyzBSdgJ?SW-S0~OYm@Rn3}2-78iQxS#rXqQm@9R|_+HB| z*YO8@@!c>v{<~K4@6qs24W<_)zsL@K{Fw|lE0YVhbVi`+k&EH%5x;rBR~>FH z6iaPZH_4`KuoGH7IEdY3B^eany#5@qIf$KMQG!)u&Lsq^AYOocBs_Ip>F?3?)NprX z{~3&jHdch~tFTruz8<^`aZ1-U%GtMSG%>4EMhWs*-#!t6sJB#t&?tddJLTAO)9tmj zO23}@-Gf0LHT|EtaO~MRXuX5A$_VwvBxy=3gVMKj54?_dYQe|La8Z9rF`?71`&bTyMd8a!rqn0$` zZYWD)GbUCuX8iua%k(opx8jj9FZN7WR zi;w4s;RilORP0`Bce9+=6KVt0>fDG|+QbGoSCLm{dGy8jjNtZ}|Dl?b3{|h%+R2_K zO9=2@&emOU@A~`+J8wp^NNLix;o=ghuX{|ZeFA7$X`>z*<1(7R7t-JiR1kzRLC+XA zV6WA_oNK_$EAIMa9wYbK5_AMRxQ5Uw+OPyO*Cy(ySL-~<;F=nsvy+wzv)M-kotBv9 zO(CVd*VE`nsI=CpUM?lsalm-kSQNer3CCo$J(`Qxpm^c>s_!m{^BP<3bCd=M$olZb z0Ywan^N!msJM)amLy+WksE^X|smif~@?MZ9dctjUVE%Ybrr*UrR@|{F(PSke*pc{9 zXaww^fB(o5<*vwL-MzFge`7;>-LZe9o5OEqqU*OtgZ%Ac{v|^YWF7Wj!{NFPX}@6g zzq;Xnn0>c1_ab5aBlYL{?_X=tqb~S!Vf_dD?myBEkuisV5PuosNWEsMiQ7=LoToe- znOP{`AGetx5ijHI_bR#_@wJgIg)#8~e#B;ZuKE2ug@lNF>%Hsz>W06w?*931SPpUJ z2P!97$x$jwJ~&TK{s)qK4;7Q_y{=!~@DH@zhrhbv*8krP8U7o|O+`d=v}~1DRUMZw zB2Xn4d|N6SUD!V;($M?!NMG+NdeTa9ubi-LP<5gyXFu(`-Hqvtz-7JPq_RE@SWB7q z#PuW1^~=dqQ?n0ZTZky(H#mt2U_PU0-*zuSSaF#Ar-}Um9pPkOk5H0$FyA-M)C_#& zSG%G>54HO{bTJ+kDbEVr6UGw{8?sK2kdvpqs4Jy$yaP#fJGxmNsCr z8Y`H&)L|)%lMy0!jpZvP%JByG>6gU3F1-}7cFX|jGXZsZkob6ospeBpKs}XEO^M#C0O=9-pEt=x zuD^MiU*6>bRfSm_r96|;)*4+PmC5*+5Z;d@Uq@UCB&#BHz(v~jtxb$L!bc%q3Moyu zi^(b1Gi?E6JL7Dy_<*jn)-z`>LQZQc>m^hp|DzYVthstF(5#oC^J|*oCsZ!?JjM?H zBbB!guvy4~t%+)Uy-4Kx^E5K6P`o!SmC4;S1)ofZ*DYTc;GjgUivZFPWQC=}d{}PQ zg`n(|@n6GSNi-E#dSr%g8J!zT6G2XUBg;X5SgpfcCWS`QRTdRsrI8%<^89P-X~4w# z703S72LHprZ)Owd9Cojfu zerx*QK>?-S8`Q7H`8#zi^RI{e=ZM|qx%CNn0Tn|WrrNP;^d8JLK0oH{uaKaD1A;)B z9T?k`e-9K0E4S3X*ZHfF{;cy~2@>zgaMAaaS@;9TG&q(g=nmRK{NSz75-W2beuK;{ z!*~2#EX?nWqkPkuCehy~yR-Spql9wpBb!b4 z4(rWvnC(o+QS#5nAJR124=1DxsZYwA!oF?B*K9gb$;oX2}MR>E8H_8TEiaR>V(O(NlJKN06F zt(Zv|$N0$;!>o^I$xJab+elD9lE+U`7ebJNfA>O3>DshSk0JaUms$7->U!VHUW{G| zCM{o$WG}5XOj47on6JIsvwT|qSB&lKz%c36?Q2%<+Iz~FNalbF5qh`2Qb*SIdu~yB zL$vT2jscJo+sy*Js-pA~o6Hdfysvpz>!pg)8_a3oSn)!aVmG*+K`@IAaR)4gNDYkU zXyPAg%!PcYsU%2GYvLf?dU)B8FX&@j0M2R4nh7r|u?Uo%c>G+kJwF6AmKwLi&4V_Y z6Ro;>5|ccQbFeANttq6{&wfnMJ>w<@_4o{2L8jqr485%H{%B(N6?3;R9~k%N$YXFxL+5i(Q*(1ufgMlRqUV_8M5QCV20oE;J3BaiHjyu0zcC-${qx&6k7bY9CLfPL=atE>OdQsh! z-iZ>xCht+ZhM^gK*%}>w%Fi3Bw9%}!P2R;^O=Ujop z0h|j869giJQoH*lCl;@-6*P0l_|j*=#l|?u=TRF!=OBOXHAyWpp&DZ}`U2|!QD<3< z>Th#jA45iH+$2#u6%qxX`kKc&qdv&nL&E6jPbAK{ct$A8Bu zM5f{_nKRb5^23Nm>nPOOncz=)EF$OwVr1{cS|(OAP>4C#m6Z%uD1n~qDe2bcqS$mb zt7>3Eg+LQ0^>e5+L^*OWo*HZQPO?j+F79Na|2YWM-vqn7hw81t%mme*^>$Jx-p$v$yS2`|D`V;b3uU;4 zcTF6QQ78HO#@Aoot5-=L3dF-$yo=fhJO{?X6C589L{iWEyL!kb7#ilNoh)Q$rjhm@PlHCL!B@ONjF6VmrCKuh37_NF^xK6@zD7WPx)FiF!ZwV`=qkf zzd$^_-nHYdht*U)z3dwm$tx;8#Hj6SjEYUtsN`~$>l0R^4BQ}G(Z}qb+6Eo)y;aB> zLD4)moZ6-D^1kgh6Mv&Z%rW~`_POjW0jrEP!aA(J8gv2^zC~Kws%v0N@Ks4H3davhnIL#FjN~c>t!1$maRN#V3?^WG5qz{EydaqVbB`!y7 zT^g?hWhgs($6ZUx0-%;EwE{dfP-f0$lrE({qI~(8*=_3xK9?&nt`CvmC^d|62D6Ew zwf?b2c{NXumc+LPsRz$3KhL~cDG6X4U1>8DM1WSttkVED7y2-8rFjCz2jy6FNR}Jq zsfL!vpKR>y145x#pD(T}GJLoqo-P}Br^)V5d} zkT33(vxmsNy<&TAw8R7{8R4hGraG$|&XHFkU|1GY`2~`GO3aV$PSd6ZLUDp7y%Wl?d&NF5`alFkXwB?I#ap zP)eA25O=Uh&+Cj62Fd07d0CImMl(Xp5?ELp(45InJk{#D2?w88jqOD!FsZ#yIO|CC zj~Ig!lD=H#Gt=7KHVc(tet!`(uM$UxyeEG4!*bZDad>Msj?(P&GSsz`0$n+7>$_9 zW+|B|bUIawo=qGGuwFxAZ-?EAZfIQ5+z5gFdUzQVA{F_D6s7&0aOdS;+SAFoD+lj&x5{W=6c2WvL{cxic zKY9XwfM>F@6bEK0dF9gV?WwZZ<~S5<>#{lhypo+0X19&fd0sueWCnvpc(7`DAF}KX zy?%Y?y$4jT$$!b2`}QFbs$0C6$F&pAdy4{~0#~9zh5Z%;)jQ(6dGJ;T&;Fayb#RE| zJjS#^OKew9hTssy>wLeLBIY7aS<8^%n}O(R`OQ{h!fElgNJ=NF8M7wMVgmukwT~A; z2ryUBVd_#^-nrB#<)CFeuR5%@ZG;l?!t)m4R=#iDUpB33qK&x2%ROqi9?4)0^OB=7 zFHbP5Z)CBn>)^#O_dfhs z;de5>IX-r#1|`l1Oj1mV#Rj;t(BGkXGlueK%b z{*;MzzmMM3~2#Dcuh0HsQ z%%o?M*k{HlIi_PDiQ0SMv>=$w6EH9*qix9>jZe^!$n7jt`AKA4r{2ajIKZ(-mO`h% z-8vCUkQ=V8c`33stTO!J(wY(wFB_;g@-iTjoJv5u%r8Srd0$XFRhIIaP3dZ9Z*Jy< z-27$T+Y}lgxi<;oNl4~YyjJg9O#!Wvf*3-U^v<}Y@>Q+FApm9af*4auzg|88Nu3)x z^gJX$S~?Vs?c}PUZilvQr_buFgWHK1oY{T6wHdoOqNeMdX0^JUPG?7|H9jLf!t}Qq z>h>VKgUrlNh;9C^1C(P+;YVyP`^kuoIXq>fYxPm*wrNDO%A)|*D?t@fyUd@*d?k-B zLWG#OA1?2Jmi61}`U-GqTn11JWNzRr8E=DxZT~m}KFuSop5A^$a#F>}R^EY{Jh|%j zYM8QkSB%1@fGuZ!;&b1$1Eug)=58$I_pCkN7y&(e+J;}MUS1COZ1KJ#L!IC`X1<6 z`6S&jl0aO|a@gXfmD0<%w6nkOB1~&7hF0z-95_nf+#z7epN>nM%RC_zj>Kkk#kHjc zGy^L2wfEU8;I41-a3uBKQEpFZ%8viz54l{IIMs7H`po25H1e;oc|~fbGiaR^1&##_ z6{EWP2gcDyU1aA`s-~k6h+at|#Y&4!xFlS%!ZPs$wCmVShv)hBp@)*!tQO=zBIXVh z2HB`RhY}gyy@{W#Us9=B^uSrdn#R}zMByRPs#F;;eBqq>Qg&6Bf-F|9aX8d4yH$;L zwkev|v|9oWlLMi+Yr`MVok+kZEREz1cdm@YOap`4khlW5Uj^M=Uh?)YUpH36HTe25 zi1#KkM9zB-X&%>RSA8Cb1$8AA;agEcm;9Iaye(%Y5=x_&GL?~23*c5jU|d#j-kj^e z3aKEPN62+DrFyVKI2d@6H?^;R9DvGncVJavBAX*5Nte0XL^cznC(acLngdZ-3wgnv1nzXK*`kbWFY_ ze!WHLf7@FTYa@X54wu4Aj~n~ErI;7F+-9RylrwZd_w*syZ|^Do*hvXT6t_@AE2oS9 zQPMt?*HPxMP8+)ja*aGh4tfTi_I#NiiIUFK+-v_l#hq~NYrAob(#^O zeDW@l^4nX&mCBrLd1IcCL!c$bj>!nN0%Mg-->q|Lxw0uGM5Y$t>Fl@9DE1W%a}aeB z&K#czB;PWE5F+KIO6o3mt+bb(urnx}q(0jjXAg5|%s8sKx=p-( z4Z4Ag6x?y)#|EjKKb+%SoT7I<4?&T%yQAbrzosg=UH`gAA_T4a z%Cc_AE)XS>`}#gZhu2;91m9OGK<8tJ$_(>HV1cd2OSdGQMqX8(&sNv~BgCTsYvEqP zghyeuH#O%AT|M;I_n*ue1luOJ5U#7ql?S3y{9?K=2Kl%$XD(vkykWmMs8$0d z(wgGvO%h3x3?RJ-%9ePs!cKAQJgD`@ui$n*eHW;T8JIABb1z(*CQEr+$27K(8c79( z*HGmZfe1$fUJLaN9VMZ}D#T$1e=5N~sgNpf3rzNm$(PINxB~86WSH)(NC#Yb_T5!{ zoYFFC4QDA(g!+hd+qNL{ahQCFVQjYjyXjML>E!|*SVPf9aIJ4+_`X9NGPB>VLYQ$3 z{Z{3Xd3mm0dKg`iGkU;t<2|W?_j5xZ{h}<`1`GlHmn0T5X}d4f8aAWLM;dIP{P?0i z{<>svYIR8+HBGQKRDwb#B?(eH(0!`wZ2o{gUtZ7bZ7ChB;A-`*&QYnFCWCAEVxiV* zsY34|9`KJS?U}spe`5oDM8tQ0SNvbT@;}QpD@T`PQ58JlIr{y~#$>{r@V%B_?&A-x z*?+`W4pp3{so>>queOjo3F4h?0^t4q ze7xG`7V!NA{qmK+v(x^0UrC*mcH|>EvVir)+a9b~8G-|73``x)Mq!-ScZ|1%H)ufH z#%=|_rr|q`xY)u z2uAg*kX21?h?$Sb9SKQi74S!wszh)gBZsBX#ADDU)0GQ_TZNj?ZzWn;jCGv1dOWpy zJhy6`eQ81LJ5K8XljB31=w_Ja@@Bc zb$qbkG*bZVHISM#G<{uWh7DG%>3eoGyVq){eeE-Gnzx94^P!OTTqZNdaniwXdN9~K ztRp#3&eDmEqr?lTh=KJncgr3{{ZA?zj4#0_j>xrOLUkf(K8Z-2#MIp`X z=hOQD+k?F(t_nvWtZqz}({}v1KET_-oDvX*leUao9Fxy|+F}9q#{sDTc_pXBzL#Vx zMD*w&BXa_=G3*3!CI8G_8(Ko|m{a=4y@)-t7`$_gQ-JB9$j z!zUgmNs4|4U)V^qX1uIz@mTvQ&i?3#1LGrz3x)gvV~9H7@F zB6L#9vOTTsu~>ikBO=-xRm;p2$js|*spyJGSak@a5vz*Dth3>nwp|TP%lydd<+_xJ z&|&x^_A-5}^F*EHaGc&AaclYlG3W7$@m;gm61u<<1o$*lW2r^u`Pk+T1c5T5rAq_V zi*Iu34~$tDM@YGV5=KFlqgHhH?qG`ZbV~v&e&*eaUlptp zm6?!Yftu+drtp#}*R1)p!NTT>amW4d?+L4;>uqxoOb?>_w&Os(1aL}_UIHLbu1Yc> z`3-m*?3Pp|$Wgymr%t&#sy$~N?A2O8o6Ng*3Dl$1A7vRTq7r zPqUfYPgQr>aC!lg=-=0fz02^sxpg42r=s#aI%ldF=~{`L_wQzoyv61HK?mF!Ab{fo?ug5Y)2&kc%m6@j5t&@R|XRyw56**|E1#Qm`!Q#VZONA4=_5Z z;Dsrzv9X5vE_bSdj8iEEu7xQNwt-cCWkGfbBqj90i@w4x7h7S=xuazyi@WRTIn@{# zDKHN-gqnf(EhdLnIy}gqN!~t7=8m?yl91BC3VqH%iE2W%@%+|LL)B-Sb@${*T%3&+ zEJ3mGySlZR&>8(y+v;klk7S*DTssGd5kGT>KgF<)UCdQ#>JV14or(ct`0$5#a!g~! z5x)Pr-PJI`nmHyTJ7^79ef**zfn~&BZCB|BzQbh*n>9bZjg}7U=q!Ds4eDz^UBML2 zK$~HxxpuhZ1&V{&dO7{*A0O{)kZrQS`k8J7>nRKzMrrwoa|JXUSX$5cS+ek!m962B z0rC_Ubxk2S0^luAG4(6G>2U}LW6zK)bYBu*4Q6p!W}(e}?{ki@RJ*Ke+xifg;}f&H z>f>rx8z>BhO`Dh@qvl9f;qd2g3Q1Ge=vlkTQ}qGiS$t>>CYNT6v}nqzReaN@T?044HD zGeQie@&ih*twt_k`@PkHngBk|%l_78LLvCxf^CWPr)BbbBdB zbQw@UeKg$@gu_V=_vYPf)ot2jC8jQAfW!;kOo;#`!2jOZZx}Uhs85Z*2YZ0JNdIzf z-w;^wm()MLE>U%?MtgYT{Cw=;8N*8SH5eO$;g7X z5KRW^q>kW)Zk^8H=q6pETmpx)i9DF&2CakGfo5)k01x_)sH1f;wXq=RLldp&JNELL zDp1_HlwB+tl+UUskj7W_IOCYp_4Cda@ajOS<Y*BV9YyNa zwqC8_VwDB2N5S!WP1i6PIv;$+@93NBkxJJcT`P`p^l%W$IG7)NiF1x(B#1ya2t#EU zzZTA%u-?pb38+_NGJLT>ffy!{BJz7vvj_k)MXZ<%ZfW9u}eU z%9^6s9<0Yv0fVc9>{ORe=bW{)7x-j`pcz`=`Y{M&L68)g>J4|rho4EmnC~NG6 zyDj|F(=2z{#XU4A3eNAAFW!s)^}Y6Y%lEVPe}z6q{P?HFSCE`*v%&QW$tj3g%0*B-cIrkr>zP>EksREo=4VM%UQuS z+LNIWYuBflPd=j?#C;I^P~>j7VMT^yREY6{;{2F7o^NOtup49ObL90eHxju}-c1wlB1%;vnRihl~NKN1S@5@Tra zPH1&d^jSaI10{_$y-xo1gf}b;M*ndLmmq7tNg@Y#ZQtE0_OXGglZ^sZVtIXAFLA-5 zmd9tyGMxgLIRFz*&Ak^v5<848RugFQZKqgA6Rz?8%^uS=^wt(rmK?)m&-BJ;x)!Rd zbN|S}O(i$KcwEBO+7;MXjL8;(X!PwQBJsR^kF0p5jz_6l2*1<(ttG7K#-~(C(kf-x zr69MR1K@A^VrisH#d0C#c?1clggKE5J9%Vci9W587+K2$kpr$-%V_L5-3^@ag7@N( zK}nDt@l5ol7sb$3u*p9_qg6@d@YX5`ItmMzIxsK*npnu@xU`*9%VwKOk&gkFQjq-Tub$&c_4MUR)Qqel(A9cb}OkLXNlqmq@l z;0zqpTIGNMC%OcxiiX$vyq%ZUwmTyIsCgE2$OC}{C8q%qVvO9v#ok%v4*U+LBG)nC z;)r3mN^GdGu|UT);9b!(karz7I1O%TZtEt8^9Ee8Z>>TG4wQ6b81dW*j7%&Bt*~-F zp|dSYpsrkwp+g1v12piP^3yfUw?)++nFAP)vbyvnoD5B9kZqLLH z&A?%>c{I*dN`#LHHg2$t(~kkRj+P8|Zv%l3cF&qCvXKt<=WsdcY$i(h6w;;TP=>4< zCfKlRn(Nn-7xL@N8%}!>pC9;DMZ8^Bs5(V1q;dgpF63c+(AvEk!FpeS>2-a)-yV1c zzp0%E_xaIf`zWpsh6^s#btelpi9D>$&M>mXPKO3bS@tI|mp?!hyF3ke@q?9H>t#1O z!QirRS)-s7N0NjTBN1L_U%0Vr;*VQH&^`-k%ZYVL^VUOtiy*j@1Qnpvn6b7Iw%o(Z z*p9{Nos~WcoOl!2$^KCQya;VXGYQpfWv6;Tr8oe3PSG(u8c#}Vvof8$v4g%|Vm^dr z4OTjP%k!3hKq6Q(g3S{lf@ktF@Q}SPJs+(&G~uA;nfs#wqMYgDw)k$P%>%`FftZ<+ zP0fBj6hcFU&_&2vZ%^0q-Tp=wy~M(Vw=?}}R+EKY;$psSAkuB7^rw{=wNZW>2aSe- z8l%LrVvCgqJFu1K5n(Jp1dNpy+;dKUFC(LaQC$|^bQEb_-!v`D+j^yAQXso~qc%BI zqTQZ7p=o%T@i67{k-ba0PAMx1=X_Ll4^((&1_H|p`%7=3j1Mi;s2#f4KBN<1mP3z| zLFrhCCwOi|_09P#jD(^5w`#)lE6QN0b8PBu!+50I?kY=jtloM@OKVh$@LE=1Y^lCk z1AJrctUfOWm4A)+cGkw=h=1caC2CFA-^~|&YL2A+yn<5V*-6Hj!3WU6dix3 zZ76uQcDz6Buh#iHs4VlZhy4f9;Z5syVipCYY>Z-yvKI2Wo0iAv1)lDC!QE4sT_*P1 zHh&P^>*>?~K!-E@EcAa4h)cf5%_5Firk$`n)fM$OW5uVX1y%S*$f25co@(lRga}De z6N$U)dwB^I&<4l{1kDx)@uG4Hm95RxAE6Y`5u1h^LtD_l>Gn24?^K_3RF1XmIo&PW z+;JalloZ01!F|mI6#$_Rst3w;t7MA$@(Cn>fri0Id)ktiYrjJx%TEx*?Ia4-2%cAS zB~c7N5;i3gjPm;%a6?Z_i~A2=sB~%U0?@L1xTo;ccmgAz!R_ZKn^*xPsio*i9ZG(^ znD7n*r#1iH>(#!gA8NeR_EzJ z1~i~Y*?$Pj%_e*kKm?QekaL&pQxEWU>K3-`(NvytYOQ}PJ+GgQ2pJUXO%;Y=TB=_2 zj(*VMEuUkV-*=b02L^RD){C`U0mD)lCfUGQay|*{bHjdKdR7}0a1|R~&LgJQ@b?S$ zRH8y)$$}nTipcLUHbUJL+AFeX6iq@Ohi2~gI0Qg&hGq=X+6P`9;|SyyOj!N<(_*2{ z>VU`TEgYR%bdZYUoS6lJEIzXwuMf8KThH#^NQT=Kb=y0vGJoMLw}*aRaVn(PKXtKoFBqF}LzDQr6;0{bmT-@QUt^bP&g< zHL^OVh$na>xL$YDQJn&7IOSVCCMAGy#9s8Z5WMW%D7)Q0YdF94>I>c71Ag;3FyHqF zoO`)$%BPz~tl^s?$^AOnb7zAaJEAQCCFUTM3Nxe60cHt)lgv1C`CSJ3u+kgfiO%d)Mp15L~fp z=!@5q_a2$I*(%zAr|C|!Imw``hX>^@OO%0L7VDTntuRMTRj|;y=8)UluTnmXz0jgD zS(NYEmEl)iALuLXVWAxzi-mc#`FqE8B(xacWsA(cif%RY_gikZXDz|joYQ6>J_;vc zHEnIz=2UvkKgPR*i;HH)e7+D4lkU<9VlX>P)hk?VGUwHV?JHVi#iy@b8e~(b-x3PY zBV!kkw7Bf6>)M&WF%FK5I2@@}@CzbXXP6Wh@3Z*ip=?$SZw9Y6RH1M>f%;jzye=)W zF*w`L8O6JjxR4x+VZ1sSK(lq_iq0BwIyuh;crlJZHb_0bUJypcS^<3$e_(y(v&><} zra}Kk%#&%m392!fwk^vT)aO+&9~DH}h;(5h4U6mitWiT(>lkuB-r31^E~4)JN7;g2 zK>tCk_OA##x31jdV~VY$5~l!14=!Kq*BFnHFvXMgNOfBuTiAOpbbPzhhPfT=kQPK9 zsT-Z*JMaDQ$~R}(&;Vug`iStqjLFVWNf-lGcdKa5xN(2!Uj^F|N%&-?*1_YLiCpC8 zw6#UwB!yl+VgPOKDeGsOM9&!J9PHG0OEE&AfZ|apXMg~42TH=(IREnq54~&_L^IOg z2mK9_|HDorD>#h>G;$ll+wy}(P>fKphe=azs5;R+THmn|Q0RWo`Alr~yna-O<})Qw zbC4zPs!1elT6}mu;J(7gr>xss1)yGlI*w`hZG>Rp!3zpYUG@Z3AH=M~m#Ll`PHJm5 z&eQePb|7ik8-Y1S(9RNxsm=B8oQsSZt;pQxC^w@`790jdfNw})NmRJAQLqHD*Rozn z!Gef1qle@+)srHaB%>v|MtEE^qd3-Rn-e+yJ_-{Z|Ek)b2d2Q^wd0p+|Ig5Rx7rNz({}H*{8H^ds8R2#UG$fr{XY{o{}vLbZZkmyIMIp4O5P1oK}QMV zxHlO~RhWOA@+Na;_bK+-v~H+uD8h&7g9Hq#`59IFEk(94SYG=6u)jq6cMjCMX#azk zy6F3O>Z$3}U|t~B6i?SaUY&!U)csBjso>+NbjTx9}jeJo2-#$yBXL-kb2=4q|H=yM!pT09em~+oKT2f z8xKtheHw-0Pa1~uuNGN)g?w^%S(CJAjtR&KNvM@`?PVQ1 zZHkL$oPcmvleUe9B8ZCPSlsYY`?fOO-V3ta%Grf)3MkU*m1GI}2wy6W>63GOdtVqU z@Q{yFkJKsRtbIVc|B4)+Gh_~i0}f(aUKKqQEQ|pcLuU2FcY5E9x!p_w z+-7XXvw#m&**QfSh0tuat8Fw4AH%7f&A}tjcvr;4RqU`aZS6;l5%8Vv^-o(PpCNJ3 zgd6i;*E5_2+|GG(w!m} zAfa?4($b9}5+W%G0s?|cqkz)V-T6=|sYrJ!h%~-)l^Gnk!(Z=jo^^a&%k^FReD1lk z&z^ht*}L4t_<%^3q2anqit{Z7?%it6uO{ot3fsyWuIO`TFH;q_-9@MN{l%?`9CI~c zKtQ#9%s67l2Q&Gp|CLR>Tb`8y0WpfSiLJMZ4E8iog07*GwCZkeCU{CeuD&0Qajy+o z`nL7$Fwq<6sd?2D)vh#XZa-eKg;;0_P-6-1bu@BlM((1{lEuyn)QqbuqPI}pOC9(^ zwI%BDJBc-wu&Kjx-D*POaglD&Q>$q*r4GgA%wP=k@-n-n5h-RNTUX`<${;1$o;xaU z6*wklP|F%02%*dflc#kj%=wVn_$Du+cxd6WPR#0r7&ahTZe`4~#0uuyMzEbvey$am zo59SRe5ZR)hcjnlL-vtjVk)w^2gksp+-Q^(yl5?DdSvH3YUTdenUujPTI@zSn#vzJ z+=N)%>rDBvR$Qn<4TGK%F zf(!567X^ZT)2uyU&=AG4?s@e^b0D@SEZ~1(lg!2 z`{t8Bxs61bL`J)5XAyfd(qLYJiODC-u6I~J?Wf0Cx#G_e z6pGA#?2EF^%!&-ESdt@OJWy@<+Invpk?Yg45md6|)fO2XMCVKs$z1wwXEwYjB~_w~ zTjoeCB0Nwb-#%*~O+T%;|I^~h)`iI2C`O|O`!}XB?=UWGX6fmb){AcuII@^1q@`qi z8i^xMHJr*(K~v?s;!XF~$kXKc(=bMP)e`yK30M60Gs>DQ$#*NSnEHJP5_`)tRPItp z(>j~wvpGzG!nM>bpGSP8H}zMe>Fy~cFN7_P1()65BZc3sI|5{apw1Y! zTmO#>^OSVadA8gMJF#0AEK|}$=K|XpP88<(pXBpkEuQ_sW6BiBT*cCvfHi@8%N;k5 z>eXh8ghR}R?TeL~U);={-T9teGkFbjKR>Qi{a&7~zO`L2`Rs%5=Z3c34!@eAN7jj@ z)unr>``XLw@t=yXEU7<@x4cVaFIL%G{BX{OeS&9vGLCfnem$<-tPPb*{<;{C{DHY< zNwmF)bJJCwJIam8D=)G?mz-4&{JI`<7R}T2M7lLQO6pMGDR!1UpGeEOSaL?7t0KED=-Dh#6PX;lz&4_66b0YE=oUjp8Z?#4f?zF$IBLd zTKNzOo%la1du+>x=85M)^TZ<;8G3^#yjw1~8mefVM zrGgpvjXZ+MS4i(w#u!=8+F4mE@cz0&N`l`k;=Xogkagzw4wAHt$htfpo`D}cS6uJtDf8OZ&)BUa>B}VX?9F!=7&7lJ2=x2S zmOgm0^`MHlOElBDjx^?D%D}amEqXo0gc=?1k)IG_L~a zYE+u?lBFBRyY4rvSZRe%JJGYcpsFX_dk~oRA7MALv@@%s23<;_RFcYnNCBHF?`Auw zdcJfqzSwn4S}c7ac)WWJs~7t;_VfcNf2TltG0|iT;+vZBaTd%=r6<0zKwbB9cwa$Kp8#im4Nn}g#m{s1W zKc@+$+^hUX-UthkK1S5suO9}%Ruhhcfs)x6Jso{KVT z-aF*W6!>V161Gk{geh@VcbROmG5F&L&v8<=T>Edf_@OEfUMRM^M`WaF(3TyC;#bb1 z$(DTv`Gj)3A9$3jx=d2yOSh^FyKlwN52nAjB;IsQ3Dt|+f0~1-HEKP$*gSR3Xr@CY zJdJIGeRm2Kug~2e^Y_JwoVsh%rR-<5wz(Hg1kQqudqe3Jvfgdh`act zfux;jVzRCOCJM;B!Um|N-Nq1{A8*@?ztjYToUW5oe9W|(zr7SbD6}(8d+n6?) zn7{ZcLMnb*T%%fG;7x!g{Uy^otNR_uj2qfFyBz)VX8IB=d6*v4UPtw!7*dhkzR#R- zF66Osm*3?348lF!H(^6f?+NDf=}VFq3pW&+8nz^~DU#z{d3fbhbv5hpP$275XW=!3E!$1600^6SJ3eQ2rT;tR#IjV%&yE2W)Tb|CsUQ6>i!m$M{tMkuf z@X?mKYRnrZeXRKPvSxeyP5hrc9#+R!`XJGVN71T$jwr1=k)5HZ=bO6H#RlIqqoFo4 z3Hk^3JJQWbGs!3#sO)*^qscaeKXlnnS}IwqCz^{k%Jmlec2b7<)+?_`^k+&BF49}I=HU?y?o!} zzh~(3rMqrg^FQM_96FeJ#KQJ5OYfQ52)^Wxjpx+g+;9>x)@^| zQlGp|-=D*Z^kM}=Ah;`>`MEGnjld&!)Pmb6ve(ojMOvwb+9TrYranL1O;-+HsFE_s z;S)qA(dmeNy-Uht{PBKsgEr116%G+98}gL?ke{!0UTjX19(QE|sn6JzMLJhcQD*Nue2S6>3st@jrC(~2}I6vhVFY0iC~VtD%V{q=rx zqIL~copvLN9!LK;m(^t)af{`0ilIVeE)DW-4zB~h5*k^7QMKawQTWI*e{p9Fy63|p z%l~bWGp=H$z{bxnH%~;CozE;zzH&WLCB4bR+ffquLXoU#4nK8Q+b6 z)6g+(mrXp|DmJWD=;HN0JFg)srL)7X_@J=FmEd7@wt64keeHk-Gb_5l7s3{_B-I2? zvGyz5tExlwxcBF=NvJ=m*v|T7V ziTn+j=)Nh@2A384`}T`>e_@xBj4doKk&Y?s&`# zWkYznl5vUWIh*Sg=Tdw=PthkVJD7i{;z6~FQC!69@T~u)u*b23%k@P9Yod(2v#~GV zcF0VKKUVUCBEdb5dNPZ7Gug+yIQ;!y+!Ve3^vWf}nff2|532VV=x$T6Es1|$)M-M` zvgmSrVPP4#QbinOTj6fNFS*e*>L)9M!{#;E*J))y)V-|6%60>(z0Hpe+ovvRHf8Q$ zKhl7Z^ZQZ~6|r7>f9YFM-&S#=GP*bQ-eFP3O?ELA<7$mfs@pFJ$#f#)$a5O66-7n` zQkUPt(PP9)dy#GYkTda4rWL=*l>-%fnfIh2nIoSz?zMmXPN``4T7*zga7X{j+WuLT zzV=W8Hm9xLU#>S%H%kmE&r`kS6gcPP_GYSp-8lo3ia+XwrLR)ZZi>GxVT)8l}z@v+G89h5uCAC&cI#gey zF1LGAE_Ky@XBP;k?E8w=UE{q5IVX@G)P0eMX?{C>NkFUPfmc+;Fm5KA!|Wr6*JmH@ zHKcuSw<`^I&$nKuJQpqAd@o})*zjd;h+9@7(&mWF#lA?lh2aK|*fJ8eE#0N{)rx2Q z-l-eOFUEy%1AOYNrGJQ@U!(dm;!1#%<)Q4^`;a7ksC=fd7k|%gXaYq+n>bwbE+$F@ z>i48h`&LN{nomJ%=oJOpGO`}S@#0sTr@c7CqUO-dLzP)BuI>giI_95OPDDau_PR#M z>zNlkQTIlJq`YQuZ6GlZ<9Um`9A8L`vcSq%nAAtL(((b@HaR~`UXu@QUq71`KFjKT zA}WMC;-7`S*N0(~gDOWooVG>?a(B=FO+p)+(hhpwM2Cyzk&Hl;CzX zP0?-2d&Ka$)F=D|vDqXvS`{fg>#ritkDm|w?TqAxidV1NR@n$#Tbc72pc=nda>>M* zf_T{GJwJ|-Y&Qdc&F_q$`>X=|5y4Ck6$kjAaJI6q{p{O2uvl9E86m-R<+0$-P?ZAS zj`HxWyNs%Jw9f^9qb!W?p%Fw$u)4BPM6rH;BtgQ5{oGP2x~P7=)JCU;`PJoTnz%!| zgLa9K!ki9gi65g5zKzS59pVpM{y6eW*#p~vewS-WBm#>83+eLBLgtSl6<$}j^=O4& z36&JJc>0;|I~58H-FL>P44t0l3xC8p9yB&1(71pq{k5k#ylPwt-|=Y(35z6(egjoy zp@DLkW5th+b-S;&Zs+QS7X})d=M(7!7cdu}>5GnF$RrvSJ6+?H;;+n@A&kRWE!;Lu z#kN2fpfUW9YWOo%SruDF@V1K+YFW`!Ew&#XOfKB#GBH$*e0&-PRpolwRM=Z;MO#v-rxF-6DWBVrU4#>7}MiO;OLIjc-_!rMKWyS|o{LX9T+)Kup3fva&yi4EH z6+*V6WsV-ineEf0!d_bZ?Y}H5BZ`$@vuBy{A1vxri03l9?E5}8>)u|S)N5facyZaW zEOjhFY@NB*l>+O@yE?JwlXGR~dy*WD@x=rPS#82bC$?Ukd%Koxm8V`aj&H__7nCXe zOr8`2%}9<#088P~b=m#6?GX&Vz4w`_UuJr|l~Wd)i(W^PrC+WYvgvy4DHgC|k%dpD zp;dG9k&J&vRJ_ZJUzKQmjIRrJHZv2Cd6*EZ)cvP5)X zhOd{bP|8rEZ0mBzWUrZv9JK!;7R+ta_w6mKIl91-AJU|lIy?XFMFDd5ndecfwuwDMJ9%|u1+=RhTMq+hxr zcixCo*3W$5C~1$URJyqQP&F_t$ZgZha*bTO!6V~}HOrmH7Sop$RP25wp%)1;o4rzU zRWfK?6n%}FuI+_MtVp7ir}wMCVTiz-LJQ5BRr|&~$$^N}x~RZ(+m}8{Aq;dzHA-&G z3l9{hg|o!R=~=>L!=CD&t@5qNQCsWx{uEN~wzKW9b?k{HV^SnGbjOFqmEct)|G%2j zc(HVP3I~P-1LtA50Xu-(R!|Xkj!=8t4!{y~sc!~Tc z`}HIAq7;g#TQo8$@v*r)yqD|yc)C!q*E8C*&i~Nd3cO6(RVy`2zSY`xc{5^a*UNKr ztJ2J6<}DZ? zT$T2p{eF_3HIz}XWmg&2=vt|NI~P@S1Aov5Z=>r)QRo`SZ^<>qGM4vg<2duXwZfs= z{Kf})Ek7AY3tXDLY9@`Rt1?`1o!$;ec1G;68aM2T{`R_VcP$=LC^4%42d$Dk5;7+W zCnO*Pf-<@@L&6@G=8(T7K?Xi#KtJApE}5~RlPO&4{c(AGxepVN#||WNo`-57|My19 zg#7teew6)wTs{8#X3(Exb4cPjyz%Q0^M~rAWTNA8jxZ@X;3f7G$iRgR=*Jt$>MT?a zP24@Mt@bf?H-iY0AOqg?uR|?dqU~|H^~2_JHHh>icM{}{DKPVYQ(cd%V<)_YCxCnN z1h)_xOpTa(=^_~+IJYh)RE~(-aJ_0F1D<>4Pz#qh75S)TG%i(V($Fu-8FFT8OU6Ki++^oaNKl{_ge#SZU@Lr z04-c%!f|=QcXVqY)9xg|5KFpe8fq=2ar zb7N{+ZG&@%(?jKmxSRcVzGp%vJq8)@+>{J|wQvt0j?WxR!cq{D(hkCn$q1DrUT^hq zvW|DaYm^h*^+#&N+*#9Nm%zDmIH1N6aZ}8{j05*;q7bJK+~;s*m5+0$K5rud_iK|d zIbv?%uL8P%>ZzabYi?Xn&5ps{_@diPdYfo^~ z3;)%^mHIl)U1*JIVhWKaK?eNydj%#(%*}&iKLTX+p7iI@kcWW4%?2r>h4?<#EkX4K z=VpQ55%<9jMMJ1RhwL4&k2jKt^j|ISsR3I}sLD5}V}rkL%r{|r z#49d&18wb3hiQmsgEv3+j~p@jGMd0L6deyXc_3aKRD_vF#I7fWoM?NLu{X5$5AuShQMap^-Yg`R0G*h}l1V<*Wy1|9%%{9uYfj zbzllOyNCXv_|)Ldia@~v>fGVkFB-t~h}kVE{JEzNy*oU6?LTtF>_XB?M&Rso#xV1U z*qfbAsG)ui-ng^Lq4?AQF7x0KgDK>mO@a(~c1W#|KQ5a0vrbomsSamHW6Cxag4EPKG z0V+N64%XXHN`0G=5rCr`Dn^1McUApB>2&vv-Ca*nqR^N5J%m z*rjsbID&uQ7Z1gJr^a<_Wn`DZecjOWzs4bY#J-MHs`dRXP(SYzc_c^7-m$jJ1I|wH z5~fGQ9(z|>GZS*ZkDa*Z6QJT#1?Cw$7 z{ouZ?HVLLj#J6(n z9_WAmr>i-}@k{TXeI!TB4s=xozvoShV0uLB4zxjiaJX*p>;3~OJ~gfbT9JbLy6oaZ z+4BVq?Soq_v#wAw@Qsc{|9I}8$epGkrY z`0xAABR%5pJJ8Y#`UD=HJ+2OB{*;&&=s^R{&fWskBVq@d#$-a~Jq8)@#!H~$Q{y_I z3kx`VTI-?l)8aayISSD7`EIj^N84ibKmbE>plh<@WzLs;#1=~;3ze?pKtqh zX#BLe4mjZq{>;i@^GJ@E9XON+&VH~5(<5RBPSJt;y5Rjo@u_hgaL5eY*JU2S^oZ9T zIO_H10N06_*5vn*95FkPHV)3Ng^u(mp67>5|F9U^abE@Gd4sc0LB*%WbwC<4I6Erl zq4Cq=Iv^Pt+|O74BS*{*q&9*7mdV41nMcG9BqIT7MnK9DF66$4e!P(`5g;8)O;VCa VMFajl0`i|9!8 diff --git a/tests/test_fit_equivalence.py b/tests/test_fit_equivalence.py index 8b2f4d2..090147c 100644 --- a/tests/test_fit_equivalence.py +++ b/tests/test_fit_equivalence.py @@ -1,31 +1,40 @@ -"""Whole TFA fits against saved fits of commit 3298278 on well-conditioned problems. - -``tests/data/fits_3298278.npz`` holds the centers, widths, weights and cost of -`FIT_CASES` fitted by the code of commit 3298278 (written by -``scripts/make_fit_references.py``; provenance in the file). Each case is a -synthetic problem whose reference fit recovers the true centers (checked here), -so its optimum is well determined (one case, `LOCAL_OPTIMA`, ends at a local -optimum in both versions and is compared all the same). - -A fit stops when no center moved by more than ``threshold`` (1e-3) times the -length ``sqrt(max_sigma / 2)`` and no width by more than ``threshold`` -relatively in the last iteration, so two fits whose solver arithmetic differs by -rounding can end that far apart. Those are the tolerances; the measured -differences are printed (and recorded in docs/validation.rst). +"""Whole TFA fits by the current kernels against the 3298278 kernels, on the same +machine, on problems that are demonstrably well-conditioned. + +The oracle fit is ``htfa.tfa.TFA`` running on ``tests/reference_kernels_3298278.py`` +(`oracle_kernels`): ``htfa/tfa.py`` is unchanged since commit 3298278, so this +is the 3298278 fit, computed with this machine's BLAS. Both fits therefore see +the same platform rounding, and only the kernels differ. + +A case counts as well-conditioned when the oracle fit (a) recovers the true +centers to `RECOVERY_LENGTHS` lengths ``sqrt(max_sigma / 2)`` and (b) is stable: +the oracle fit on the data changed by a relative 1e-12 (`PERTURBATION`) agrees +with it within the stopping rule's tolerances below. Only such cases are +compared; the others (local optima, or optima that rounding moves) are counted +and printed. + +A fit stops when no center moved by more than ``threshold`` (1e-3) lengths and +no width by more than ``threshold`` relatively in the last iteration, so two +fits whose arithmetic differs by rounding can end that far apart. Those are the +tolerances for centers, widths and weights; the cost's is 1e-8. """ -import os +from contextlib import contextmanager import numpy as np -import pytest +import reference_kernels_3298278 as ref +import htfa.tfa as tfa_module from htfa._kernels import max_sigma from htfa._matching import match_factors from htfa.synthetic import make_tfa_data from htfa.tfa import TFA THRESHOLD = 1e-3 # TFA's default threshold -REFERENCE = os.path.join(os.path.dirname(__file__), "data", "fits_3298278.npz") +COST_RTOL = 1e-8 +RECOVERY_LENGTHS = 0.2 +PERTURBATION = 1e-12 +MIN_WELL_CONDITIONED = 30 GRIDS = ((20, 20), (30, 30), (25, 15), (10, 10, 10), (12, 12, 12), (14, 10, 8)) FIT_CASES = [ dict(grid_shape=g, n_samples=40, K=K, noise=noise, random_state=100 + i) @@ -33,16 +42,39 @@ (g, K, noise) for g in GRIDS for K in (2, 4, 6) for noise in (0.1, 0.5) ) ] - - -# Reference fits that ended at a local optimum (max center error 1.75 -# lengths); they are compared all the same, but are not well-conditioned. -LOCAL_OPTIMA = {23} # (10, 10, 10), K = 6, noise 0.5 - - -def fit_case(case: dict) -> dict[str, np.ndarray]: +# The kernel names htfa.tfa imports. +KERNEL_NAMES = ( + "Prior", + "ProjectedProblem", + "ResidualProblem", + "bounds", + "factors", + "init_centers_widths", + "max_sigma", + "ridge_weights", + "solve_centers_widths", +) + + +@contextmanager +def oracle_kernels(): + """``htfa.tfa`` with the 3298278 kernels in place of the current ones.""" + saved = {name: getattr(tfa_module, name) for name in KERNEL_NAMES} + try: + for name in KERNEL_NAMES: + setattr(tfa_module, name, getattr(ref, name)) + yield + finally: + for name, value in saved.items(): + setattr(tfa_module, name, value) + + +def fit_case(case: dict, perturbation: float = 0.0) -> dict[str, np.ndarray]: """TFA with its defaults on the synthetic problem of ``case``.""" X, R, truth = make_tfa_data(**case) + if perturbation: + rng = np.random.default_rng(case["random_state"] + 1) + X = X * (1 + perturbation * rng.standard_normal(X.shape)) m = TFA(K=case["K"], random_state=case["random_state"]).fit(X, R) return dict( centers=m.centers_, @@ -54,58 +86,56 @@ def fit_case(case: dict) -> dict[str, np.ndarray]: ) -@pytest.fixture(scope="module") -def reference(): - return np.load(REFERENCE) +def _differences(a: dict, b: dict) -> dict[str, float]: + """Differences of fit ``b`` from fit ``a`` (factors matched by center).""" + perm = match_factors(a["centers"], b["centers"]) + C, w, W = b["centers"][perm], b["widths"][perm], b["weights"][perm] + return dict( + centers=float(np.max(np.abs(C - a["centers"]))) / float(a["length"]), + widths=float(np.max(np.abs(w / a["widths"] - 1))), + weights=float(np.max(np.abs(W - a["weights"])) / np.max(np.abs(a["weights"]))), + cost=float(abs(b["cost"] - a["cost"]) / a["cost"]), + ) -def _differences(case_index: int, ref) -> dict[str, float]: - got = fit_case(FIT_CASES[case_index]) - key = f"case{case_index}_" - C_ref, w_ref = ref[key + "centers"], ref[key + "widths"] - perm = match_factors(C_ref, got["centers"]) - C, w, W = got["centers"][perm], got["widths"][perm], got["weights"][perm] - W_ref = ref[key + "weights"] - length = float(got["length"]) - return dict( - centers=float(np.max(np.abs(C - C_ref))) / length, - widths=float(np.max(np.abs(w / w_ref - 1))), - weights=float(np.max(np.abs(W - W_ref)) / np.max(np.abs(W_ref))), - cost=float(abs(got["cost"] - ref[key + "cost"]) / ref[key + "cost"]), - recovery=float( - np.max( - np.abs( - C_ref[match_factors(got["true_centers"], C_ref)] - - got["true_centers"] - ) - ) - ) - / length, +def _within(d: dict[str, float]) -> bool: + return ( + d["centers"] <= THRESHOLD + and d["widths"] <= THRESHOLD + and d["weights"] <= THRESHOLD + and d["cost"] <= COST_RTOL ) -def test_reference_provenance(reference): - """The reference file says where it came from.""" - meta = str(reference["provenance"]) - print(meta) - assert "3298278" in meta - assert len([k for k in reference.files if k.endswith("_centers")]) == len(FIT_CASES) +def _recovery(fit: dict) -> float: + perm = match_factors(fit["true_centers"], fit["centers"]) + err = np.max(np.abs(fit["centers"][perm] - fit["true_centers"])) + return float(err) / float(fit["length"]) -def test_fits_match_3298278_on_well_conditioned_problems(reference): - """36 fits (2D and 3D grids, K = 2, 4, 6, noise 0.1 and 0.5): centers within +def test_fits_match_the_3298278_kernels_on_well_conditioned_problems(): + """36 candidate problems (2D and 3D grids, K = 2, 4, 6, noise 0.1 and 0.5): + on every well-conditioned one (see the module docstring), centers within THRESHOLD lengths, widths and weights within THRESHOLD (relative), cost - within 1e-8 (relative) of the 3298278 fits. All reference fits but those in - LOCAL_OPTIMA recover the true centers to 0.2 lengths (a well-determined - optimum).""" - rows = [_differences(i, reference) for i in range(len(FIT_CASES))] - for name in ("centers", "widths", "weights", "cost", "recovery"): - v = np.array([r[name] for r in rows]) - print(f"{name:9s} median {np.median(v):.2e} max {v.max():.2e}") - for i, r in enumerate(rows): - msg = f"case {i}: {FIT_CASES[i]}: {r}" - assert (r["recovery"] < 0.2) != (i in LOCAL_OPTIMA), msg - assert r["centers"] <= THRESHOLD, msg - assert r["widths"] <= THRESHOLD, msg - assert r["weights"] <= THRESHOLD, msg - assert r["cost"] <= 1e-8, msg + within COST_RTOL of the oracle fit. At least MIN_WELL_CONDITIONED of the + candidates must qualify.""" + rows, excluded = [], [] + for i, case in enumerate(FIT_CASES): + with oracle_kernels(): + old = fit_case(case) + old_perturbed = fit_case(case, PERTURBATION) + recovery = _recovery(old) + stability = _differences(old, old_perturbed) + if recovery >= RECOVERY_LENGTHS or not _within(stability): + excluded.append((i, case, f"recovery {recovery:.3f}", stability)) + continue + rows.append((i, case, _differences(old, fit_case(case)))) + for i, case, rec, stab in excluded: + print(f"excluded case {i} {case}: {rec}, perturbed oracle {stab}") + print(f"{len(rows)} of {len(FIT_CASES)} cases well-conditioned") + for name in ("centers", "widths", "weights", "cost"): + v = np.array([d[name] for *_, d in rows]) + print(f"{name:8s} median {np.median(v):.2e} max {v.max():.2e}") + assert len(rows) >= MIN_WELL_CONDITIONED + for i, case, d in rows: + assert _within(d), f"case {i}: {case}: {d}" diff --git a/tests/test_tfa.py b/tests/test_tfa.py index 96a3c85..0c7b82f 100644 --- a/tests/test_tfa.py +++ b/tests/test_tfa.py @@ -816,10 +816,16 @@ def test_ols_fit_with_a_factor_that_vanishes_finishes(): """With 'ols' and a small lower_ratio a width can shrink until its factor is 0 at every voxel; the full-data weights then had a singular gram and ridge_weights raised LinAlgError (seeds 0, 3, 4, 18 of this setup). Now - the fit finishes, the vanished factor gets a weight of about 0 (the minimum-norm - weights the varpro objective uses; measured 8.6e-10 of the largest weight) - and cost_ is finite.""" - from htfa._kernels import gram_solve + the fit finishes with the minimum-norm weights the varpro objective uses + and a finite cost_. The vanished factor's weight is 0 only up to rounding: + the eigensolver mixes the zero column's direction into the eigenvector of + the smallest kept eigenvalue by about eps times the kept condition number + (eigenvector perturbation with that eigenvalue as the gap), and that + eigenvector carries weights of size max|W|; summed over the K - 1 kept + directions, |w| <= (K - 1) eps cond max|W|. Measured: 8.6e-10 of the + largest weight (macOS, Accelerate) and 1.2e-8 (macOS runners with Python + 3.10 and 3.11), with eps cond = 1.4e-8 here.""" + from htfa._kernels import GRAM_RCOND, gram_solve rng = np.random.default_rng(3) R, X = rng.uniform(0, 100, (30, 3)), rng.normal(size=(30, 10)) @@ -835,7 +841,13 @@ def test_ols_fit_with_a_factor_that_vanishes_finishes(): f"vanished factors {np.flatnonzero(vanished)}, widths {m.widths_}, " f"cost {m.cost_:.6g}, weight rows max {np.abs(m.weights_).max(1)}" ) + vals = np.linalg.eigvalsh(F.T @ F) + kept = vals[vals > GRAM_RCOND * vals[-1]] + cond = float(vals[-1] / kept.min()) + bound = (len(kept) - 1) * np.finfo(np.float64).eps * cond + ratio = np.abs(m.weights_[vanished]).max() / np.abs(m.weights_).max() + print(f"vanished weight / max weight {ratio:.2e}, bound {bound:.2e}") assert vanished.any() - assert np.abs(m.weights_[vanished]).max() <= 1e-8 * np.abs(m.weights_).max() + assert ratio <= bound np.testing.assert_allclose(m.weights_, W, rtol=1e-12, atol=0) assert np.isfinite(m.cost_)