From 9e58542fa30ab7764a463f6481cb46abd5c44600 Mon Sep 17 00:00:00 2001 From: daharoni Date: Fri, 20 Mar 2026 22:51:49 -0700 Subject: [PATCH] refactor: cleanup, dedup, and optimize Rust solver crate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Correctness: - Extract clamp_tau_rise() guard and apply in BandedAR2, fixing a case where banded and FFT engines disagreed for tau_rise ≈ tau_decay Deduplication: - fft.rs: extract convolve_impl() shared by forward/adjoint (~85% identical) - banded.rs: collapse update() into *self = Self::new(...) - lib.rs+fista.rs: extract compute_raw_baseline()+update_baseline_ema() - kernel_est.rs: extract adjoint_spikes_kernel() (was in 2 places) - biexp_fit.rs: extract golden_bracket() for 4 near-identical phases Performance: - banded.rs: fuse normalize pass (eliminate O(N) second pass on hot path) - kernel.rs: compute_lipschitz iterates half frequencies + pre-casts kernel - lib.rs: skip banded.update() in FFT mode (~4k trig evals per param change) - kernel_est.rs: hoist z and v_f32 allocations outside loops - baseline.rs: cache MSB in FenwickTree (avoid recompute per kth() call) - filter.rs: extract forward_fft_and_cache_power(); iterator-zip loops - peak_seed.rs: combined median_and_mad() saves 2 allocs + 1 sort per trace - threshold.rs: move s_bin ownership instead of clone; fuse PVE loops Cleanup: - biexp_fit.rs: BiexpResult::sentinel(), has_fast_component(), derive Clone - biexp_fit.rs: fix stale comment (0.2 → 0.15) - threshold.rs: sort_unstable_by with total_cmp - indeca.rs: apply filter directly instead of throwaway 1-iter FISTA - indeca.rs: extract interior_slice() helper - lib.rs: extract state_byte_len() for export/load Tests: 109 pass (3 new filter-path coverage tests for indeca.rs) Co-Authored-By: Claude Opus 4.6 (1M context) --- crates/solver/pkg/calab_solver_bg.wasm | Bin 416734 -> 407544 bytes crates/solver/src/banded.rs | 41 ++--- crates/solver/src/baseline.rs | 20 +- crates/solver/src/biexp_fit.rs | 245 ++++++++++--------------- crates/solver/src/fft.rs | 60 +++--- crates/solver/src/filter.rs | 70 ++++--- crates/solver/src/fista.rs | 17 +- crates/solver/src/indeca.rs | 204 ++++++++++++++++---- crates/solver/src/js_indeca.rs | 3 +- crates/solver/src/kernel.rs | 37 ++-- crates/solver/src/kernel_est.rs | 76 ++++---- crates/solver/src/lib.rs | 76 +++++--- crates/solver/src/peak_seed.rs | 40 +++- crates/solver/src/threshold.rs | 39 ++-- 14 files changed, 509 insertions(+), 419 deletions(-) diff --git a/crates/solver/pkg/calab_solver_bg.wasm b/crates/solver/pkg/calab_solver_bg.wasm index 54fec6b1df63e563054a3385f9f48ff63e47f6b4..828044928a9639afbc5193364253e6582b79d756 100644 GIT binary patch delta 59963 zcmcG%3!GI|+4z6<+WXwk%$&IoGhEi*fR3Pu+(8Al0s_EG#QiOf)MqN_?{t{e7Rc_daJpt$v^X?~l&ed)?Q1 z)^l6Wde+*{^$mFQPi=+7qRT&QB@zki9d*7OE{#bz_V4Xi z?2J-YyJhpwvg&iT-H@nJN&?9gNfK>$*tSwm(n`vc>{ezV1(c(bgDjP??W9skOQ}pE zrL4g!-RPu-XxGX%5uZt`v}#UgvT2L*DLbL;q*beGy>8{CviUS+hdNd*IhCCuDchiu zPGg#Oi8ScjA+CC80CG%b}TGR^54t8>KR`Amuog>=$sZEb6_ zZHG3K&Zwl5I&yS5NM2`MU)HI?)|2IFDc4$AUX&VOeYt#9s$<@f4gYb;lBJ6mExi2F zi?3X|@QTY8Enb@VVD9pzOP53eztJRJws6r$K63GYExhnQF1}*vWr^3UVidjj(o3#f zdf7!wmn7~gznD6*ZG(ynT(tCxg&$mW`O=I3U9 z{o=!pEnIrZq7Pn7XYMR-%nY%-9e`J40Jz#IO ze{J1lZL#jK9<)o_?0d=e5V?M3O*{FNwd%|EF6)X*F5aoO*w5L&v_8pANM^+mOV3h%N0VExrP?ML<#_T%=G>Jhs59qaej|5_{P?ql|&_D*}Z z{Xf>XtQ)Nx=+0C23-(8^{ErX)#{Qf=_t;~8V!v3j?^h41zgvH^_Sm1bryf1!r}j(s zeQKL}%le14%YNCu-kvaa%x~>4*hd`J`7@~Y0ril2+uCnEW&h6pyggJ8`8mYAMSWkb zQ14iOu(zrosQ*^is%Pw1?eD7XYNh%FBmSd(k9t^rN_|rO!hX&Ep8BD>PJLSa$=;+& zkEqY8&!}J8d+k+rU2V;;>_6LmYKK~-KBu0w|6<>(9#vmZpI5)OU$-}_$JAi)&8cPdrAGy`ki`N{Z_rA{=oBP^^$r~y{huxw#)W6)f)Rb^>^}r$^N4Ky1H9! zwQjTjsGd=eTYt7T+kyR4HTEWZqy0DaCYc^r_glU8ckG+(=ha`V@7kZX-%+2mo>1%T zd+Z-t1uK)gdu!_}_Um6L|2p4lef*v`^025o%YUmqeoATcB$d#PpXt_zXyseoTh@8Q zi*c;etq(0;Yi;%pTcr=vc1M`E(5;6TuPuM9u4_b?F(={vxtzS;<~zkZZO^G`@g2Xd zcz5}=xc0q#`6a3{S|RP zmpG+!#eeYPzmkxS-}8;eW30AKzu1;249YWrR}#g$ZBkcSCiU&6$E~)#&yYG$Qde4G z>RVyi-?V(oYP(BPJ0$h3sO+k+>>X`8t+rQwPHMZPu8PV&ZpxNdxBsMVgCq`+#K)te zd%}u8-f@Q2w)jw|ETAz8E#FqXQcgk?e zm}KU*yjCC6vPzN~C5b5(0F5)HUe|ZF z?D#234U#l2Ogh{oeaBr@Tq#NQc_~;lktcaWTekg#$fD#ghWUq?q)&U7w!HNeNp+Gm zEKF)LNq-zVzGaIf)k;!Rm^9ENoiX&(mVG-(s*$9DVN$^)Jvg+jWuqh&B&l$4g|9wI zQodh>l~Of#NK!5;RX%UnNVWa3@@I#w2dzg9Uop6}*&QLcYHfB81D~?pUS^h4(n-J9 z>Mnk{{KW8MI}XT_*EO@8>6P?D4w?2()2}^b%aCgMRqxlvzytHVyODpM;o3u%lrQh} zC*EQgR30KyzN53geTAort8Fh|V*$}Bv5XdoLgX`@>A-R{RQ`v~52-yrE}wSTSFJnm z*>~73%i3LjbmSlFgtl~}Z&8tNfpQNYemf~WM|{lsM|sZ?vjE}9QC~W!YqqZK51{8p zmB@Vlksmq`Oz+AyuYBa_B~1*#lAfj*fE^8>{PoeZ`XyL8V>(n^UAB&zJ~9+=6NdJSblcg3F_G=%EQO|)Al`g0PB)& zBy^Un8zk#40WIA&{e!vLVMh~w!e1Fyb^G}Hra$wSZ3#X2xOeIn)Fj&dSHjTKfB0#_ zgx`!(zx)lygs*$up5aPe2= zV<+k`MR!Z&X!)6BQZAeLkp?LvF-7@lGtv7d9#(tZKH%BmZ{ksWTa`yn8WTrth@yHX zjj&dgzclG1)>q5v$;z49|r`*DC_UNI_ zOea9@l`Lz8Gi40d_mAF2UUkfk^`*^gLg9e`GIU%1+A%|F>1Wd4ndsJ@ZZFE@ zbr5uz(_|}^=N;SC^sd6z=2hjBj`a%3yuV3iYduYtubtCc{_IH&qhJ3t&@_C&-w;ow z|LQ{weUQIRM9|wzVe#t}f#MvIO;b;ym!F#Y!B#<(-3eG?OlC!bpxirCKi+cF6Y!%H zxInRR@{|6KMER($#(UOHd(DCXZ}4laqe4os4J9aFUE$e=lxNly@;ce6ZnE5zJl9yR z?T0*i=W(Z~U5}TC9=|ZP?*+Q|^YV4aPiyH#fco3-stjf2Z%F!O`P*}%qz&JwB&`XP z#>_Y=wLxB=E~n=lUw(8(KEKHX&2Y=l&p3SkO?OlZZI6?5K0;=HmqdL{W^YgBZN&<@ zG1TAqEvZ*XA;x$hJjo3H!!wU%@a1pLdgK_|6SzW(@0A=K~vW+lf4lcRjr z-1Aba<@M*~`{uSW+NbB9qEyu zqsDFidZo@?lHb5GLdsR&Cqj@i#?mOFyz?Z#Ca%))CsE_z`CnAm{pgoMNZS$WK8K^`5&{O`B|OxaS+I{K@w1R{ow({%){?A@0CI;2%dZXj z^N^V}2(VfT4bhgaeJhus7s?$WVhv__vDX|IZuZ~dpIz#r=v%y{yDqDN$-6_OG|xIy zw$tqfTIrU~%X5TEcwNG<6k)5|YnVsp zEVowI`AGpyJ1_P!kJ0Pe4zIR{JZ=+}hi|pIi>}SrEpDAXo!I&N^Hf#m7rSW{ZK}s> z8eEge=V`T8w{}Z|CTFcK&Z@B~&{h8IBZEq&NA!fCdt%hOGp|PdTdkk0Poxj0?k|X) zM)7C!=kV9BUnlER>DxaneaimdTKiwN_V2O!^;GHk!L@r-%h>r*pXsv9X9FPL_#h~| zQBbW`xAhRQtrB6zR9?2DG1H|_)V25MdF53*bjjq}S;>`R+-|FG)2SY}jalGPuUqs` znc<>@in_t8)#?NqWoDH$40uCf3(zADtW#UoxgSK;$#fUb3$Zu)E$O^2T;?6?TPFxs zmOz;ig0N^sQxcy_@w_hCts6($8y)%GVt#ix#jjerQ6hSilkHcCq-M_&K`{f^nSc6& z#S68?+^OPx@)o8f?j)biPf6V2XsdWeyKd}MD^i}4=f96aptSGr%0%eUfm{XJBA%tu#a|`T*oaq9^ zcOdG7TPO2hg<$#GA2pWd!v_r?i@{Y_1=n8>1Q)>rz$M}S;F7R_%k<_za7Fy6plbzO zwN-GHfR}))wgYx0;Hv#+;QDI>uJCCW|Ek;*V{p~>2UlGcTy;H$-@zUPT&6?>t~$U~ zKK8kE=|s5R)?`F!Gkjx5Y471ix z#2^|JOd*tAgZ%_>2h2r5;v(W^u!tKg?;fsmGvUF9Tl?$XQ>m|kzj}kZSzRw9Frwjj z5Pquey6Hx*>pm}E{_B&Y=H-{o6k5=QWe>PrLf~C*Szcb^QpEBHow2gH(uJO2ax^P~ zB&F*whR;EObwXykZhi+1e;~Laa2iXekZbueuSJl7iTJ?&Z&}`aed_)xAbba=$YoS8 z7-Z$0jtOYqvL^^rCK0FRycHC zx6BgC7>smi>E@@ve+^w7BuioNh#cj`PmTLPL$})iAd#u{I%l{$xF@s8OoZHBgs$bU z{vHy>EVgo$DI=p5Qc;dtnRBx?99;g{Q=Qh}@;gt>D@sLzX9i-mdT{xJKN(l@YTd(y z@vlyG3nPbr!{iM=_P~Hc`O9eqC-|6f#q-O5b?o8} z_jGD;=IBH9fF8d%=Q4MyPWxjnfwQI=z?h3q=oYnLT1gLQL`o0u>ToAA$*3=rC`i(D z9`-VKD36GPj+BHRcd|Y_MD!Efrl|GZ(+~hDa2|HQWxacDD8O7l?%X@>656Pp?;5D=|f7f82A5gM7XXovfGsXbc(JV8p?1{ z2|kkoK*!D0gJ+`clrFyv0zfe(4vBpqLR=e_*qiBD0P1CfmCmJuTw=sGSgxun$o|{4vFG^4h81O6v5S z8fE3nxBkp6Nerq+y*#CtsHuiaxdSC;OVY5z`jCKziKHOs`o+GYsMz z7!IWIh{Vd2dwkfEJHv>qdRW$(+1S%VfL&?kccv;3`=`5NMHW8i=0$xt@SiK95G9&;8S5-Pdb@+E>pOrc*d&t5rOh@^u>M#r?w;~9$IueAD-hd$NDj)S+ zqt#WO{!D9W_x&*otcMMSgxeL;En*U=DPg4ObQlU{G9?^uPFN&feD*+ZPC3F`Qsv`8jqlpC)+`b`N1T<5B zE=HRWznOi1Utgum(vU^9wR-w2lUauCWre1*OnQaP;Ty8(qN}wc4rIMpFJisEB#P`2 zFwf!o2v}zu?G5a9!M2FG55!yvc8JhFl0OS;5Ls%wMJzGKNKU`}t5N4eMUm~>_L^9! znruL!Ny{1!P@={o+Muy&)? zZ-SwS0kzTJgq#8P*cqMTpG?}8l-rc2IU|Toqcjvwk7>0*fYktfBBmhQlq=5?&AUfR znK?B!1Yr5iqA{S;hHVsAv)Lq~*DCHrmqC@(`4BXW-^dwNwZ+Nze#s5FgBr9(eXOJ&XM>2x3mNZ1CLmZKAGu>Gd8qBb~ zI!P!v(l_O=5Vk7IQ$GPyQvTkUl9-cX@gF78ku3OWy|HjKgTKHxVw26_l#W>?>Ok|% znp(p)dDpz>{V{X~`98nH?T|_A@@5K7@z%2i8I5%To(M+^1@!_^NrdZ$RTV^=DBttz zrp#pDmgyW&e)8Aj2ajYBOqv-)UJn(vkF6qvtR7nSp6j&EC_nXkUG90VrL;gFS(C&R z8H@gD?gBHKMn4ha@NuFfudpV9Kn%n%-wrX#XZ?&=HXEQD!InyX-4Nou?huqWBdO3? zRxi-CTh7zG>UkmJ zny!a*8sg+XEmsl+@jVNgB`PxuD>mG&QGgHQX{b>#Mt2BClMY3O*xqUuw!i_Sh`}!C zv(>GWiG}Ry=!rh}LRr5wk!onS+t5q|+y+BSUQsF&UzZ7l5q; z+#cCZ6CzAfrUgQ3@~>Ie<7SO8hA0J00Pftf|AJ?oTVC?Qh|yVD7x#EAW-5!iMa*J( zH!IkDE)bTTt!`6lDRV)%{J;xsB>*7vm<`cEH;Ay1oMPP1=xk`{0pI-lFyV(bCxd9R z+&Bo$TJzL;uEDX8o%t;d_f7?hqezr!1&+Zb%Od$gXUOB3}G{5KK#IksE3 zb-SIKFmS;qhy&pQjL zTEOX^ZI~B8K6`Vh)Jr&AUWRhM88(??irbBT&LUH!JnP0;C|tp&Gwfn{^iNXd55L$x zxZSHYs(wAHKG@a_sWO~8VZb?ES3d9et>`BAzc{3%b+ZVV6W4@hQuL5nbU-8A=dF<> zhol=P@u>Ihn1;tu<7`96b@~J`bfSmolh)9N5fq-3y-^p#dMPY207+hy4AO&HRZc%S zA9Ft-lmy0O=6iFYH?^FxAFYo;(X2vnf$pRdkrfzZ%eTGWQnIJ()6C*yt9MSXdz!E~ z;bGnxy|B{L^%?q{b-k?alDwV4TQcr}(Hw@~K@@nm=r|I+7i&b+pAbPff+DUzyRCu~6m3J3Z ziqVNIQ`o{t+CpBaN6)Hp=z#!iD`n0o|L&zjOJjgG8|q^1V7FVJY5H=Od$vB48sHBw zUOzimhzH9C1lFME%<_sIh6zI{4L~UgqMhxYrHf*{$?LQAOzePSN)t1InJS&9&GJ^w zO4MWYy-9ucj4$6T|4NTeTbK71^ipcUT77Qu4dq?9HXOu?+0)m%7fH|T=|YwKjsXQ~ z{V^wYxRzmMD1xA$uI*zQL`cw?3-v{Sx7VveVzB|dS`PKgioVUYsz5ru{Os>@XPrs) z9o5Jx?A>#OG`vpXNDCkq&91Nr;h0S7Ad@1IA#Hm|Auxscj_D=1`wY8Lh@})6&ms-E zq8cM%It>AJ>W&`hvl(>&0)pvNwH9q!Dm1D*KM62gu)XU7Ga{O^O1 zJO=hoWWLT2ygpN(389?dk5FnVBZAS<)@kdh)#MZkIDfr&T4fp+*Qj=V8oADZI!e@a zj>#b;#8ys~genE13iI?^U;q@@>1lv13nItk){MJa;xqAl!T#io z($M2>EUdltk)xMjJ&A_5cN#S%*ROTY5;By{vb$m;HOd8~w{vT(W{PC4&~UiA@nmzZ z6%Jxw!CDG;5nWG&HVZrnLEnf@7HEjhlB~^-1h_2zZn5Je(HWL3G4&KeEA3^Gsh5ga zc*gnj`gnJgp5Y#*>)k=@fVzm~#qf}Y`V991?752wFCncU>|<$VwIYn>^e~i@H3(YP z>3V7T%)Nt3QM8Gb1+T~ta${Rkf(AU%!{mV#o$; zm=3eg!0Nq?PGWzDSlUT$rMReB4WUu~O0SelGN^Bv6u=cPP_N^v$ zT9ZxWb=t`5x;ZszQ3O`Vjw+ezM2l`=l&=h%58_15{7aFZvb-V1^vz^uB%AS_l%wY?wo6UP81W~a?2)-Xk;VjK! zZjxzdx>`lS)T5Y@#pM;R4=%J}75X8_A~Fr)4_8w9 z%7o2l4@Y|-MJwznGT66(I6_c^0Hy*)%58!ph&xBP;1rX@ppw5^92R1ize_NRP$gA2 z@s}2?hJ~br0cJe`A^Us%#>Hlh&R*=*vn5~@K3rD5$WqT>kZfLJ@y63Zra;ijnd7>8 zfT&y%KF~yMopwc%pp?fmH;&l}x48uuswShXWHZhMdKO}Ad5$?umi*6O*LZ`ub7QlIuaGrJuqy?^o>YdjKJ!pcO!~+BB zI3CRt)My@dMGkrdDKcO09JVrjEkJTgC z+ZKgn2yd1AjW|Sg4@-y>G}cw)iR$Qq@rx%feH}+ypR=r|`o2Ee`s|R}+r^!5sJ}u{ zEdP9X?4|E zxA*ON-g==uweq(h;4k}r^mnVK(EiFz4*?;Z_I-Kt0lB~LWzu(8ed9XR$$IP`hB;m- zwLE#1x*Df04MQ5nOz6}ERp61EpfV9NmW}6}aN(xlh_sqEX);^O?qr=4)XL#+iVm3( zI-V^5?XSsXbNNC|ZC!%%P4Xuj(IzR_kXEM-Z_=3ws)-iV1l0yf;E7Gm+L@pR5^@O9 zY|xrfBT6=V@UnvjA;BsUBfMl7hp~?$g221HYk>xNI?JPe0yYGLTar1PeH_8Wnu?V) zrÐvsQ95V0^Rd=nNE{<%|Wm-ieP%%S;5W<754Fussb-cw=v#5i16ZFB9r=bXu%Z z?zrIB8P!&3gXy!T!hc4$1y)vdSVMzBS=D~}*!@j_XdJ&&Xrz(mux{)>kNuK`ddmJk zTL9{H%R73%o7l{@6i;ij_-52&9(tc#0O#o7^I0|ge4}gO042mIHM1`7W-K{fKMN1= zp=LM?fa*d8RgYV%$F4W$F>=7Kk(pybgGRd6Ov{7JfY{BNSCiH)U_cL%Q*}Kmvd-|c=0q@-_9mMo2K?xO&Avg}TRYNZ2d{*T-JzV$@$2!8+he^d1 zpbs;M*sqCRQx4Ota4*PLxL00}z0Dg7wa*LJ0a)JR4)2}jj+M=#b~wu3Z$}D0BglmU zH^g{|fUz8eIAHzVbxMzF0lC2lluC9T1tZv>6fOfG!3kJOL@FUk1n<K5$!k5Gv>BIn2gMVq0LYzD~AE%`YKBZD9x-o^w>yI+B*#c8zYbCr)6f1 z@XUyHwuG#cP$JAwGQ=*#Q9F}f*PU<~%4Ww*FCmlGg2K}g`T!(|MwN_%(kAqZgH+aZ zgmZ< z2C={A3aJyt(pIeMI@ngPnoGlgS;PtsVOc`S(LAmZ-Bq^cfI99cBQVo7xkpzZ10V7x z5_1e*A2sXTNxJxDZvx?CndSr)qHrAr&E+Zu+;PH2^1?>YGGQf&h@?x?W7_UV+3jR& z7s~E@Z^|C??@{)s_oVD}m9mpMCmzmV$a_-um{~PqK^gTgDZ4YK>`tLjdqmlS z_aMblf(Apu_W?9^NW?62#`FolN6ftxdIx7jNyr!)c7yZURblcW&_d`U1Tzv^7)q)v zoFt!d_z{MHc!)VuSy&L%9qdl4mT3nvmcfRx2sfcmaB)$ug%Sr-S1ol9a~&$fc`OpD z#{>^IsnL^%9mt^(>NFsXm6C>S#J#Cyu@Y85Z?@`*!Kh|6vc!Vj(T9eZZr1qWYjS2y z6AT(uL=j}q3wfHbFc2bYm7Tm5^1mVDbf~Zl8wxbyZ zh=(S0I3f)6f%De~uQaQvxHTbh8b#u?k!g6|aMGMULUv1%$Rki4uP6W^4AW={JlE2h zj6|6BaukHkKt-`W(5XjA?smy4JxcQ8js?0&y5Ji0W&cK)Q#um$!8ck|M{)|v-<067 z7Bwn$?Eb%6=qhiwsFuiKi#uWsg4wH~+Kd?!j2&^XBt+z9WY?a= z9KeCTXLQxt{!n)z3*;lPh8qL(3Zy{_^o68R7Qi9VV=|(lnHTv*liR=dORcJ@R2}4S zYC#f+!~nSd?r+5sBg2VZ76S*fb&L{Kv0_l)%wTvR4LuHA%?Kz&os*}7kmtEJnm}|7 zCNTJLo640&96Wik48VKM2!nnbVKwMDc)&2+yjx#gOOK@o(p)7*gvNM7jc5}&J$|z{ zO$ah4I)TVfV47h*rf&7dpQp!Te;W%AowikvU(3AJ5UycX5wb=#I?PY9!WHWZ&H6iI z&`@F$aw=A}d)HyV&+Bn9w|%fzkJICVg#)0yaqHchhzGG85PGkHPfunr8AJJ=A=t?% zZbsSmro^-_Q}q`M<$JH?O#uN97$C|4DeeKI7R7w5(aNb!W?d?rx=?D2Lkkb8^q}(S zKL@oGoY1buPH3TqD#t_mha9iQC}LYdYvFxc*RIoJd;9e=*x0Uyltd;LjDZbYHF^uD zLd@J&MgQP9CT(On%78f#V@!}u=u5I#Wx%Z=Ho|qomnMH;H%t7cZeOOiS2P`S-l@sQ zU>n7ka{lUm2}Mc3<4_8$=ogRcWhH*E%(AJ-0n;~1F5_+E=nm7apNV?nJ6M^7K%N?| zvL)*pR;5Lhvtfx?qQZrO+m2Sl)T-X#;iJ`YYVWN<_89dQbzN_8$1!S3VZ~P;N_3+( zSpIci4PHA&4I93?_n`!WcOJTqSxK%6!>BPRD2Wd-4L+4jJXZZpDjPFZ-4si<4fO)N!~c~~Qre^XeKn9cp`qFNo_Vu<#Y9OM)CC=-upcPm#Bk`=#33^j4bp$Q>tcgKhRCYcmPbK|e?qWvnyQh`JUC4qGg#V}#d^Y*ZecW0 z=ER{4&<8#PrYnDO++jJ1Q0Xw>JfOqar|EF2+F?hhOoy#%haD`ObXa0MC+aX=5(wQt zUCk?)&PW&NOn@EqSb^qgzFMLPZ}rt?315Dkx+)#_?~UVB*Mhr+74++l0)CtbraMlx zI~WB_cWnPns|rA4hBi@^4v|qhwB~qqUVUZca-b{eujvic4E2S%y#lHX2i_8>a);Zs=S#3i#XkuqzEyuXoG z|D7Syz*}D-)_`cAykYN1P)(H6U&YKDxJvXI2?#DoK%SM9TP3hMN?HBYQono~BZOmS;@9t0$akAl6k4*dt8#Xg-PZa=ks^; zPuV2_iIL&$>>qem0um!B&-4%Mk$}WV%I-LT*JZ+gRYDRafmiwm_DMiuBxP^^z*`cK z7)g1ve_-Wpq)3dUtmq$DB>{<%lr<8#t`vo?>z{P11Xf2W^7JKH3T=?UO;JFem27uO zpf^g9XC>Pv3EUB-$g`5|zS|xuF-K97JS+LOOJGZsBF{>;$0e{WN|9$J+b#+0h*IQP z$@Z)Sc19`ktYjQ8p`$mW z6nTc(k^uGA-eARS)m*ote-`L{!)$f9x}!Jv$!s+yEhGtD2KhOvT|vBu&QV9pbIBYv zdVml%tPK}BxEpRRJ)Rtl5uWisbHX})HAjt@CB%9dyg3JpDt7yChQWz&b{i3h)XDMw zFUb7OuLNCl)j2sBpB*y3VB=iXGEM}D@OZ=M7-@ykB|^q)qUhTS5EJ7Ij0#5g_qnR4 zF5*gJv$0|E!4+M~wblk->Qd9Cm2}j~?yhPp%CwT~*GkegS2CC=cCf-$kUY}Lh!fS= zL!^}o158Aqk~NNp6<0t-GVuDuiRy@ncSt~7N=cDTuP|AdDkiIWHNuR|V8ym2g`>KI zk@M8Z0r>HmWELLpP98NZ7P@d{z(p2?`#LOF;(bG7$eXL@sY44!D4Nazl#TP$gc)X< zL*fm)YuLX|;v`BpgneTyzzl+*Gxh1HK5CAniB1N0FGzOu%{&R)xY~YmF#lvVM&0`5 zzE7O23K{jx&A|<4s%EwO=HQMq)hxn&d(Y&$B6ZW3`!=7WZfz`Vu=4qD4YM<)F;^WK zTzj|O*yo?8>KjrVw9bF3H7I^weWEYVO-C(#_kUgOY_MMM`@;8By=@Ka>%Cu{nkkw6 zfWhH@K>m(Fc1mkzAI`DYCHUF3E$P>mDyn>6DM0 zdn?a62hFG3v&A4Grb)>oE@4S67bAoqM`7xoKxbq4N5Rr-N@lq1^3R+-MQ83BB;do&kB@!L>zi0(YZE-UWS%VtI$9AYD$8|BJ#KD}|GI z@Kboc7=bq$DhJsyg^4)~_qU`@jNnc+t=Nigm-m$l)9gUmS(vuQ>$uS!tOx3jGW!x8 zn$Qqwd>2=i1Z`*7gHIqX$~ovpSCeH>nIjN`G$B_;Vif#@a0cHNcO@6EwI;3k#SfC- zdfm=%i#god!Oj$Omjp_%EjH+66~YT=l$_wzbL@e^BWKuk4K2EX3ko=qjMO*TwD2346WG=ZD_q28r$o+U0rRqxVz0SX&@%_ z8qY>4nhsc8UepqHCygd!i$M~1Q;Flvr~{US?7<$8`}3NnXV9zeLD^RfHU< zOZ_oYxaWIeO`57k$o_!;SMmklo zKAX_lBNE(La{}aGszpoT~?RvmQ0vV$gUHc(hb7?8%(gsd9R`bi~DQD5k(wRO2zty#zgZZX}oh{p+w5vYy{$(fCPS`|=>3=_RYdP~2J zHsB#B6Kt2YuMv^ym^X7i1(|o5TLNVOx_Iu%**Q_7MK=~7QE3S)r9r}IOx2pe;y7R` zL;Nu%z$=a=P^bi=EM_MPBWgGXC1WVbPIuNH9om4?IP>j|U08iEiXrU24+E`uvH2_H z?Pq2YXHr&IQdsMfb7vf;!p?4f28!)c6ct<&4^vGOZSD86d+m=TO?Ugd?%tN@^tVZ= z+rwHhVf(#tEwnG;jbS+0|FF7D+%&P&w)$C?!CXlh-CRmki=QKk^Lo{&WS5|X1`ZYx1tW}L_@N0?UtY$T;-mfF7rqn-<;~T_r;b1jeonKEBCupiswSEIp zTt!)ps_`3%%Jz>c_)SD{;H;V}?>7^b?jMya`7Ok;T&iZv`mIEQa%S~L`kVeCH;X!AxZGf>lZCcs?hk{_v$0+BUW}643tSN&I*=3sxE5V zDEH#66(6>wk21-ei4R5bGr^um)fHss#!;4ycLs|d!yoF-;OfWJ(A0*BJ+GA>x zbxYs#k8zJSNxywuEpK`DRZzt9SNxG2=3+A#ZSfz$e?7s?y0`R=|FJsFvc3|0?n!lj z?Rv@D=|9ectikC!)i3g&;W`1Zi)F@Vf+Gps)ay+2nU8W>&y z2;KW|6uJEm#Xp|PVAwP2xR$-r&=G#GK;6J6gzoZZ)YL;L4R%OL**5_#WY>{Z!4;4K z&T{cg?0!a_%dLwtVJuiN6c}$sxfr(ZlwYVL23Z^WO8={-h$rFqlJ?&$YewIuwB6TW z{Uvy-ym1v+XfF{%5v*Y4&=$ zqqFOa*Y^#dV?SUG-p3A#9<|iC?-u>WS-QdHu$ZymOGQ?feI7*@c42Q?9IWWFFY~Uq zC~@dgVYoF0!$3!PN=_e2&CR;5Xz~pWCZ1@|r|73ojEfEpesrRJE>WHHOw>d3?4iNM z^X#!i^;V*OWGxClf08|t$Gs=n(`fZ?C)ty&MSWxE+p{c!A3xb1Xx$J5C)*RPPXXWM3uOoqpNDXI8yY_F-P_BB@?i@~@Ws{d>a(ovquV`~%S@$&Q;~)SH_mdyQ0Y zhn;8JuEu{e(e2+ZDYbe)@oT{or`U(&?vZeo`!&B7>_5exZE&W>Um=dF+R4A%pn9&h zN4C1`D-2Pj>Bo?gQO->Pq()88VRcgBYCpyBLNl(-fcA!ulM%vFf ztj!wNxy>^Y<{xJ)h~4e1Cz=;-R^#_gXWhEZ;_fDHJDk#BYgcF#4<00~c(-6*Op1fM z2wpGm7kZAy$+&Sd%OHCT>n$?(kb+qNI-i5@s3j+nb>fJ`TS=jM)TN|rG@1L7t`l5+ zvuZ1$|5!*MeAfImeNbBD$IW?Epe9A0NqVYF27$66ni)~MMC2!mCkW9Ukf&(!L^(K# zCk!O4d;2|6-BEZ`6h``l8CIIc1N~*yT9Qj1!iw)OjYoQu^hBGZG_FyshYf&(nEM>YnlgGV$~XFgR{WH=_QDN}R>6?XlPC%)lnM&q$5N7&7Uuzqd*z3%=FqrNp>k18?(7;ee(fxEmKQnj~j4Hk5TTzs2wub7IV#s zgah-KyL=p7(@o|Y>B0f1y*kOEZLfWe*A{B>9DQcKM%vdUF$v(b8-fvsK;l*;H+vmx z=wQkWtbNUm-hi+a4pW?*5N#7I=N?2|=*GzA$OE4*SE=c?8JIEXq94l;nM0kXJ0(#? z2j~utON(=o%vf7oa~jQ-ihx%S6)6$q3jd58Q!TpSidZ+4Br5e+S0bRS9wK@x5wiQb zOQ0#6`+^hxo-i~&3hs*&V@%(4ASvvNu^}dncxuAwdnYM-s~|#f>^gO1E<}Qj%zChF zy&6A04Br=FflMMZ79z-?HDy+e#|yvV`Zxh-1wVGLQ%y6ZAC*pTI3Tp?fY9~>Lc8M7 z`*y076gS)z(@fLSrGL^ruT8g^eCKsUuPdgVZ&t8TXJrP2?oVKTQf7XD@OS%necu70 zm7BNC={K%5G4##UY^a3L1v=F}V_CzQU>`U?(%u1+j*%GO8yL3ny-Bum{lMVhs~gnd zrXf*!0zM<-=T8MEtW#cTTT*rwtv8aq*B6M8&kqU^5(RO#6?i3%B4*eA;q> z#2hA61zgVTVN6(-|8(E0*V*05dLcODbM{ENMcLdQX}+8gT>m+Hl)Cld;K9$?i&xlsKyV>#p+5aM6{JgVsnLBkgg=n9`#2G?_v~+ESb^;PrTzH=;qZBT z-}%3_hi3ZDeZ&57ZRz%AR7~e;>l*7Sc0G3{u*0F5CT0DXz>V#*>YAl2GY% zj!HCY710LQ>zeD1+?FgR;98iNm#T8zxyzOjFXMFsiAVVbUF!rQ`_t$>({h$xw>(i{OZ#?p$ht5r6RrIb01bYMv#E6$-8MC~wTm=Q6M zQAP$Ws40!ztHyfTuVd0Gmq~}ckx2(~W{Ac{df7#M zLBpSVIbSWPhhxhuGft2~h8qK-0C+RtgTb(m-Z`;c&7}DqsW#7_+GsyqF=iPXGbf*1%@LcrXx%0|(CK zETXYI$T1AF$-tU|5V!ttn___unSqOUve}2ua2!+C6=`7-a!Edr1s)mBrs3o*_>1GK zEKHg;@n)sbTt5DZ#~(iymdUwKp)R_?tcrS=apl4}rUeqDDla9Gm~H&GFhYb1xh!8G z60=5WS~#C@L zV^k{40*&!Fab?ThC=XC~=n*YuxmPQ_sb%+snd($Cn37tgPbgKM#3Ktq1EZ4*=pwcW zyd3mk0oxoe74pWMWD(wIa61iD**g@pam22`Fr`13xeOvi`&4kks2C9SOWWG2Vy3#N`?nv`9(rtQj1ydv^QgTAYuunNOm00HCyYN}Xo)SAn9XI1!@?w2*85pA#*np} zUWjtY!Fj>Y0^=(P+?jitkhiofoB>0`%y1WhIZGd8quGhr7@en&pm`=LBfBYdJJk&# zl3AD*UA3jwZm*!ThON_G!?QZUJ%c;w7@cZ50N^ps7{_>-L|Jj6GCul*J57Srh(n8s zCE3OHAgLSPI{jqDkeQV4HnePwN9o1@9eIR zWIBckaHMx!b7`x7W*)@f`2&FL&v)PYgu%CjQC`}xHqWx zoY{k6{tZxhZSgfzc3@bxjXa7xegEk>^K7uigZ~ea z0w&QmI7Fc{;RZDO?Q_rWbgsnm-}A?b>tWZJJdB~fE$T5Mv2n~=2NY3+E6k8ZgYpU z3#=sI+%TF{u=xwl@VZpAJuT0o*A$%cylM>gf5Evh$;BwT7`*Tm?7b=0F%-_oL(WXX zO>@HVSf|q89sT;whJDoc1G^4`_~%_G!~1r<{{N4@i?1iMo%XkXM~DIzWx~pTRn9dV zD?7MA`ld?y8up3fedjt-HpD3~a+Ywi1O;J*JDYuGXQIueR+WurC`2}z{WVfb_L&jY z%nPl3NI(~2(h7j$##8#a?l{(N|O<1Kr+Uz9)`A}>uR!joCLY} z)lCP>Ry#*dm*pH})`Bk*-<5N5PGtttZbB%6Q=68ML+tKliHCwL;CNsg7qm|VyH`81 z+%Xy@07uIdebSgoWN|1E6bh8VgAvi0HnC)HkH-t139|2#u)(wU$DVf~&mJjsdm5;H zq0TM_yAO3n7T~fXR&)HNgVs)GcuBNAL_$^TgYh_CGFspyS>T9nVOGv&;XFIm1Z^oM znqXR1A)*8ID_6v@Pz+;UfS{|DT2-tpVqI4ToJ}Fa)_awjYHPH%_+#QPRrBM)31f*!>1I9Rd z_^KjP7AuIRfpFo_BMqTr8UG(qwcr|PQ_>1E&k8sc??U%nA03(T!^(};B~DtVa#0Y$ z5a!=HP(f@p3L;774s=V*Iy6m7j7*vlwM{gL{%cgT@Tsm@W!Yi8!SLD9=KESaOlQq{ zQfiS3G7%08%YXz-E`B~#8tR&bD6Me0a3KH@d}9Q-o&>rq`UI<;Za+k9i~9TOs8h2EGX)t68%?tP68OE%44JcxT-N zC3kdl`hit$iJa0Dm5F{VE20Ei2lUrawf`6_HVX;*!EPF#tYHZ3p-KtqJ5n4qpXH_K zSq%C>{MV_u9gKjzGA4@49~g#Ug2{V9nAw;aXvVo5&q?Z*+^N~1fv7>YFckrY_Mt~j{M@9wc+$NP3t zT)5?=#KFB3bj6Yq4wO5QMBIx+;!nOlC;t6*Uvy(G7v_eSs|$5F=HqlO0u>YvwSBWm zc_kbk6sQpbM*4~13%gqMrwUIj&gBkYxs$DiHj5o%Lj?8b=v24@CNe3@A+ZHeFr+>4 z=Ld#Y9^p)~b_d@+!f7dq$BAHjhMJfpnF7o}a*v2mJ9;E#ilSD^_X0wK&4N+jHx>^r zm@<%etY$-@3Nq+@?7Dig3Q$?SBj$gw1M0pfF%qlUn9bSgD2jdqrj&Zxr#CLxP zYr~z)z>v$+*q{LOl-UIbo`z=`&V^HSe=k*8xWIv`j0*V&qm~&FM^RZmCc!2ipi-p; z<4bBY_l8wz*VmjK9!RjLhFcubX$Eo%!Fd8{T3_l=S&A`VOu z9I^5qt<{GIs4L>Dt=nHMjw-zZ#oJZiv*Vp|Q_T$!0TNOlraOx_^evz3%&~J47e)k6 z&U1d?t(2vRQCY<2W=$Nz5S8%v{aVt3O(!{13;20Z+4P=aiToqM+b20!;>YyS`Ob|z zGABD<ydQYYZ?yQR*}JU(%uc?`P9Jo?O| z@01TYcUX%~%;WQykShGW-(x(HCMw6gdk+e4JSe>FpzzLv!n;kl@9b`;Hp%|Qe_ZP1 zlh%gdqgOb;A0tWw0^F9R8N6yBIuVaH!?24i^#}|L<1}%OEZ1)MsI$l(f@hbpSvs=t zv-lWxyzrBgdgQAAFUBuZE2{-Lhr!gI{c-0hB~zOTBZBc)I}duh%=8JdrV-zJ;*goz z?cvn6T;oiP7p0pd>gm3-uW>S_{#>>GOH0nP#5K}v*543^O#Lgv z`VaYpGrpa}ZxFqp4m0scxGLWFtTYz1+~_>H^b_VQ9dfZW0)annMJ|3$_H0F8>=#!S zMf;$*Tg1?MP+YG$vVLIMyMno8=VldH!L@6g=HTyT*xBjf`|LHmPnY+3>WfkOOKY4@ zS<57@DOmAEXMSnZo)~rRk^xNgcfBYB0Jwe>vFG53eFsOZlsftsS#xm2hByK^G>Soq zkGP12+!UPjB}bg~F8-49T(bDbjEzP7YDc!hxV|qHeEvG8wO~H$grbjQP(Ijroindw z`P(ZVO}HeZSo6Eaxj@(&XB$;FGQlnmhcb}gMsd%o>DB*a)$lAkbvi9)1Je?Gg&@Xl zSTTWTP)*Fr!p5K`8arAXZrm8k2(gb*LNmgT+m8LXGvGVUyCjMzkZ;y;CJg2YDHWZV$Ee48`b`b406owlVR`(uhUHU^fw0u14@ z#gGTSM?%O5-yzx@+QVipe^sT{F=Q#tL zw%qyE;cdV{9C`j zO^jRqD^av<-h4TvE`%kjAymn(ks}JHQ=ewwM6lf?NrEAS%zfi$BNY9;|I4_4$K?}$ zH3Px4OZRBTo`Mq&eyLl2K>E`Oe}!lCXX)IM z4lg(3Li6vnY3stPRXI=71BGV$70&lBY%Wivc%GZT!osBsFO<@Zc@@%_`OFeE2d1pBHDHp5X}bD7CVQ0KvcAo{Q=ME7*|>e!s*S5&U8!QsWIZqyG0tTzhk{jKtvZ zyPctIdz^f?bFrLJK%L@@Lh$h2Y=+zvymGg*k0T?${tnJX9dc44%MhbWom1gvFiYv; zKl(np$yx7M{}~Lp*D0m8eEOk8@l$XnD`4X1)y6UaHq=OA|J6h%6J8hl#u1CYx)7K$Gk=4EIj!;+@Qalo}z0glf z<70%&3=2|rLQa`wqvE?E!t{*8&b|KoZ*=a(e85L`@H_y2(VsGr)AbLC->3NBOF>Px zpEF-`DT?Y*4~%7Z|NV^Shi7xf+%i_CG{Kn-9NGeQ&kfpehp)&%TVq4tvim_tN)oo$$iP zo4)lY(|6Nx6?a@z&Zy__>T^!|ETf4#c0kYZl;}l`?VtXxj>``AKks-cC(DtTf~4oU zgJ*zC{d!HOt+>oblf2SXcL{TH3i2xEvbQudW4Ahowekf@Zk?E*mWvGqeqQ6hm0-s_xO%JeY0gq5?{hA) zRs&grG?4dR0_=#H z0@LHuKVgCJ^bZB98`@C;cIHii*%UZc>TU{4$nHEGJn26#eKLoaq7neolsGwRq7W9S z?9IO-1#B4v1%hwg@650HgGA*dDnnH60cUaR8*#P%e_Q@1jAp*9gT$2&u)Mh?xcvd= zkfFkP^Gpj;9x)QxzAPd{T!az9YY#Xd)G%t=%n!Vk38tCYIsAf zTdQl=M=T|gD8A8GyUnSwtT*}wKIDAT9=q>(R2?yURhI7(AzBbEI{up7QL+iGG@-s- z+nvGos7O8B@N$&j7ulEeH%V9&Inv3Q4YN;#&KIhH6xo|JTS@W2R$x1VcuiekGke#3Cql8gaDa3# zMcK^PBx1ITMnJ!frzWzSH5BK9M5@PcX4}ognd8Msa}n`~f4&JIMpt&oL-B9E+$dhB zBK}i||F}9+0P&xifl!6hXTv3_P?nM4&StM>lDf$WWtn#3K_!-gjn-x_e|&;}3kj&r zULkmHv@_h#bBn{H)7N}a*GS^DHO2o`UVair4V@r@HEAKtAOU~P7d_lB#q2?})O9|7 z2Q!Xx#)oD`A=oXpPWanwItkPO<2F0EV7xQ1lrq<>q$G==1cZ&I|28uh9wmw=^s5%d zS91g<4P?x{uz5Q{bLr@PMN*Dc++@C>6i&FALlm^Ue8-zMlMIizRm$Tl^5&I;bE!Ic z7inO!d&ywVcov6m2A>%3Od2XSQsd@lJOND^!KC8VaB;(#@}CPDcR1s$Rl)2XPIu(2 z4Gi$!eqe_)cPzf%VzrQ%&>RD%@nq!uLfIIE)tI)k1 ze_Q!N#06rjVv*Q=Z^RwK zxei@tb|q2o$odQhJ9)mEM8V^^F$pP~D$uKV?oD#@Ep8~H^)XS1WpQsbpO`mT?G~8_ z8e&&7O9YD(x5%nUs2eV9_UJh4jJKMX6GI+n6rpi#%qBOoNsqmxhm0%;*pbFm1Eu0+ zCvFYm!YZzEE2Q^kj{pw>*&&jf^!0~s3~Zl~e`ZZ3(Tn>x$s$~v^6$i`(RDD-@w#-wI5O$As> zPdg$S#p-wzQFz-y;a(H&yYD&YztoZJ-B1Bvg7SC87;Qp(;!up=+a==ZV95(kx3#eE z$rqgKEbAA+nJ+qXd8~ht;4gy1_rO_h2^Q^f+IcMJF;Rp!0NW0vWlqfryjW_7Njob^ z%uSpWJV9YEDvVT&3d76d98$R06Kwr}Ga?7qV>YOKV^Dm_8Im-&Sp|z;3cJB0?#50l zu4$KbKuuV)D@nVnNx>5oo)dLrUnR$zk|XxPly+CBxL$l>`X`Acz1S0WZi^k%zvK+A z+hUuGy240=ibN)nGsQQ96)&4E-@=1S01V`IW1|frw?{A}(s!N;Grrm{<4V@7VaDsk zG#qxX&-<;jDxKQ%GLyF}*z}q+TXFT_+pjq-atUwCpPa*|trg(K|Aj`Q*wC%w>m9{E zsdBrqTqTbVEEi)c?l!pgPtKtKQ`wimM^P;Q&&+O3j?J;*zF88;0Rcid4d*}*P!Ld_ z55)siL=F$|LO>TZawsB5d<7bOA%JoTQPJ=`10v#$3J8c81O){pYVZQ&|E->x&1BKv z`@iJV+g*K6AJtXW)!m8C9gM&Xa@p^?lm5=*o&$nwfvY=;Ay<+gwJAR_RaVWA2E78g zD^o;Pnpz_uCKt&pd}32IAqL!QazoPA=b8_ZYj+{SQx`Uh*ei0sWM!ql{s_ydKuAqU zxtdWe&&u0fs4cmKQcA=%4?g$HQgLsy4;8JxkkW5u;n<$s=KFl5ziovR5I7nRPl z^0Z{E*T8OfS;>>zYRad!@`?5M$RTm;TY7x(Vic zpk_w-!O+P?Pf- zNX{&H8i)`2RSR(@(OAxorB1-cV`Awif3s;JW0LGPS2^}%`@fK8|U&)P9X)1Oo7N^p9Kxb3wvE20#Sz|Z0 z;IJ%Qid4r0inR&_R)YrI35e7fATyjq0GuwUnF51VG@xOmaX1+C-M@VB)x$^B|UH;aH5EwbCb!y9;fe8oOCBU+|Mr zf?SgklKGHKK1hU=3Bpq#Td)#IU?p(S zNJ&8Ds;!IC5GdLdt}aV&v?IoFi;&pJG%* z?gat)q!c3M>W_GEMCd`c2q_yy#7hy9w+JBPC07x@G?92qMRmGSrKH0Z*GcyJ^-4TrSw zu(bnT7}fuj_zG070+UbHHSNbU4*5+n7O(^vj*iuM=-$ zh^x&n2v^OM6Yn>M+CpvWc{_MyONh=*lh6Uo{Jm#no|UN!4|2au6aUc@uAsKnS-#4TG> ziwWg+ScDzJ?$fK1i_i&(yal7%W|#r0?6E~FBVHu&!rPRiTGm-G3*X`sWSZOQQ8qf= z&Pv?~kMcIQ1|Np4NU;24JJ>xga+FU_Vt$hPeBYx@ZP{vPvmmfoA6K za;ayGYQI;DxLrr;W-dYdKy_zQN7|wrtCp6v?@YV=#$264-KjB$@}QUI0^_{E1}|_b zkH&(cAJv@(fSq2^oicqOPs@^eQC+!*t7UJR2J7Fw>9+8~@10@xqYrml8@Tre8y*-+ zrtE6_aWkhOcmTWg1W6tm*oUb}FhGVOQpZA61YC)NQHGdxg%=cp6Vq!#aG}!Q0VvwK z2vyQaSa*pcCU#&$VQ5~7CaPK{IGD*!ZYI2Y=$t>I`}A0+ax$` z17*_c#<>esm-dcB6Jq7T2Y}m!fB}|uz|%Z#1BlF~x((oT10u~K!UE|y#hip7L6hYN zLN-}8WJU(KK_41=iG%Q_;IhRU2+zfT4#IQspM&sR{O2G%7yl(G49CE5$-J)*9qaWM zSPsu1+gtqq9n0}Kig28kd(YnO6y|9+H?z)Ah+Cu8p%9-U8hTNcj) zsi;O0R?>-9yjAUiHFE;+@xD2{G@sfw+B*k?FvdCB+KQLcbvRls0ziaUV$R_ zCnwJ8Pa~wbaDQr$v2P9t?Y0_rchzAuUyTC&5r|*!PZ{<0sS?p`AWW)bmYwcTIl8$^ zNwT@a0Gd;OTeMVrlagxlM+0b&QCw8^%s`s0$9?JZ`NR5j{a)_9g6g)d*1Q%v@c%xx zpm8|UKTuazDu0*hPrI<@0RrKSs}B(1qA{ET!0-V{5;;KN>f+jS#y&rQ!y|Hj0MtD8 ze#ObNv-}MpNi)Hf3N?n)2nwP2CWeryQoM((N#s_T(GD?$*)3C2sfGj73Iid7R#bXX z2NA|7YD`dwv#`sw6H|Sd)CkrK1df@$*y+BozpV2`bdApk@nBi@MKsQ*ui`(4(WO0Q z6l7i2_6K6Wi8L?!}Mhh?dECj#gXncU7q;75kj@xL10N&my}e(6%`?q%`=mr?fE z zCGK7P^cdu_hWCu2zIqdGcpYVho8T-2kaunDhucUu84z?FG3Xr+IeG-;8PNTT{KO73 zJmkiWh=*i=OW2(Fj)qWm23VA(euPcN7WvI<`noCBGdwH}bQN2o@;=Tam*5ZoN_0T627U7vyD-!m|T4TrnI-7X)& zWCw0KaQCMw-@yJOBrI-Z*)hlH4;jmsf`ZY(5O|bNjQmHwf3f)2?WMg+u5T>clvv|^ z;QEx)CAN2WP{kq7FaR^@KW!)f^XNIkH;qbngTKm|h(NP>wGskiXA&@4%f z@#t2$+EH%3HBDs#{MhfW6{K{qOR2FXk~=h`uyp|Cf)W$}yTqv;6rw;dlq|>T)-Q7X zUB^J73FsLQgF&H5UDGx*B(=oZ?#3XbT!khO4)){i*!crOiRc>GdOf=RzPsUNsiZ+m z@IrJ0ergWbPD!0KFvd0!=Lj~7}lj`g01EKClvm%)R zp`rUF7f_|ChPsvFyc{GPjp5D0tdT{?lqS$6o0z|ys8m5IGAIo3gL|3FpQ-%CUjkFl>tGDye@f1#XX0Teq zoXP3TV6}`n>lRmyr{eP$s@Tx~LWg}afy(3wxQ8dwwJ~#19ShIc{Wy_^y6?f24VXmN z`pk0SDE4#GWXcKd6t{4BR)yj-_>#%grfwzN*0~4EWrpKU=_)(H0wkzY@Pm_}!Z39E z{lFaVaGdJ$gSXKG<}L|#@;>55wW`GR^k3#(a8f0ZykiT+jmmK8%@iv8W-^hi3n%WN z%drUbyOYvW?A=4bMC$lLZ~;%dlS1LQx1qkZt(9t7z{k91;GIr=1~4y9{>hCea?o5I zmXEPC9F~WECc=bElH*XqX~+Y#kQr;KSPvjvkwZ;!P-#YNjMt%!ao?2evPw+gVlXjr z>M|d!6C%NdCRQ!EI})33V(_X*5DN(8@9~&`Du0kmui$WH3?e5$@{*{INXw-zI5v-e z$Z)VEB}uFj_2Q^{;1>6S3OOZjlw!Lv7OF5nNkQOmpnA|(5WGkki8yF1wrn3J4-l$d zq%3k1hwL(9tfeSVkh3Dl*|`x-taLGeD3Ec5z)E;0J3-}kD8GFdsNF01x4Y;W$N-D) zrXj7i9d;uMv$Im5l|tsPm!|J(Z4!^&Tx%+wG^@UXyOd+@p~ZEa;h)@BFGx*(2Hy&4 z4R5%I=D2#r)uNT8P;kWm@8X1#Uofzels!8gH7#2^ow5zL+ynPgN1O}_-bZI@c&*|K z!!$gl5}ZAlLRo?2b`QMbxH?Sf?iAtOVLGR;`Pcn4fSb&sJbf+ya~9>x=8-T2p%|Ob zbp~9)H|nd)WXg%*l?2?N*JhgBpIZl#fo(Ds@P zZfv?kJgQ!sMX4<`OhX8RP)Z2M)Y=&-#Rw&;-2gBksqzpx3XLt4`5|%!QRKnJtLV`Q z?E%_V5G0Ig0d5-KgNcBn;DYf{I5MtqGv>XsX+Dt5*f}&BNM`#SN>3JP{0D~N0|rUM zt6$uvoZ7n-?dWo988(1foeo zKlqf67)!7bZveA;v{($dgAfJIC5c{aYq5w%@@x|Yfmxhs%GnaMN+=&$$lZ+c`x^xZ zi@MBgs03i1d(wJhKm!Jd2Ysn{3{kT?=I2G_)F|-)mQ%xds6gYBbEt39BYG;3-UJJ{ zYhv(L?md^9V*`mt!>(Vhh8%jxH^EAh2WbsNNzmDvixIfhj#T$ABk{L$X@q{Vtk*nz znO^1Crx@^C&_YazId=8sTp653l~ zEfCOhjeX}wXhE*NLQ06QGX}yLoT0~a> zTDb_@FDJQj5eoK|+*80cjPFSgnX-pF-T}?nCMUF%u-0p)dQmaiZoSeYGk*kZ) zY4*iW!YavG8^gZ7M)}6#ls!C=1GGuvMCHU#$LHX&68k4OLT6_#Di#=kC^0N&;$kmO z4S*&`c`Z;LNIZzdR81m;pkOLL7R{)uaaaM&4 z#IUTN$i)sp8>xAQMg=_~OSWR|Wz!2tyUBzg<{V};u0V(dI{_;Qq!ya7WF;c;R$Mq(RKL~)YD;GRW;nA7xFVtP%ChJ zA3Q=rw^>D$+kTgWWuP?ctb+1s);{qDHHKx;*2KgK49!}$l)oyXk?tq;Sw{Ej z1Nrr3l-<}yZuSWSv;i7NmK23Ln4Y4XjU zWmmfmg0G1fhoaq9^>4D=(z+L8%C6yg#R^rtU5p8GnvWG%w{Y9#^k|#EtVd}Ch1vBg zw(E6tIrXnwjMGH+6&G^Pg}D5|Ic$~CO#?8G>;Zwp%Fn5Z8Zxqcv2q!t%uSWbmLyf4 zkrPf>d~DH@p|E`5OC9T)aW9|wH(lU<#pkO)=d9s+PgD0asVk@o&nU9*8IaqmENo(8 z6mkADG!gNw=V(;GrHim>$MsgzZR%JD8cLpOg@TJY!J2E@=TCyS@$%KwzV-nFI|jk^ zsJpf8YYdIv!DspKYWgPQV;ze%@D{FBzzv)LV>CrRMG~kVd(NW^P!cy;4+O7T zYH%Z5@!}V3C&zVd(UBr3rQr^&$ADYIJvY*&jmFw{-`D~uC89z!6FT3g`O#-1(V8#S z|27_``E&fo2FSdeET{Z!RQ_c{b(XXL&n%a7@fwJV7d}n3GWMO{Q}BxJseSzCMr`0Y zpMTlkfBqk(suCDagPZmK+i0?i@w(^KXGzpX_x){D{>Sr?9GflyulP3&4^H#h7wA5{ zKM&cCU5b^w?gg6ZLaX?)O*FRtLZ$RA7w?8}yAp=3EMa(e=c7$D0k7BOy+}Q-nhT7K z;cibFb*rN6NU4KOufBsJo@pI{WuJ$1JgTBS3gqn`_414pqtN!?PhO-vcbwOHiP{8* z%C;Y7h$5#kLU_nal!jqH`6ZfEv$}FCw^QA!muRw^%fvEz)~MLZjW*M@Vo%>piyAq) zBk`EJ)VRAwCP!-d(an^~1OG!KC34As=mAnB2)BNjrgn)?Bk`e>uD)3JEvo8# z-gonfR9i=#_r9#^eCB27hpgctuh7-NoasqBPv;IhsS(e48<9$Qy8A13usX_yUC=>(l0V)>X^qr9G+0Pobz4!B7fmmqiKN}sAQm)6V217s zb51#Si2A}SGR1adr3c|8yf;6%n+C|U0&++~WE8j*f8TCudNEXp)Vp&e?~np(BJtC9eJ96{EYflLs##m&iY|bmLYz3FXcvgb=rr$GrPTT1*hnm zG~|!+$n*zA&6)R6UX!yn*$k@K5grMnplIYX^Z5IH=shGHwh!(88xPx04UpaV{dCxE z%S%6}aC4CaUyY`Q+yg2mU==~^fU1;SfKwji%FnU8R)7^1q`ORo z`bOjvH9FNLd;%kR`N3+LN~F|T?0+@rHXm&7u7WPk6+M*Zk>A` zqdSaRIie+YQww5eadz>`-WOVA&K+9p?u>>A0N6vb@j8~Ac{S5(D5g~p-gweWxJi0 zi3wu&4R<(7t)fLhHyovPcr;H9Cv0>fvaCDJkI@ng=FK=vkXFr=#74Lhp@y*kW-}ud zr}^wLx&?hQ>07KgyLmG#Q%2h%ohuJwS~vZUrW%`n;g7zf`ueB5@jIv(@3O(BpF+mg zOV^!Isbf^?R$Z09b*7Yk$FZTb+u+ZSQ(D4Kh@&xP6DC+YxAEo2AsTsrXB`J$Hit`& z(^Uz3mw>QW;d_^G{0a1O3D3r1p!OwTGQcl`-$!^*AfcSHa99;Tp_#bR07s~VtV1~5 zyypbn*HNNn`4;LBsX`tp;s!97nl}@<(b)&dJfaPonVN`%bH1m=X)7eJN{It6+G?HwP zm;6Mj)djr&6IO;u4Nm+-na~Cd1mYoU2x&n5p`WQKBo=`{EntvyK;jozCtv}G5C2)A z(n&wllJE&ADLF?#!M;U+Ky!SvP{L;-UVy72K%l8qOIS*gzyQsy{u4EYjP`7w-~&FD z+G$nc(XB#kj$3KgYS|-~F<%#{$Ai2J@dH& z3R++{zgC5s&*8mQlqUEx&`^SlP9XK;_)+yyNxA+h9AkUN$k)`LZ2fZoQ}ng*)I*&8 z3-r2o^QvEHE3BJ-rKj{29P=Bk1N6dgWDV{WNK)l0zS{_1R1-|D38vKq+tviTaG&2n zLoef{ztfOp(Kf;NU5FE91QilK^*d(gat{3gqQ8*G{Xs+J2+ig{s32hpqAIb491j_I zng+xcsux-diQ{`#pQ5^C3Z!uUujWOH*4Boygj1*Q>%%jgh z2zEO^e})$6MVx<@@)N6`M%(JZ4Z$Kkfqy+ieFLi>hd-Jg0&lQGBvG>o5zPy zdAwdTj}6E9#Xl)ym|YLaquT^)5e7vhkDLkCJiw89pkBY%%mX$1wx-X%sD)93aVReTt^SE9Cx3i7WL*}g@D{HRdK067C=Tq%56AD4a1v!ffTfg zkQ4a90uh8yu%^#w#{;4(;eGVrXi1}2r;@o*TiUh)y>2tzB4U+f-dIxug57W(p!_tL z6d0rHk6&nZ>vB8z5bM7V7*z?Y_-MIg3bmKVH<@XJm0(VR(0`B`6wH)?D&aWx2N0`^ z)e-ndX1R)R07RK5!pw&ep$#J7Vcj&{?Aw3uvnap^Y%LQUiE3-siiZ$qscu!mb~3V- zNya~;MCAdLc20=B-d6brbjfaBtDDU+&ED0`Ca{j`W^Rv#qQ*_qR1jxv6}GprVteC| z3EDJH!&J4t5dV)f6_9`c=@&itG5^Cb^Nd}^yvQ)?$DR@FiJZ99FdHPRR@JQsAv(p= zwo~~t!@MY5a0*e3;8!vMmWbO{0kr1A9efE0+{((+YFq2wpxV$HvV*KT7##MfMSr6m ztcM}|jEq{K;0FLlM(t^&asZ!IN46#~*1K@A&y25IF1bT3z^XIBDz$SD7V{fEvwky3 z3Pk%zFeF5OvJYNI1eA(6e)|+x`OLBQek?@Nqy1)cWA!PX;e;(gUW-DXgJYGp=PVReB9g;RN+z(eviu>VCb+I}{ z>S~(MZ>c=RS*Dqt?lhJm%aEJq7nP^|i4`gtPdClHX6KRqIECaP{y#O%T##yJj5*yX zKFSZom>pt25saCDUVSyjya9_!i&(R;_CfJ9FlZ7Xg~d3LUf&UG4$i5<7{_G9CJ!vU z(2ZR~nd0EA(2iE&(MYe@pLmq>4?Lt^@|OVPhlzo4 zGyE8CselWA2zLt)ycKS#m>d2Y+}`*v1D1Ak!#BY#?PkMa`^&05G_3%E-Y-}H*jwO3 zfV~CGi3*2R@0GmJe(PlMSO_l&=7w58DGaR$=Y(D7c3Q9t8IlQLwatTfj)T zb3L%cl@D^muY|j!1BazX;xXLu5O@?!x(A*Pw>QHXfU(R+GQ1ydH$zEZ1T(-BF0gcg zTY}_%hPxXY$*GRLw=>?Yb?rcgjG_AP@4ghW>a5rcAy`k4kFafZ%h?~JU zz#TmBDK5)2(;~As7s0(f!2&n-z*-di31AsgZuVbB!4Ct@@Py9;+{**c1dJgXj$|Mn zojnm|0KUuvF9s})>=r1n^soy{d~f;$!14i0fs8*Tv> zgzLr^Si+5XyTA|NZsLJIhg-(H8~z2{;gHAU2;6Nv@HcRWJn$j71#h_-9E4lOqzmtb zTkw_({|UE@aToUMH7k&$qy6P8-3alB;4L5qaL5yG0`?Yo24SsW&Y;QV%VmXX=z;S9 zgK~<5cLvnp{w{n|6kGzhfhT-TRQSXEMq9J~ z1#e;fmu$*lWJJ5LWZ%*QpMYCx?uN&3`*vnV-JK~x?FM-7M{DVtwQ?&cj5$29o!Pct zrxrmqC#2Ek4@gKG)`B;*Gy57V;`pa_<`iRf65rU~%;nAP&ArV+X+d0Px*a*34ebXZ zAo+fdeC4b6$MS|8^DYB+dIz(0LpS^+!lgG|_#UZNL%zQQ3PIfM9n8F1Zdd{0${c}z ztz5H<)fpq$l!xwFUjE4&ae+;@Qq@JGI_1oH)Qe;Vz18R z+%B+BWb$2N=VbBAV&9R)e~P_1i~DwkeI|=%i`}O+zbke@YfkD0dq-=&RP2N{{4cSG zwc-6@7q;PM7r@@zhOZO5Nj5(%_ULT>PVAN0+%C`Tq7UM+d1eNzM_{E40{eC~8X+ap z7;oj7y;EK9OW`eT?ZU624Q00NOyTtI=JVk$*`>lz65TWcy729CKcL<0;j?*h7oGRg*U-1 zqtk`gz+Kk^%Lix4QsPF~54TVdUf7$#E`)huL9Ei`Zu~=Vd*R-I<%8UC`5-U+EYgR2 zdLjrlB2Dflcn+{&ZWlfV*c<YG3 z-xUQAUAPZmS;Jje;(Nno zX^{2V4L^za-2k7$T0Rc;Q5e!RA3XTGK4u4_{3tiI%={jiAPz1@;2xKht~j+J0uoaPPRFc0B@0{H#!%Fah+srOJ27ZY!`_WCe!0w4crCayM^9uvam-R5;GSFPC z-%(}^f&@z+!S`JZ>Cy;({bKX&@ax@!+IE;dFh^m2fuSBjtv*a9OedJ0FhgKQ!dwTl z32|r(o8>yaZDY^Hxu8GsJABU(FW{G20s}XY%wR=2(3w9~xq&87q^^%uCE9U&>ri ztHNj9hiR1#doH&cZhmGw6Da$0xY;)*<$g3^EBMP0ld*k255K}J%bGeks7-^J0V5qR eM*2fK+3y%rvp!#1*5OL?D%~h8;zviBL;oK)-rIcu delta 69071 zcmb@v37lPJegA*YbM7{C@141`&t$SZ=VaeQAS+99AOu1Tn+A}5i^>cxAt_o7+*q-q zfJG+sqJmB)s92$*flS(<1{*d01c@32E7qvkQrBQZAymdfa#O_j<=ww==0P>=h_PbRcx z2yLpQuJms5Zn`pE*|cfXK!Z17Y_l>EAkDp{t~8aK^{?XHw0*w!QR8(DykI&8KI_}c z{8ZZiwf{>$Z%nu6`TTfZPh*Skr%J~7e%kZ0B>BGaJ%0?T#$>#V9=*nRjb5H2CY@;= z>zTao=Q2t*<5G<#pEcfiGa;R|f!9EEK4(l@F5i&zR6gscGRDt%P3PxKlSyYA3c1iX z6VqOkp)oH{QbUW$q;s@f_PjF98&6eJDAVMZOqw5`7Wv7OE9d9@Nq&Kj6($etOZRxK z1HVj9>YUP?PUrk$E|;Tkxm+sO-k$QN4ou3NFmqZ?RiyH5xso@1=IlZ~ONL@DlkVu~ z>}1@2I-N<+&7`yQjvIJi=2q{+1I6rk?>hrCvw`=Q18cH9*WKrJU-!0a*Ijef6>qw3 z;}ut5x9(lCMdy=~)lS6!c)*P0}4c>DD?uDfF64Od;WF14<4!@70XCV`XNT+$6! zT=kB3T>s}(as9@1H>3`F4rG+AB6}ylUgR)IWQTN#xts zy)(6CU}65ao_kDt61RaG)#?>*efzcRQhz@1-h8k3xq7@B@311{aX#nO+TO5n-8-(h`l^lBye0L!f!>CBv)|{pUU9`cuLd}8UiY>u z1hE@7&QJZ!936Xn>emDBY6vTT<2MYZy;S5C2F<6uTm0{t2fPQ(-}zhopO`(qeqS;B z%s+Y`^uFwW!@te{g!!fSy!X%M8Sh?ipZ_)gL;ii{U%h|ve&}s@+I!0j=Ig}$-h0UV zqW1;wr~X#+wErFdpS<^aH(&KFf3N?f|A^V`f64!e_i^v-}?vrXT1y7T)2x8kNJONzT$t`-{-CPsdvEpxv84}Lxqoe zzw!RyUGk*=l>cr2JLa48@e%LW-b>ys^!Ew>TmEVzwiHx|7G)4^GENH_XGcV|0DjQ1ttOIq zHeWY?Y2MF>f9F4HzF}@PA22`i5BiUpz2-yagXRT)N5y>8+-`0&&-%ajKkS#AN|f6Rx?9p<0>7yVu4TjoylSLP@FzxzAQ6Xq`S5%X{6N%NTZ3I8F#@n6iZykD8; z&A*yon&0qz-u#RCh54;1e8Iolf5`l`f55y<{(JnY|L6i_GWYwB`d{~oUcT{ruXO#= zzvHjGbZ4b-b9z($rirO&pVuE|CZ((?Z`%IT0_^Q`1HWiG?WDi(b}Tk2n~w5*c8WET z*SB*hoLsI)r~B;G^839V;q=?>bnExVnTvgPa{2uOS2wSondEF##OQ44Lew3l%b)YZ zbg45+N8RP^15Y=QCGGvNWkIo4(siJx_3Yzm z)ny{0mZtnuY`xZ64mIX74LsX2zs1F-tzR!UFrl?;O3pU)hx|LIlqwW$A$sf^yHnAk z=;5F2PL-cFsRFf{@}YsQw##}37578M8QWMs#7p#iN<-d$N83qW_jAwgP8G)$82a<6 z@^+upTRfNg)AlF4?m?xFRq8EXoceNH_G_Jg>vcc-BT{>n`f^hC_PFeA-TS=mdz3mx zskbL(pKxUdHjMdh_x^t*u}6tdBt@T#E55AvLa+NyrFJXzxg>ReochGrabEWmKP0tF zsr!@ExZ-Wzc(3bLC3cRkc$$VY&^_VZoqL`ktwU+k;=!Eck{+1w_0G*oYFE-pagxJ| zfm`fJosT_DQk#-M4gu%{SLz4$vz;&ffTUI>oe(F@a!Gdwx0fGLQcFPv%P#V`aAN1d z?-N;8{&Jjul1u6f-`@G4lA4t?DNbs4N&6-)=zM-ZNli*>kCVo_q@@$jUva;ZN=h0V zCl!?hGN+e=0{+re-lk7|@8->$%Zp4ZI`BPm7nQqsL>njGtE9rHHeR!$XTD4B#-yTw z<&$Qb7oQw>>!fXvP0QrX6F~4_hA^M!bkMG$FU+q@S8OI4^!m#84t##{$zu=CQ?R9# zbh(z{BU6{Xc805_>};slFc3`b^jxju>T%wV+iCjfYX3xztEOB#aPG8d(S3f=6d*bi zrENpxZwyWOEpACF^ z=5LxscD60@s3`Kl#-UlCBIV54@A3`}JU)9RaLdp6>m&Q-+osV#`O@4q1G%~LDY10! zTaEzFQCYvx&JaBHHOmnEWC#O)F>mFl1kdJOm&$hyJUegsl6cT@iN+|C3>u9~IOB<- zh4sH}n6L6V5#DI-6&L zn`qJBNCo%#g~$=67%xmCp-B4M8twTRm40zx_=Gde?ccd;LA30CMU1dD^&({($hlt$ zKM-u`wdfc3neEG)%-HU9TV*3<1 zPyZIXlr@Xq(W)|vDeGUii+*a+^roHv2{?PB=XjJi`vcj<^J`K2lc<%8XLx@#uzvA7 zyupE=End(6H!Qi7|6f=#A$#B;#U2>=!IBvO*;~5OyKmskr61@2b4w?7Fq^=5pR#Nd z*UCGIF!sdV{I9?MwC9(2gk&3_NsFpU61Q>9YNW zH!$Vo)g4C{hvWvPo*Xu03ej^vrPWUjJbdz^$@}%Hxr<)j!&t{f2mgtP&Ge}JS^hy| zX#`dAl=JA;<)>WTCB!(8fI)bfQO; z^z6WUS0zcee6F6fCr)ZyadvjUUY{Pge8t&?=iDm;D}Q!ic*U$UUwWWk;Kf>!EhK0e z-SW>Ad_eQCrT(^gGhLV%-T8UdDF`{--9Gd53q891zLmfDU+NzpKjR4f<40%w@AQvP zoN+AuX$A=~FA2_jGK^IL$ZPA9hSindSNgkIgyoOC#PGj30RKAG;Tg7e8-`neX`1 zOUS1ta*Vg0ZF;#eMQbLQC8-4~yguxL27)Fk>KfSh$Ii-l<)x-(SHXjRNClql3TA zK4h%6sW)uuSB-6Uf{ovr?Dcf*T-zBm`^tOnp#oJkN98`6YKzkSVH;X;0mwAjmOix? z^ebnREw3#36mU&9YrM@GgnH0)h7{>^fgZY^(nEXwK=9V2SLrn+^FQji^8HN+ul#a zzIw!{ktG8=!$IwO`&Qeu-S*x6qGv0v%z=8DEU*Y#Y`e|&2kp!Zk8ay($7s^a3YBe7 z*kp|W!DJhC18j+9Df19Z*Nn~X915if}Ir{UlF!; zL=R{CTxMnimVk7=7EL0L_-<06D77M${FHC^Y{xXSCmR?&cW0@F=VRI6cs(D^h8~sk z_DI%xOlqt{*yAIV+s*GL=z5u&g0*M-l6kDT1dHW(1L3|d8);rA7XC@N2~_}5o%6k zUX7i&T5xKkziq2}F<<*^1wK*+Xd-yq)z>H1Vf$?ploi+Gx(!W={K_I&2TnmUy<7l( zCFs*!361J-#X_f62i8}VrS=z`TjNbNur^ml;QEKd;Y#o@xGFpvt_r^jToWckua5w0 z=TTs-yartVkia#5`sF{Epbxlq{!icvSeuUoYhMgjS4{%e=8E9j>=gEL5wN~0sy0g&I*=X#wo#~RYM`qe z6HV=F!|s2HsmB6c<^KolMnHE=O#M@!bLdz}#b{l_$2#mVkB&e`agsz199$M$R z-}cqLB&#qS3Z22k^yl}0}Sl_*627?PM!H*3E7uI!=U~R3`vf!eT*Nl%6DLj_meW6`DWHZ59 zdx6ai4Te)QRC$s9uI&vkTJa|L^V?-Z*9f_>SY2l?uuawP{3Y)03^sL0`3r3ZZ*X`q zJy<(9L^aehRBw=~FCN-Tri=AESWo!M;Lui=ll&_Nzxt~`u0HdTCD+oYOGb)jww3?W z1eZwuZni5!PcYtW)qS+^3Q9Gh*ExpTVJ}?b-;=c+_JSq;oqCQ8nL5JE5`R#!m(=A4 zu?#dbv~7w1Kvr7R#$Y2@M{(DzU0;jhP#sB={3_wNr*va=JNFKA1~dte#Uo`pLw<2$znmyZU+j!@RZ z>+OX*g16cWb_8!D@thPt2VV)^jC#K);WRzMJ+s`=z9gemkY)C=9e2lSh+P2$t!$!m z#l!6kaD0k1IS~Kc>az{rl%-RoK0C5F32|a4kROYtFA>&ZqeNMwogwO&kd%xC_m=9K z@m_EB!bY!r_6(*-=&#Hk>kZDvX@X}qqc}M)j2NoIJVv5eGEh{0R1}K~ebq-By(yJZ zTKJ5b7H)G+3%7DguJ5&TaF)SgY$;L0cV@8_I%cPy$ObE=bnnduXT-{PFf~@kgA??8 z47EH4;0^IZ769?PM=9wUN7mK#KTD7A2R&M(Fvr?;3)*#25&DB)Hr5h~YkG+)qIjdJ zhMm;e2y%G6Lbal%N&Ga{6H}?6AM$FJ-cqzv1vhcgpj>_Y1|OdK8k)-gLuvbJ!nQJ^ zsWks-eft>FR=Un8y>KdH6z`Kh4UcU`YXJolwHp18H(3Knm{{IsWC-_He^T|(^n#v z8{A-_2C&;>a%U84L;w*rR)5mu^;G7+T9shp4#gy*E65f=HI_;| z#@`w}o?2x)%aN*LVDT!1Qp?fi{75o71fA-el#(g{y2gHnR8X`mZVBiV5ft!Fj14@I zym+(qW|(vcXf6rMdQXM@72AIas+b3c^k2rfi{!2-ryYy^LtVCCo~E(RU>yTB5!dfZ z#sxFZE)L@I(rr|bVg(0(`6%cP^D1({PvnTPAg%O?Jx(hf56&d7&>G-wro#^Nl;5k) zD|=%qt-|BJunj+#ErWtN`BvKdV7-zm&5NJC`)N5073{?=3;jMavrN=7t+GX;ZqgFM{u%y;8 zp)qZVlb{d^YnzpTSFOTmO?2nui9;(&!Q4%>Sz+-bD4-mxb*-q+EzafH28bmrq_?C4 zh7ok}ce0p7Z>(E%on9){-a*0o`^mf)qWN^LHJX1;Z!pQWo*Q|q!!w1Z;CcT}Bz+{~U7`J#(LHP`Ab7>~bGbL>W0uoI;FGcEB{ z*a_8(TfJ$;Hzv9WV(Oym3uC;F>b6#|v$EDMC}mj6sB0=Hk{WzRl#2Pwsj`1Idy61$ z#w_fohqr6r|+5wADt)Abu$Nb1n5<7C^tT}AK`L(ub=TP~x245J1w&jY5+_s+M z2)C#ng8{pOOf1^tWk<(A${%yyI9N@4jItJWkZyuC%x{zHtJDs<(qtz&_1?MYAlpen zJDqlHDF$9=LKb<@vRK{O=8f%Xk8rEPnd~HWrP)sE3&%u$za3NkO`A8VQjA+G(|CmA z7S~S!=ufg`8S7XW03^ZsnSKh5xDz2mlz`Mv%+`hkKJci;rPyn1hkHv=$Z>h_1-Zu-i$TmJS!;B48;foV}Nb>-FpM8eF3qTJ-D z8%`%^Mka~wA*&cGY`Y}53Ih^W9}{a&BH}TzPUqSNSst!rWf&ZRZk+{My4>MfEOMi~ zhJ6a!30mSZMrQ0Xgx6v;!^(!fV7WcfF12&{T?%e5>9eTYGcBq8JA*5v74!$IoPraa zg)1Gv&V%ps=(a7Dn#-mLYeaQejRNT~y0b5M9Z}96Jx(0DV0pBo#JA{=b{0~9AW;WT z>JQi8PZn969Sh3!;i--n1S2fq(i**2f(w~^=ZH3|AZ`R!Gm9`yfH}ArvOz;j*^76? zwui#8`yf`kx8&EKo`^^RvIn|!60Rp%tCok>wb?XQzz1q>-2;ehJ*y3}-bD(BZ*~F-p+X(w% zB*R)eQ>rQzw+GE5(HD%oj^^aFdVvGkt}t;RNM!E{uN;JrUu3VexCDZW>>GKzlD9Vs z3Ie_4MNsYS@rpCvXs;NO*)LQ6qjf9Zv73wHqzkb9Guy&Ts^|51-EX*r7O&VQZD}E} z|8<>$-e@mSJr}vmg}Ubgd6Qu{Xc-KvG3LB>q1L}LZgO@}=921{dc3KX#mr;uq5v)# zd!g&YrNPDaLaI2$@%gCXEYau~Y{a}BV_&y2EJG;pd&GJ;x(Ocs9vc@Sa8t%s3i?1vpoO_E8^AllR@qN%bRn;*97V6KBCx(I-9|_Dmnk1hCvlyK9o5M) z9@TobrCjw^q;TT9pWp(9f;#5pZX`BZhK&XggArW>;3)i*94-2T$DuLzld2SW>6@rf zDo)%rmS)tXynk8+!glDB3rFL>$f}FGm{T-$@E&#{s>DQFZ7`ry8Tu5CPsjf?0 z4pCBhG$qvwBo!9usbByUSi~M#u<{&UMxK$Wm zVVI;HHdV2bYy(Gt_5#V$F?8wK8gYkf<<+|qQd}Ekq@to4pl(ME&T|**$0Ks0t%jF1 zOhjUh2#gja2E9C%v}j+Bk%Qch@vtqy!?sv%$g~KHvE1O{>KKKGqvQsDCQ;XDq^nPn z*;EIXWzXsuFwNK{4$~U@k6^{G!gK2tfHx=s;N_A&G!?nEJ>;^i9I%aP(mStpM`|i1 zRLj)Miyx1mN9_MyI|P>2ETZMEN*Pe`z)`fs?8^4UKH79b_l*hnhi!?cy9Q526u^z0 zG};wjI_Q{w4eKueKQ9$q0xcvQE&*wU)(VRPE~-8}!5cHxvHlBS=|@Y2OC%L8ux}Kj zub$oGjjJF+jv){bJIW_4&ul9MZ*-i?Sy64a7ygVBz03C2JT5Y(uqQBj8f`D*ZmT60 z?8(j|Vn0?La=-}9r?z%ck^2sp<{E^`G7={8$A4#9%tcF_iw|k0rZ~j&PoXv9(|QYR z)9OGdAT?@cpO(+r?J*D8Tx$ulOI%t6F^_2o%A*t7;qfH;Cqr}M_{|6{Cw>Phl?Yns z470*9K{U`+EX**rW*W0CLThNb`aA0_sfZqvd}z$@xG{1di$xyZowbso_3qRagw05u zM7YMXXRU45jN1zm2R-Tdy`(hC?T}att}{*QOoiUmR&yaS=o2Y+SOs+?s#i&jbMuZg zLQK)2U_Fl$hKZ3P3sry^4Z=nzF1+;>pb}$!oo*aYjtMHYOT0^-{s5z<(RRz(;86g= zNW9r{P(v6AT#85P}~ajrc;epMHZ7 z8^TVSJb}K znAX@l8`k%>Ij>NQr1m(QeH;FTtR2U;97^YK3^m%Gpvzvv3|-CKU4^8bhSM$#ykY}} zq!1Rqrxhtyd$mvmPiZB&h2q-Q#;zzLy&J>InEtR!hO;BaaUcfPOz5#@+sAQ>98Jn` z(qze&gJ0ydfE(skmGZKpHV6GoWzH65FSR4e`GW^iiv3pID%)ov$~iWt_S;~QRHQ$) z%#e%SZT{5tgmeATrq$_6loCst!(v8%kb_qukele3bNE;7RP6Wu>e$Kb^qK?vgbZ6; zy=03yNx$PHcLhAffn(E+Efuxa{XGl1u<V-cEdw2n1DuBWH$j`6 zlL*B6Hn(JXnK`uVRQ_t>poGPMa)8kxOdIXAiXap%v$Y8U28O}mB}(JrgHPHX*}b~t zIc$5l)GGLNg5_<0`J$n%r3^U?2?rcSawj8Ck!Qs+{Jma8EhLWoOvwX2$jQG&_sUaMCBo^EN7VUzT+W9nE1kYP?R7w$t^pefQU@ zi-(_{=6%TXzBhdJ4DUk|yw42(W|cSNI&Zkzz1>^n{b=~y?cVQNz0Xvaf8ASnSnU1& z@azANx4M}kpZuZs(=pjcehI!kJN%a~dZpspd+vQ42q4%ay}<~)=zgW|^-_i1)6(NB z*-e@E7;PM8X{VDDb8y6krYIq|(B$hpJ7(Ezm=U|~tj-I}vc(J7SQ0Ei7=lTSECZvM zxj4NOF6+=J$TJAXg`p8sy)HDpl?F*aH;ohGq$$*vb$lIPeUJzB+k>>t zgPjc={{z9N;hO3`9d#fb;SWzc&*a>X)|G~9oJ369oDzb0&=F@`7$&~|{xA(3=_3Yj zm-eLehZQz|+OE~wvy5HkHXmc#(F_&CNYZvnAFa1R8EB|L%I4YQLWKg{+(Yd^aUTB{ z(m$LNWN$finMobG!Y(|-ir#`l|K^2@4h5+l?Ca#|?a;=<7Sy+TuGYh|(YZy{fr(~v zo_*KqKzDWbMAO;DDZ#i;l#i_qJ8!FJ&s0}WG7G&-^)Duw4GnW`CTL`0I;wx1WG42+BlYM{ z#?1&PSq`Nj&;`|{lTA=z*?=tke z(I4rPw;ox;w)^;Xyt=zSfFUwyxJv^JQk-oj5gtNvRWOk$4Sa{&X1O+(JI*`Sp$)jN zf~&%*;+~9y^|(C#$rd`8w~gv}xFDGno^DbOU1Cz|aY7JyfAh=}KJgaOMR%CD1!ol( zY`$MmfM{@lt7_h=mrXGe&G19T2$%C@1gF~ZbwF)A~$x!UFi_7uZHeJ8D9w+zLz zIk$W)m5favEM^ZeU# zSFp4?ZYtC@ZDlEAE!v@k{Zg&^u|HlcK3)TX9KkX#Y_;R!SwdCmV+3O)TE^I!Zb6^} zD0W4x`Mg{J7K{=i3MlUcwCJ-@V+IZhZH^(oPXh!iPRxsr8|kFHbJp`KqE;V_mIiPeaikc`(YMrEI4qD5xK}Bgf8gXhh zS_ix)I_RpZu(-*YPTNb&fqLZ~Qn2d~1zRVMHhkv--8}5FrZQ~~Tud_LxsJj@);Wd- z-Ke0qdZ--p!X&M`Z)T;I2dw*kv3ux1LdZ_2#U0G)bH z4^v`y-IxzOv9_Bc8h=PFs5Ys*37gt=g!hRM8`ZdlJ~-=YTX(ahu3aC@a@;oT<>4bJ zQ*@#i>xX?XH?Z%7M_(|N$Ia_mvNe5yRrWp@Tb+%Ie<#d+IxL7{vw7SsQ>B;L4#u-v zs6ALz2Lf^z2$5u#!pQyC|@IVQLiL#oim? zyZ%vpSBGHhd^h+1l1g zJb17kIwlWBDmv{3h&eM=SebRKpwYrNV|L)fEJm2H2!1j(AHjg@pglrWPUyrNtxL$r z6HjYW3$Ux{K^i>6dzt8Qi(>|N==UCK-=TtEy>*remP|ifZAd${0ja2p9M26?(lb&P zvXMZiKEc=?U;WuE6pC3$T2AFTs!fR{Bt64b%OpskCsBP zQu_KLHcL#lo|NI^Af$1QMB40eHbs=1CbWZ9?9^q9V~}cTx|Cy#W%Nf>_^GH{fC37v z-Z9&pQqdj*w~G%gq!F4=I5mofPy3FGH)^k5qht}YFiFp?bjC%v_7J@xT1#CB&9SY@ zJ6RcpHf_8Fa}91S+MnU(ft`tRSj|0@k+Pkto;SzzW)|QFTu@y<$IPiLIrN_%8uot( z_+s^m9>*hZ7f@nEmxK3cy-L7=-gGM|Q^b6mU6QNRM;urlFwjE^5>*!`qgU*eNooRB-SZYK8mfIAe z8|{|t($))|>CbZx{Sfao&}iKh>TI5~s_;kUGmqsh*ty*!oFV8?!H(!HgG21hhxC%( zEBm~fi22H(omUT)mk$lq5|#~qvO4=X(@~kXEhxc_rS{+!5PK*gs1`8Xykiu278&CL z*urHp2enTme$$Nop6wY-TeuJ$JPZVAUZezv!9}LyTxT~@n+sJ^_Any@rzlMK}WCgAapM6y9Irn39$>zcHC&|kexd?x|chu<7Sz0>uP{QzbRrs zs~YKifko{=7n@j9cA9GFk#Wik>d%xUNotrjfI7!YkHbUk>k08o-sy6WK4i|d)xP7+ zGV|EB>gSF(r<%dR>hF#>_nDpBsvkbVEG<5~_3@NsrWftpTHSwwnKXI#w#TIbtSP~( z<$riA4x=BH8f@wU(2 zjvDPuYQ>~Lw)i)4;ehsXx_ESq^bJziO!ThxQQSj6NQ;9hhATi^9E{tj!FMk*)79Xo z7MYXAs|`7FQjwaL|IFQZpg}N_>MM&(RH=1JM=i)^h-OD3x5!u?Mc9G2u=>sf#>6tihHT&ds#^b<0SX2<%?b#opQ?%*%YJUZ5|!CT>-@?<<`-GEea?` zDRiptY@oiQRvnrNzZz|dlk4dNzt>O z?S2LBNmBHzXM0eA!6Ze`dbURtcpyp9GtM@h4v8#bQf#G#Pf6DeaXm9pqf0>kA6jl^ znP;|DcP=;cbD~hlwEDtwGsZy6rO2G9=UI`N*K>4*I!7rMG7!aeJRX@DD@D7sk6`mr zuA&FxAQx#7#|vH-s4$9&JRB$gvRc#-5rdW1#8dmW1g{EiP+ zA3xQ!#nLHIc7(-&Q%!%f`VR|}k?A7`^jMyKaSLdf)I5&gGIuv-iH?dtYzvZSz0v6^09I&1_bBhOhf` z)6iA<{P^@Z7HM)fW#1Do(&XavH}X%slnFCPEeUYhVb7$)33kFzFfr0*STv9+DcI{S>FGGRwB(`I&JTbii2Di}}yw6vI-oq9Qg z!N*8;X95H%2LUsIUz4L~#i`lVlluKaC36lciu_O}+VXc;9zOs#5cio?3_PQ9&tyWg z8YeC^rGYDIKhN^Tx{T$p9K9|_ExjswN^0KoXi{?$exx48Yo(@dsWt6-mU=Hs{aUmF zP|v4>2~snt_I^eXH!q7?JqxmWK9C9O&e&8i0nEl}=B#Ku)v_g71=hyD3D(Daw`u+j z>>3ORM(GuO)XDiS8R_Z>MsIw1PE8 zWM$D!^)gev{;mF)3X`b*{1LV{OTJ@IU3H`N?v;>FrxZ}qB?g&_D1<*jaz=ue9*+{ zDEC;`H-;tQcNwFF1WYRb9CUm(4X_JUF%J@hCI~elIr`w)ROb5hYcI@Vxs&B32yhw~ zwpuclTN*%^-0T!Hc18&oMT`)8y+Hg-u2Ql-;o%&>>U4W!$;grbroLN?Y;a9vS#4uw zzy(&V*8tqKGy5zd)<{^WKel_HEq@$4)6@aGXD!NoI|S1G}pSQ=e%%M_3u1uQD zv#RU=(R5d|AfVd~c{FfN12&)?pgNmL-_jO|NrFRka9D_*)fluV65UyQBEuqw0m8}F zfJ)(%EtX&^oT^|toVLZC5-@hEC1tXmw#6Y^gY|&{WZ<7sp)mpUVK<+}cw*?NzZ(D6 z$Q(X$)Ep)=cfbpi*}TO~_O0A63Y$(#=@56GtrwH=C9%LOUzUgjy9X!XWR7#s z!sBJw+QO1mkYy5lKxYK>+UWx$#|m}p;qI5~53}leJE1X!z)*%WQra5VLng2mlABC5 zBgeB1IXmjq?{Q(MuDz>dqFc6)&~i)kNIf))t{i6X&48|MRoV!%9qopbMZB=EJ2SzA z8s~+xyTckp-^x`@SlFT>^@o_vnH`#?r2Ji#xCUf_NG^k)ujKCtrqVO6IKo)e!?&Hu z0AT~#XXkP1=Nvdc!*lF}$-RqFJ2$bNy5+98PROF4D&Q>rR6!{uj8JPu%KfmJZ^E`6 zm0&XXZ|rz+KBI+4lT1jrnh$(88NQ~%&ZP3BG}O~Al#lkP&}UT1i=(Q`Q{DYXR`vq$UDlG4YDT zmI$i4<{kdr##q=qjC85q_I7_kB@RE9NE;FR-rD=jEdzGO0Q8yCJu#!@ked~WNhz|L zvUIdXP_56v{fC7fIV`kibf_uXUkkogA03l$Jy8!|qqnfS$u8ez+z{qoH!r)~oxn(B zsPGAN%c$?!r0>rSkM8@6hlOt0xqH>9fo-WFlN{PTXtOrBj4c@bHo2KrK9Nz*Z)2xyu1And%94$}ZH^X}iK; zYZ%(8yX=ImLsFS$+2ByS4G)FLD~O);l!G1|x~8IxLXEToFA7j?hFu?Dtc7$%*vL=y zsT=%EjdRAQ<3o zNs31mci(An`0^4)w?#u2V@tyK*TQxp6Nyi|wyAq?2>s*6jsB!X@ofW$Cl-Vy9?-*` z{o}OYmorW+m!T0fu|SHJGAYI8M{oV``a8{cTD^ZCK5d(6@x95zYd>wy&lf)9Ba_ER zo4vkC)l0u*ZqQC*T)tgV15QpZ;db|VyQt2nktlBM9*OFR3PiPxj%tsJL~&R5NUpZ1 zL=^t&k*LYlymt!qgs#%11p!rAJ5QqA^5aT#e+) zM!iHa*CSDxh-(=%W8LiW8HvVdOUW^sx5KNyV&363oQVR!GIjaP>Ym-^&D{rd`R_JO zRVI2EWdb({>r~}Cs;50}Hj$wjsq^-O$hp1vmamz$p7+V>FTZZS)bzdIk#$=1 zG7qv=AN_{;Vc`z09|5tr82OIsNB0`8s;+i_)6DTcH+;r7%}J)|o*z@DJ9?Oh%6#No zW~_H+^>4pr#(QTDKlv?_^9X+b2{UIdJr~G!#|55=3qTCd#}S)rcOFM~IsjDXKWVP^ zK3={1NplvD15cWjJf?ozJ>K?hvzpM~f7`6^UZ}qEZS$6~@huuGXhd)%-9+xQIGjvJ z^?gs7Q#%K#H<}&ot2Y$G@kdXYQ>IWFJWxq3_JIJht+7isg7i_k4)bOE%%$8Kt*O}u z6{W0On$e~A!w>H>v&VU#9{#7_n5EjU+H=VK!ShxQx4mM9TPwe7Go^~Jt*GVtEhA#9 zuo*AX$wKP|4I*dMG41w2qwbIQfFXf&$Fpz-=89U*b;*%##~=cqRoD=Xz20q%V4l*f zoaYUjN~7H=aU_|EZ2KrbRdI_LDaoD=+i^+Ao(R3S+xAPs#`Q84u?ZE^X*(|ob#EIs zAe*>4E+=b;?YwLQ2huEy@-BsYk;u>f*2^~Z7T`?UGlo9ODM;7AVirM~Z3^@BMt8k# zCe8lCOszet>(Ff^s-WXCygS;_he4BJbe?qaHpYPPWg&{S+IR5MhB-F;$t++tv=oUw zdBmvQthw3~xYoz4A=m1z6hl|xyQA0fHEl|0H93QG->%ha-e4;=owe; zP70*}YvQ7SMf3+H0}k5;3?(M7h-vKBeCf7I&*T~uEl`(qd0#Pzm+PG&f4CNs;-?yQ zq`NA~T76!0z}mRRQUDOFq?WL=1rOcH!Cry50rouFJf}LiL&19ecou4H96kz-&{cbijkO6Q>nKD8P05X_fyXs21C}Ob>IkGF zBbe%nKsF;qK+|ItiBUBsoeIswW#kB@YYtD7Qlt7tM*^uH4pjkSjH>R2!;o$qHS)v8 zof#Q7h?-R1;0DgNP=?u%)Q*LklS5Rsb3jtKihD)iJrYgb20Lvlar2C;*deX59*wkk zCbPm-hq#4)2!8}+XXii!b{E6wVwBA=zwkZh^*Bs%el?WM2xT)uS?GB*%4X^)o9Tz& z)KC_(Oxidm$~x0V)d>QHn7Kn*BDPG7p_Ld#b;q*@dLH}+64_~+PxBVFMF$k|v@Wfd zY=g*!@;)!+J|N4Ji#VBdQvn_g+-C-m8={;}q(bWD&*|ZILr5jDrjM8?O=U{uw?PL_ zxV5HyQ!O8sE%i%(v3=#6B5l|u0?pVN=r9jW3Y8(p601G(3ZbV4B?nN75aluJ$pd`+ zpeg_=VBFz3q8^!gynZ>BW0fpYIlf3EVTJ4v2fAXejKvRf?$RrEnX8Hpi2&Xsx$VYK zGVgA@Qc25>fD#S_;ShW1D6UylIh>fZs1bcA5hUPo$(P3V$ z;a}!Ho;4pNuOzGOkN)(+#+tj~Es>&9LG-Qj>eUfWJ+z!bivGeTQ+5oH!81EbGh< z%TJ=GW7*XT?gp8@FK4kJCJrU`}1pgHKV{x zyFLSn>`f7H;a&}LMyY6%d)Gd7Kv0}?CGHv@08#Q-~S$ow%@%7yOqj{$U z27Q>Swd2BEPLtG;(i<8Hq_mQ}h@ZZd>q}l)7L@Q<*Y>$>IZy0?sC4h{?^~;9=lwq$ zj#~XyJ*a)R&GesaNdKoP@2Fn7z+dR?s17Xf-{9R zTH>EW)W464TC&u?l&HU6>T=CJ(VtlT$EE)9L^YjQi~5##QT4x1^k?#zdXm44y5D?~ zzr?#}cdB}0<14rQ1571f{s8Q_{7E0dKCcgEk|}?H z_i5$?ke0_h04~L`7!4?U!>H`jl)W>_J}b^Xt)Bf}WlxXFjs=2aImtdLY?`X-Z9V_s zk-fni=-ck{2ddsOzjM+9N-n87Kkh7h*h|q%slMn@rLY0K{6O`@W&V`Le^9uA&$>KN z?O*1fS&^<^iZ&bR`su7Aa@?c9k<|Uu9d=H>%MQ55sa5-@>a6jjkc&GtC zm$YtN56a6P)^tDKB#1(s0}O||ETl6au-1TbVm9y%wZrO+`)&qa=cIpoM*6p>-oL#1 zS8@514(`>8kd`;Z#gYzcQ3ncd6Sd-6y1FH2VKvCY%JOySg&_6F4Y_XPcOebSnk zWZFod_PIVy5m04n9RajwyrfV2U7xD=EJ#nPjE0I@tw9BA{xGOKQ)~6Pqgs7-1QO34 z)#!7zMh_fgrqc1jOBOzonaUmoiPxH`=jc;|R^MQ%<%PC^B{`N>c@&vqOd(r!0p{Ycg}aKv?JcPkrWmV(hovd#k{X2pS_Oo|(@}f@2@X6DZ6rAQJR=td@+H9p zNAWOfopqXFQLqRe3I(*F=Ww)!0#i#=DBo7ejhq@19QD!3M{uA}o$y7DT{Wut#(MKF zjx_(`QO&&2s5kCz@}3~Wp=kVQ=2k;u!sIZ#7m%`_da zrKX)~0d+a;!re^=5h!1#iVHW6KvoL;~q;Q5i2q}DyPKZA2 z7y!7BRAxGJ+@XW0o55Ert;c#?i{5WBTshh-N9VFFVbfiKZW3x5!06_Vr@J_%M$ff+ zT>ZdyGyV+XlAPtcf_Ac$2iPo8wuCw~Ubnjizc!@m>sVL^i+kg2JOBAy3jcdJoY@vk za96GO_S!To!|9z`)*8JsTfJnb8CzZWh-q$}V0-y;3u{jyM|)UTwd1PmA2mJG$LMPv zSF$sg7mpX~>CNcP=*`!gPl4+1-sN{p*2$$qr-fqJC_q0;Rb z-FIM4)833kwRXN$Hyo$Ax4D6|Y6x8T)~2gp`ET6-xlo1oac+R28{qsguRcJhit}}^ zWN5v1tcEt;4Q*nN8(P4}VBFByG{OrbW2Xo&3@qrbKJ=)Wn(wvaRE*o#1(CQI5Q&Qk zB5^rEBrZ4#BGuz}n6?hpa9Y^I>9_!xoC{2P4xL|p{87`^tUkp_6OFZnDK2CIbv;ajDuuFYHtEb=}6M%;0TEb@|T;NgiIK<9J?^vw}j=(k>x~fUIt3F zZQ5F$Y@d~Ki4%Kx# zP4m)%>;Tj`@{^}Ol*WsShg~HyspHi+?kW?R2(SKc$mHZ#qmYUJ0}9dn^}ZTyoPRWJ zR7TOptIjC%5zi;rU0~U0M%^SQvqJ2vPYSoU*CwUu7z>gWw}^~7leL_-K^{AsE!I5p z{SLwW<(91}TKHxkoclHd(07VVZFUJcyIFenpGL~&5u|LkZyt@6;*-b1O5MW-VrhZC z`50q}PkkbA7oyflD94&UL&`>*zl^reBJ1$(W&W=5ID0o;?*E~X(TB=zAAWYD-)Fp^ zSI>B-KU1GVYtXu{`ygty@=kw_d17z%^Y8Sp!GYfPE`NRL{T{o_I$8D!$Clq;{m8rg z=?jROpeQf2eEc8?P)q;9bBDXf@}gr7J?@ypco$+AmY3?U-sSIaA|&@^D8I1kq^ZAB6P0_z9OMylB?b2ayel@&VnbC?h0r~6HAOQi$Joq z%d+}!8&JzU_6}Cu*RS$U!axRY??xA5;?RNkM)>;n4OkM}$7&230I=_i_IWysxaJO# zu^@kE8p@(RY-`&-aN8)vMv;qais`B%r)NCUHeS?ePMowQH)(^h8mRNYr?sdig)=Ap z+1i}DaCtWxd^dU)G|N2-6ht5_=QusPcUGj_4<8kC_bB59Ht`FZOn@R0K;l6q?K}op zbA8ykKq?aoV?wt;(1!0*N*Fz3AXg@@u`|v6hr*2 zZMu7vT~t!vaBqM`A<7o6g)Cghs~{5kBFCxXvhI8+#s^)=$M?eLFlHGX{wDTd#`Em3 z+IeQ834ew~W&I6Ax$PksVsTV@!{&_+ei%@RyC!YBD>B+t%AHj{yrig|LPcSmD+H`l z%trSjRych&$RUBcEfOjT$0$ zCTFGJ9C;g~Q#RlhXn`;_eTG%W0_g`9QNoK)yfenE)n=9HXc9aAnyntcuD>?d1)&bf&ydD$5~sfkloyT|Pss zzSHgq9Jd{ga3^*ch8(<(-e!%Wx5^~9lNVSWl{ATylB^j*S}E;?^~7WTA(m2jt!m#S z6ceV4HK^!d9#dZ98Fa&4_r;MQd6S{2Y3}9v@K|`g?{J?KQ<#+ z-OJLPJ5|UCrg42yRu^sOqDg(_{fvYcuAy48lqo zU_Z6O4&i!Ufqi%kVxAznoPo|0L|xPgZsPhdu(i23zY_hAeiq%4tr~=Zhi(8*({t`X zIRo&Dx?5LT9~@{50HRZ)Bd+;zi{Dh>MbFlIR3d$>nm>VQDe>83`Dt>|GE8trRtPiI z^M2v?R0=WUlMmi8K8ON$EO3M`3zBK=V8@S4Bz=sheuU#C#em_F#iYw0s~B*k8prtx zfap}(W&J{v#~lQb`!t7Jxq{cTf(Ykl)bX_osr6xF(zC3SzVZIyF1PNOVNa)xUn|+T zBwt<-$mj&tgRqhXI5x9u4*)A;&eac5(+*RXg&Mf6Lq>~)u7)p~l8Dao;#fc|n8XRl zeuQv#&Se?huv@-PSZ8sK<;hACz=1xCTrokqC9s$#h1pr!Q!}8gs@I!yvt)v_xR~n{ z{u}Why(dcq39f-Q@Uggz<60Gl0J%D5CCaAq|P^m6C{cLGVJFF^faT4)00 zIl6O7hQruym!Ny9a2!hY&=?B9*JI%rx0VP*I?024H6^2{goPzyShXhmUe#3;8PO*z zVAL}5!L}Nc($L%h zor)Wxe^C?8VIMakCXd?W(Qk$#M#g%cR6$fb_X4lZ@d@`yv@yKnZ~fbf*gKE>t$#<6 z%`N-B;ZMw0eppDiRm@Fn3X(&7AMWHRO+>5`pzH_HF8pB zxw$bWJ@gcOF+G=9w}LgMV#i#K0}+!;lFU-4+$EeLFu*J z8#U$k=-5A zn44;6qe^b9H0Bs=LI1;XH3rK)5Ahz4gu(JbDa``oXGuk@r^MG)gy^yZ5X@W~>1yw& zHN`FKD!>v`ziTH8OdZsul8Cu+52w^o?HVDUxE7sb^)IfUL9u_Tn3YtO3hCV-lGTw3 zO+B~pY?R)>vQ;6nL}eq(Wq5eop|Dg0V`fVTxT$glF0=$wi+vAqmH*91VW|4X{)cR@ zg75$H&jck{Hduj>5UqB^dsec#T`Xm@e16znE1`bzp`hr)rmt^y3jevR4p!kQ=o5{8 z?er89FHZlj+1DhI>R@f8TI z;?(DOL{iS3dsul?0VtxJ0H`K}l8p39gdtYCO4llhTXC%_H^7GtfJx&_sz4HubE8`k z(TXeI9_zK)+Nq{2SD)pmOl}fuoa}OY(j0AQMjfd!0&PU0CWqNHkxVfPjTp(xG6@W zrX*rPD^y$UWjBR_?$Onqfec0zmPK6vjp#jd6pvc-Ahf}}_5%16ZU_T9R#U#$dX& z*>A3XyDdGlNQIK~IzgAOX4=z}Ph;v47oxtHO-mU~kFwZf<*M$!j$Q9|3XL-ft8Ge1$e96NMA zpRg@vhwkP~_geL5WV+eJB0NSet1?GEe!nUO*L2|__YX_Clr`tW9>n-s36I$7tzvwM zN>?-0#9+{+9F|GB2SrB68uPicmQwz11sTiOXa&9x-bx!%Lr`1j*2QuOsB#-(3uz=Q z+XUoP2a(uXUgDu%JDuKjH)dh6w>D7as>2_UJuM6C`txcCcIvbHYvY8{D;X^X$f>b`prSFB{ z-oPj)dfE=g6fj8b#Dd4N@d`;zcdS30d?n|bm1W5X+%WgJ z@{Dx&HnXorHg$S4Uhujq$ZW>|BlC;tCBY z;oL%U(1ct~ngSqDR=l#6bsUPTW!c~hRQUsH+`Ekx0NR?wNgcDTjF^`k)mkx2jN*(` zH@WQ)rCK-xx6g)uAaq zG3l;#+J=L}B{opF5+KDSh(o#vsgCh2Nk9xMty|L6xAFH9Yv=zOF!7>6z!g0#nCblA ze-!=3N^VnaSAIkFiBF{`p4T9SOCmH~zFOku`S>0J-hPdbWu#5+49UoI*PwHT8ih>0 z#3jt+0(OcVq;U1*Pp7+!@^JXTnfOfgs!zkLG6Hntid6iw8JO3^gaEQ@kOJ{N#(Ik{ z8Pr?+8H+)EOx{R`4ft}8`0rNm)`p^4Oqt8)batTqP&*xboXmsD(j30EQTOz#Lg z7n^(1iJYw2t*O*HD)e@QU8kv|R}e5e!tUyEpGi-Sx>)l4=JG9f+fF4e+fshXsFE$z zM+6Iv-Ejs5qAhoaU5m|~YAddb59d{@pGhyMxNQ~U&pW}pXt!_%Faf^ZzQb6?N{*ZMg6$Tx8UQ)^*aP=gNlmpbMmh-M&n))t4X&&ns>#i2nl zC&$v#eDFMd(29 z_epe~L-o#)-y7P>*ITYRBeST&wQr7xU2cj{)Vk+HMAmyPPR51ccp_0MV06T0mW z0YpA}?q))l9cU+IJ&zbfjNL*vM;mAhu{ic6&|xXEuC&_qioc%r+|ogAeNbG{J$uC) z{yXSd9&sqd>BIeRjxQ6xO{oB9YcdyKdnQ9Fsx*O#;^5%IYe54ZNUf=Ti)h&N~p+J8aLKiYs^<; zw9)L^fcdHNW>xE$H;Bm;4>IM_s?4>$Lb3Xc58?pfY{g2b9*HNT1O6M&)e^qXsjCB> zHO+h;NJCl7q~`2)K3tznAD692G+4B(WMXsGk35iG!Yb;w9!M_^MhH#3HJNdnD{m&` zR>+4h>P-){mj?7TOZI9`x6{f4!!ONEuksre7e`cAAD{k*@R6S;G7Gs>iuTk(vaO^0 zU8p+mg!Cyz)@i9s4yeN3>bp-!-^fX?7fwju#p9L*=}+*OxX^`m#E(^r(r@kJ7D0(v zOf-4g5QTI$TV7m!a#4B%C8EWy#Qw$UI~dy4OVTIt_~epwkJ;=E|J{=G)P|K0HFGCR zX_NbGz&>_1dXwPvpL{%pWgai-qB~=ondPTP#gVp8R)=rlYMVRne;hb- z_ZdGqTt(5XaRdS~+EROs?v399jiA`qn| zr`eWOBTv3!+`d40RxidXjRyZ~cPjcAv#T>O=2Of8gY%Xb(hM$HSlmS(?*rX~m{DefXE*mtR;m zWQ!{{Sab7|Q@;Kc&{rSUZPFk7u&KKKFVkI1i=*~TlFq(r9=Ey=YAmFeI76BXZ>T7*8=3)wVVf+ z+2S@4lWs|Xb*{ug%$&kdjTmaWH&m?q?i^s+iK1x$)IA?a@2y;IJvR8FN9wsV8|0Zu z-gal9z|l9G7e_d3Nt9OSsI(v@&BR|DAI(r3j?GY|FTzkeT@{rY3zfgA6Q(utwJ;o` z2~^;dTJWq*CqWe3m~BeNnOGaa&!&|#$2MEVHrw2ee?I=Fk6aap!Y*;LF}uRIBz9KEli@;5BO#gO%5eHSrxbb&Y*}m5!=+w z74#u#OtY^x%COC5>MpA8vT9|)75*tSU_SYLZ##3`4&UTBM;!|WaTc=#c7_xRJcMGK zPD$gQllnYpzDvv&QUl$;*2T^46FmCtGzqA2xg!z;M_|SV2&c+bzk6$X!s`i1V?#5R z`i;v<9x61|m~5S!=de;hK|rSf@Q(Fs@m^7&+)2Easw+O2eoqVQsRuIz(S6ZT(T%J7 zKA66|xhd&XtX!0;Cx0kydy0v=*D`v90#o(I52a^tcN#45h3dG|(p4Uho|gV7k8@VI z$G^sp_nofL@VJ%f2fVA!E^yX6g;esRgJK{*KoxUgFB&`|{O}Rs-A9D?9T7g@!oyR} zOgCj%a#*k?UC4N!u14ple>opVPtGlM97Ipo1`zk?X7NSdB^Q_!$64ai#0_8e`t()) z1kU-&xq|(IU&@P$b-$FBZOHZG3mCt|SYBd`IHKxZlK#r9TWW&K zP3M+c$ceER4*&{x83#5G3tiQR4*FxO`AgGv7us{m@$~SNOVfGRT;oV{Ctj9*!ZoK) zXV8AN<;wI|!aemW<>K653%TYVaiQw(uSA#Tk~Nx?)MR23AKv}0F=pm{jV4!V{{Kq* z68I>J?EmSWuE{YuCP28tF-d@g`xG!Bq&ZXsMO<}NypTgVR&Kn&2@(|%G$K%_0Rlu3 zK|zeh12t$=P|zs2u*xb33aAUwMZs0%|NT~XO*#Ya`uqHoPfyjW`+W84)vH%Cu~iM@ zoC6CR23FW?>(2%M8AdoOX%sfiH8atN2?H5zcYpm#-&GJ5(d11ye--8-2*nc2eMFAu zOy+Y6nP5%wI;2T*IJE?;!E+bvzq8(R$FMuc1s4c?y5n*mmIpuPS615SfmpFM9-?#V zfY#SO{MJvo)rf-tL>tK^S#`$@%mhu~iQbUDoZ3kxVLTHy;JJt=E9RDtC<%I{c1#5G zsXu0vR6C>kSaMQIzyeMl&{h4ttqRSQ8jKK)8R)rdw6`!G&1V~C`U?D-#9Y!UyD5>eC%w$Yq@*&M=z10CY!-M1s_*rwO;U+zdLPG3}wuD@MIL%v%)Rk(a*EE%CN=oUmDe=&53GM=Tn-V(@G!;XsrbGz5Y|p}A z%NT7-Wa2FgUkp97EKCTh< zhZP&EH>f`Z;`xC}9WlTZ=IRN~q}shfM2n&J!h=f0xLHN9;j$SnN^Bsis1J@vTcxvy zvR%Ba^_v0<6rBlOeZH`8RrOJTg+z`G${C+QWCXNexS7R4a?1gZV`+7j8&zY(qbs22 z5$&$n^h#j@@X#XX9ye{!Vdry5I9|Xay+ZLKY?5#pn5fUmH%{eiWm;#+8@iO-ybpul zb&6*-05wMcb?6q1KD;tw{ z`hbL-Zt6-jwd{MV)h2=#uhojU{AA@jtSvuo=%^0e?$5;G0$SC`RIQ2hATz=tk3D-h z&UrY^avn~D3oj#}j7m!>;_zN-5JJMBcu{x=CWYOzP)r{yqFAobSLk86;wnC6sGJ5+ z1SFXlh>I-Yq52eeno~LG)F>Q>R)SDlyB;YG*pg6mz7F41k8(c+XUPV|Dgau4vk3$;zk#!MamsG@2j9sBa zq7!g)lT1$RaJ!qhr3F4*ZAEV6dOC{SHZqn{=L#rBB!|+Nca#uWQq_iF`wXuM@ z9;j~S@y!$jIjOSz4Lq&fVdj7;46G=a#e9}Ns{`bIu*bo0P#YUwQRn4O zu$G~*0M}RPfl&$G3KRwyLHW_oR(gFkgegY|vQZ6Q@&IrtA=PS`U?{ndIshwDtcYRA z{8TF=VxajDI-UkKf?y&it8uL4{dw8i9LFse(g7cvLew0|P^t+OWS~*0Ta%JF;GGdE zG=REBns|s+$>%sqHL=4C!v|Nh9HpAc;Q(eGhJ3)0{6^h38tNHcuCBl_xC%NLN^fdf z8cH>5->*b5RQv8TvLMO8NKvJl2&m9RoQrG|92AiN6j2|j0ozN=IUbV$Dad!~SlRh% z$7$lSM_nvpY&l1cBw}hBj<1-{W()PRN=oU>c6P0A`+d{ zR8df zD2sKdjCI(h9Z*R~wE^wXfSw3>$7qX+spM$u3sHV>qP4^yIL^RN$%3iYxhQOskeSo0 zbc-$J*w{Xfjjf_=#wyutn$-|((PNtR7TV%J(^Ojwnr>Z;W{@kr5#QL47_IOQ-~K7)2-#BZ-d~#z=99Z<{9?h?Z~@; z#9wl_F2q5CWL%WM(QadJ*M9b92Hpu>{85+na0#t0|KLvR<2W?pvU{v|jTCoZUGSx< zk6b>_x-j{KIUw!Z_Ol<8lfiJ8|DI=MLvW?RgVqFUc|bSd5K7K}5aKM4%MA}&$Du&2 zbUtn%W$t+MrF*8RqD`TS$n%^CM76%Fkm zYX}_*xnL}U08IsMNC=h4u8&(+Ms{sME2K~wHQ=aA04oHn*tY(jS^iy%zF z{%^9(nqTLHxoTxNWwpD#_4z_4%3ob(B^$BtC>I39hbX~xZxF27xqEIMBT5DXAd*u$2u(*U*wbr`m z9NgLC-5%Y;r6Qhw#_CSK8Wt86r0+T__l#1tXw3Oh<#Pnreyc=(;n{IVhhO=| zb=F{^?)E*u&bo+f!a*-t>nuo)oh}fKWXfx*7dyUYZNfPZoVq670+j<;UGUtbY`)og z(REIpYee1s!C;xJddq52-s^Si2D2kLYM51tSk+1ox*da&7iX7gNVC}%(P4M__P4A% z%#Jv{LM$p;EBso=BD6?DI@WO&OHO;o8fT0w|MDGcu3>yGZ+XwU2*wNV0sLI{-ilsZ zCMR#T8pD_mqX!2iNRUdArADGs8>~m8<5d;)!fOHAOa2p)v)srKZEz!lZLC8ua_Dce z(Kf41B98q`40F&*GQ7=dE|eMrIeD8?3>dl?RfbM!7m!`eEJ4)0QLlY9UhyKyi=ct(Z%*M}#!|J-YL@!4UnyB|94j10Q9(yp^u8+c zIv7wCsakFyHhG-GYJaJVEcZJBfAa(^1kQE>&c$eQuPu~!`M{bV7g+Qy`u6X#@)PSq zKERQ()2h$6om%a*+QqGAEH@{$jgv6c(KI(@Dil8HTMu?ZhEcNCJ&jYdY`9XSf`(UtI@mW4D zbNclqtU(4IVY-;FIbrEzfhYvW*FTVODW@Cpr3Bc5B>5;2TkI?XucCZ;Pn>RH_xsuho#3-K0RQhMlKgJWqt7 zZE8LBqFyn}RN9c#>n~Z552^uP@;zPG^}(c@2z{Rd5PWjK4ja$X@fE-z+qcclwnT8SXS;~Jc;Evd3|KN z16D(XP-b0kNK6Mtm6H!xjg$92#1adh(ZXrz0c*H0#P`4jV&$cEoo7n{45 z$j1&@Z90K#O0H0`X_r0$(t-;F5R_^2R?03!;Z-q+nu4_uhrJVGZTSWbTfJ27|Hf); z{;*V@`o?M+Sj4zb%J3m8J7M3$DxiWb`7 z|M@P#{EBmWHqJh&F%`nUyJj#$Y6Jl|cFM!w$=Kf@`TQSWEdk=~0dk=|0Q$Y%y zQa24y9~=@Vf@>U&bADz+KM9wFRkh@hqT~vOLz?%@7;3g~!hor7AR2b_+=k*{hzp83 zBo6S97$TK)hQ<|cgI8q*-9zFekV25SkGz%@w`wI)4uE*T?$rzjVc}6cIxTKh31)aQ znh+;U`t-?1;Tr6Ffw(oF(NMBA={%)TnKjaxDcKqv5(||DunN-_OqCPBS_KI(h;5aC zcS{a&JkGpu#OBtwx^D^S$~xf&(8Wf~S_fw{U}<^}Pel|6I#1PL*$QMC6oN#raX8@6 z1O~gY9qo2E z%&<9uXpgKQv$s20bpDr2Vn z>6ld#VS+lj&`d@Oz6GiNS&q6qrTq~F!P*5yp87~YBV;E>`^wlQR0iyBZ!an5e2B7H*dz9 zF?pjJe=2XKH{nm@&8(qq;*7lO{eR?5p|d|zZ(Z$cbiG|AjUlw7Q%kIsY%C)KR}jWQ zQo?BtH+k0;yi8IQm7s-`*ih6)eTn63s;$kYpw^8tdHOr6Q~A<=TW|V1?b<}riE9py z4Y15sYS=!G7dtd32xPfkrhMf6C#+@-+}AvCbIfJxt4vNyknMpJ^x^HBa@pWv#)gqEfmOD80CPlNUM3kc0 zCJ-qA^)<2$RllnsY*{M* zR=^&?d*a-H>O_>Yo6I!|dEpdJAwIl0#=|VXpImhHx)>%tNovU8vPC$x`N$YR)`;ph zU@#-uwzWtkrU25F7rkX&k;fG5ifFSSA-hULZn=wepDWn3OB`5c|8Cu9+#sL%-MSrI zuR&Hu>sN9xYFg06@P$26_Eabl)dJTo!kb;r{QSNeZK-ika!CdrV3^srXL?E zI8e`@9`QRjfI&Icg-}UBiL8zj4F~P}uENI?06sYG%qw^~f|NM7H+?3CIUQ^>&*!s` zL&lOc7B&(-`!^TBb#J?Rg&YNfgcN z>$2gY4KUclDH!=EyquFLUZh-7fXBHZ$R*(#Ts`z3j8IfA2>@l4U{F?xa!GPrP~;I3 zUkHhgNp4eR`B}iDA+fP3BOpQ5)Myi9vWm$nq0QxMY2v&%)dJtIkf+l`&%jnR7I+%x z)e&0^^V8y&>WN*J`Jf?xXe2V^4-GNNG|;n=2Z(rqcN>Y(U|Y0nEG__JW?EyBOSP)O zH!Q=rn`BQOa4n{#(H7F?lu(A9WFn`;MUMYz`m*HVF&WQ zIjB5;-r3lI145PFz=Jz#6N0J|hJF{T6=KBfY# z$W-gZ4K4s=4`LKrCKG*Jx}FUOpY0dEuaI{V69xj0BwjFum(Q4%gW5QQDGV?W^0_QK zR|MSL))E0KWQE8|8VuCsM^psB%R#O`hKP;v02Edfm5#H(b>kxy@ZiDN#x3X};u8fY z*(DlAwd2plyO}`cU2erdrSxc*2^T=@4L5HS)#?M>0u9NEV6(e~vPnX5e*sODl*B+g z0K_H(tQS=qKLsI34uPJB3syFWR)J1XZRP4c(27|=y2c?iPmSVS*@tjy50ivNu=@qr z9Qmk66S59ehe3?b2TEZWe9k5h8b{o+;5rXz3wUEKJFbK-SXEStZBV|Vu}J6hB`;=( zJ7IKcA|}OE@THTd{qpfS9ohP+bhF5o(GY9}q3% zbuX)GsV)-@e?R8aj9K+k~X<5Q?!C-b*A{S*W3re6hQTD zKOR7^6M3Kv!>*c-UAJ1px!~7wpqvNB+yf8FrrDx><7ff~f>F2v8ltGe0{)}2MVAKN zlo-TKT0WaCF1J0A!O%f;1EUIO(E&sUAEUhwyHcfX!XJ%@?!odIO{QN9kyHNxRSBFS zVwTi|ERVDhos2-?0aSXEe6gjtzrpdZXlD<*tLH&F?3RCNCH`(!ttfAsBknWG7q$^& z{N__ec~Kj2i{DryM|KvMoyV?L_vKQfUCpE9XwUOx8H=Uvd0wu=V(Iq0{IRq6?rfc! z@_*>mjIN@Wmkl_nt7wNehnIB~4gU-SaCcV`%K$XGiF5cKy!$S)%*sJC^-@Y!q*_!l|IU0F7!k80pZG2R zED`a~H10Cr@-mZu``q8DSFJ%Q5F#*eYHq-Xw#RKe=Um`Z*E*QIa;3lwZK=|^yY?hC-FoqQyK%8K2awOaPqFDeh@%}CF5gu7VaK22?c}# z2b?Y;eT*Z#G|kmQaTM(8XsV#?Kw&HV1)7r1MAcE|LO2<3(Sak9$@V~r)!{x^dy1My ztxWagaDrB=!7!8o5K3&oyaM!~U)2cHl9X0B6eBahi{Q>9kk}HadztF)xZ!6-xpxi= zB@^US0ZPol-6-#L4ivoy9Hy3>)>BNvLFU&z#UA4e`Np}TF;vI@<6O})QLRB-R8N)1 z=ZOw};}O~8JaLV2vs`$d=m_X72Wr|!w9&h|#7CJG-hNRkf2h78 z>c9l>V6+>aqm-0>Gy%%%aeJp7Fr?FXuf~KSFSA9*dd?eDl#4*U5e|I&g>tbivUOwN z`6&68Er#^GjtlNk3>#w#BuCZKlZ8cV=oC1aQ@~GjoWLVsAb)gylIQmmW#^r-xZqTm zcnCA`e-q{9J^G7A2AXf!0MWcv*_RKpZia$n^uojK~tuK4e_-?wCKg*Ka0*$~W)U1Q$*FBya>40cXk{2GH6MrA#<4 z%*1^F2UIGZLXVj?z15EnljZN>2c4^f)g%VI5Dd$|T_{44*m!=PMH=3Cfm-px*d>Fs zV+p9W=Zt*mlNINUZ`O*>xaLgEE56KZCY^tlB;uF^o))bsI=6K-Edk9jI$P8f5|7oZ zWxqk<>NGY$*g|Z?Ma%BlCf5%Vkyh?dQA?LgrMN>yEnn_>C5K%kO3qxM;=}(Js-xpz zQErx$$isuh4RJhsv0q*-FBu|+#%(6Bd3AZ|5OIS)P)!nQuN*&2bP5&x!0E{)xrnx> z<+@>_ZF&jp_UU^oIF!QU>{5E|mER#kpiD)OZ(b^%3KafWlh7jCRYIk-ojsv4CxS|- zQbj1wyiDW}xeU5OTnXe-a)oG>=A6lLU@MuP&?9p96;P9@&%YIxL6Mi4rZp=LLLrn4 z!^zulh4T~#QYmgCn1sf*lt(`ac5~5&DV)8>Io+|>A zRXk7d{^M;h-aM>EBn$a_6;?4Sqv^-NgKTc=dT=>`g^*8TMwU`5rHaWDaja4qU{dl0 zX~p%_oKB?KNFWhx76i-8kZfrMvU_7oF2YGbj1yz1=ow;YB{ZuZ{Hw>rInQM zXbtp37?|0Q(JWiuC|c>EEib!KoCa2$b(6S2J~l$+L6%aUzDe}wfiUS(XG;;iG*jj* zOWkgQ@)s3rwsHM`3ty?imtQ{;+?%+27GOG+E|B&ZC$Tu$+>2mokVlu{mXsoyqbeyn3SOX3jC>LlZ?$G~s&_MGue=rzhgIQt}Vgo)rWPmdp;5#6zJlSOOf{mT?_o|AVYdFfP< zP~W?iL&$Ymg#lxNb zr-`d^t8UpeF&M^?Y2sq@2XIuTp{p;RF8&I)C#Q?)FfwL{5ly(~RU197G6$m&5<|l+ zogvzwUT@D3-BW9b6HtSoOX6=6Epu42&*RJ_Tyx0w*3wSY4l^Ep>MS!5#(}3gPNlj;o}uC_6PtA$Pbi|>Nj3qp&y!y_>&w#HMBNC;9;H{0 zl|Sx>OffChezifv+yp0<92H!A$lhTuKCLDEG?nn6&X0dYeD^tH>aQMWIs^v4@P5P{+}9|#Vv~;8Zxtg^jF|L*9oa}wO&Wp~8bT?xY8>DasNBJU z4l9fafVGjlFbb&P0WmZ{00)sgfHh&IER7>h@2KpB1GIROH^@x@MzwOl37*Rd8yLpy zSGx=qhITY3#w4aC`W(X^lr^;=lkbJ$4_ZK#hsw+l6gw)^$7+A3Q9dT;}!jSgm!N{BlnU^pZMrR#r~5|k{Q zH6>2TFcO9N-<0&>;WsGh^M@f7g~9O~q-B(nK9(M|lH#P8r-ToM4?0@<{6I~P6}U}B z<=8{)8V1Xmiu=s*VPpgr6i7=)6COGdk(di;dSG#6h8AMy0_~G2hUv=SKp0NqKuH#8 zWug&`0vn%KD|8^#KzU(whCX4{IZqi5*ND)(H$nJ#ICgD{ zd5DKOhi6&XzM;+7j8OF%KRioN)hVcxxa$xO*M_c}O34*IXG9;86p$o9H)O*ngM5_j zY>(Q~3ffJE?iJ%v4||>%(s+P6#^50e*hs4GYOeq>efibBVv-WyH%St4f-Vh=v*#)e zr(s`cvi#e9qGR1J_~NRP4zrKrs)F4&v_WpTPkfWJ^URijQ*2A@l$-At!;NjS;at&6 zaBHR8fWreV6sUEnmfs1)e6s8pVqE1ET#4*>icuqbR= zdG_X0N@?&TJ3k=0HFMR{K+~)bj~6H{HCUaL!a?Yvk#impT{|yDRN^mn70=aP1EikD zeIny)-8LA47~9B69(n*IuhXK1*v7$9Rd*d7@FXC0Kyi=qi-Jn>(` z?I;oj#s%_{N3fd&H${F_B<^-$9H=y1C`LC}s7M@jYKcc-Fp4pQu(rs@7K(Akh0=Od z^u~Fv{itYcf4?0gJJ(*c1Cy$N!@zYYjdduCb*OY4WbvaSPY>_CkBYW9>aER>l?@h& zY@n51i^PQ5MH4kedt;H9tW!B>v3S<3en%c$EN-AVaEW+C%MX!?whNv|p||4){Lgg}q_bJ2gYgk9wt*`KtCp%Vzq=HI z9qZ)g$HaA@Bc?y5y2ynZKQ1nVU)mE0a!S5Y3%cnEar4=k&i{NuG*7Ge9u?@Kr+Fa!j$v)9iJZkMS&hu9r6#i*vP)!?BP|zg&jB zzk3b<0;~mYX<((vmPeM0CUVkpArNN9@|rMPmWv$YYI-wva0O1kPRX7tM4yIVfxZJd z>FC7d)06~pj4dBpA@Wk6HPpEvu6SZO#Q@vALiDfuETRDOaMlA|W5g&Ja< zD@Q!714JP=g^g>G!a&>s(g9jZuag_9{WD-ixEUzB>z=c&#$DN6_vP!6YK9!Nz9v<# zlN&1IS%m5$E1shn=ys-4bJ@O3qzNkY zt{wcc7ezi0&*>M@y{F`BFJgc@kRSy(R;HGTj(#ijb*(%Wy(GFAUZ*^_LCyat38O}` zq>%_%A>QB_Q4UV(vkT<1a?#LSzd)9kqcah#7~RlVDG&t(q8lEUctO1qcgU_Uivzmh zo4g_-*!vZ~BAVhE>sMdFf#dyHBC&tn|B6_pQ%TaPP)UMA^vPwfiiQH+8@x(3*(^Hh z1dqKc#yau6X)xgaht^eZ6CH-4ASzE8(yC2r8{#e_PY%_-viX~0kkLf;Wz?cK5%qleAU!{%JtFs1h*nGzy4q`!Su2Yx zbn3wx=S$^{6(IMOZ!afqo}nB&B~j0^TL_hVFDIW~PDyWbO3=%x!wZHgXO6t`ZS>C1 zPK{U4euDO)5I&PBm12Oo%rEmBL9P^i9<=4`O*rd}&catJMXv6~cDfsLYr65*O3^Z< zHp^Ugct@-QTSjqexQF25C+p;qcf?W*tG9|4Nd;Rg zNJuAv9NCIXc3|h;D+h1G9&MhSuuWW>G#6|P$A50I+_eoISt>7mUu1SH<>nK&EB*c8 zA0kBZtTi;@;3r&=l@P+sN~&pwKkTyPeGqu1viS$1X}XG=KOPalNo?)Ev>+=PA%=Y* z9?gD|k*a-27rZi+!ae{9G>JdABdaIN}3}do=`-`W zgADynq7;Qe^YH9{0xntR)?^cGu1l*#=0a(fHy!y>mB@t@XE2xu9tU=act*2IGzEbd z4Aud*n1h|uxppxOoUmnl3`vV+virwkX~*J|6+ZXXIqa>XfXa=8?E0t^aDC}XdFo@) zsO_HFHrJqy1$O9Oe@(%h1;(x{L}uCc9vH~X{zNqEFEJABA`VAx?yyVDu;UoD4jeEm z)XNp9G?7yyIB)}qAFwcK!w*?mt3DCgM2IScz5N(cdG zQ+JB5%?}sIpLdGeF)b$U5^uri@~K#EJSji^RICFu>+ixIk{?V_Wx;M9|DhqZ!BB0m zMQyM{ZLo(-+zo=bM2_4oE>4Z!&Lp1}_c-6&jd{3Qe!Cll`V!fB56};8W8Skz6eexO zLAi5p5-xTAAufntV({IZ%GGkjKg1RGzSF342fL^M4XbuSC5}TC49DSrgpeoxA###< zDbvYj`)b+dGjMyS$(f&tMw^0^q0ylV|w)xi7@tuT`0md`}L;L~%WPICR{BDW8y z3S3V!!`PXPv-cjY@OAY0dil_Th-U8s+}h{spm*BqAPQKzV14W-(%7^T04y|j>^bhk#nh24yMs_oV~@HMyGN1D}bZrK)Kce_LKu9 zd#tw01IOggdqsBcZ~v;`ma`OjI_+xMZJZ+r!)~iBy|HPq|0yyf^MKf~<4`)eVBtLp z5>w=ybmBp08k)r&%vknAFZo>;*IBuQ6cBY86?7cb3}B>DjAMA_113*h1h=tT?Ej}| zV%#Zz|EIX9JD=jQ|L3q!soPZ$h(ekW-;sBJA+m=k4zhwFnZmR^@;$*+fVhomcLi?Y zG*toSgT3ygfE(!$w8Q;VyDb6~_~pO85d9)AKZhLrz{c*}t_uYinm z!v98&$^qzV0|WK7OE#eWnOP8i*}SEn?F%;csdv2+u>P|_LM4XcV1}9ogMEEwzSq+a zU-|r(qD$|J@6e%w7zcc(AHN5Xp|-XONEoe~;KK~I7qHhq@+prbcql9s+9{jw6L}XF zmO57meD<%%H3R%8#z7aaCu&tm*Ti|%02&4!VcK^hO9XonYpvZdbf;XoPh1o!16IYk z{6S8jt+ewT%Nl3D414HhUQBA8=c{Y))n0WWQRR5qc;|;ac+5Tw4m7R(8W{h%dgcIQ zzW_Ko>etl(ejB$$n_zMAb(iDzi^TTTOr3{V^~TxzoYX_B%uMu{B6rE!gpC^wvKvF2AzOf(>Wt)*(+n z+B$BXu;QWOXvxqjS#V&Z)VVgK*TF0r=M%*-EbjQil$BrCd{=pC8NJ?7wP{?<^&*wWA zUmbj77QYj??VKa9W#CJW5Bbty{tF+*UlcBfJ<|if3VUnW*$8B$wK2#at_hke7Z`#5 zO$XEs`4Zt=&-D44!PIbL!txYiMz$g|1}VhX2A_r>f}NT1FA1N8FAksfp9eeRY50EF z8D7J)VQ=7p3u53MfRSG$8sT%tAqsy3m{rmV?uXr*;TM3t8SVzm_SE5b#)Pi|obB;{ zANB^wKN?}1a;Slq0Y1+YVI}M+aWwpf7?|DRE%1wgy#?L`*jwOrfV(08Xo72EB9sF5 zmUuZ}Z&Tlhs$BVItFG7=mbW>-o*otf}Jg* z{YS&zH3lB$fycn!$^+lW{`V%xK?Fi2n&1N1y$P~~F7WshW((;8JPNxvy+we%={*kE zo8EZ#zc<1J1n}1AZ!r-FduwznU~i4C0(==hT_D07lp4MmcAVbS!0dl-lQTe;Cjbk_ zM05g#nSh?jCT4PT4}9A2`5JrR7XW+1UkjKmqT{Ds4TsJihof?~2xLdJe_MEShti|HT8(1@^p5TG|!S3}R1h}cke<18c{yIKkhQBG%DZmccn|R=K zv_OO<(g7O4A?$&h0p=Lh{>@|jGXN)g{2K!17}eqH0w&(l@IlzU`F+Fw4|@Wfj0xaJ z0-SX^f!|^GX80@Mu*d&DuzL$|95A~?$0yv$10QAoduwV=QQ_aGx-{uT!0nk4w^LXgmYEs&_!1t&vYA2HMxJ%nGSV+Y;XVGx6=q zl7|xm{Y-zNY!?hnG2>HZNifhQJQNhs4e-ARA3IgU17K&T%1?rUJQ{4sf5D9Jq=8`;r5BQU5~tb7sT6T}q(4 zy%-g`7VwpBx*C2m20k7G_e1qqj}& z8fJxCBmbzwXRwp*(E$j16WA342hgN!avgp@>|VGpU{*-`vqE0@S;UW==Lzs89N6SK z!t;QM88mztus8g#G2y?D3BMyI{FgD|KQ_?*-VDYffVYGsioEdD75ZlJvd7}PPwp<4`mZVb#N z-J3vv4a=@g0z)IZA$A~7_(by^3V`sBE~@VXw8V01I3IAT2WEJ$KbHos>N@-}gzpLX ziqw$rW|)WY37^jo2MM)ooz3cRWsjzT0lnv9b`FN;prAwl4xBP&;uOGJ)0~Lsz)rfD zTFLmg=Nzne)WQ+`tbzPCF_6!7^OvT9XVQK^O6)^+8K>nBvaDI4hw-EQzFA;*;tH&d zBjK9wEZOqw=7Bm9?U;(vC%{v~1+a4*YB&%BlP(D1`wYY82AG%QL+zsQ6)`XuPo~d5 z9cMU8>Wy?0JiRc9b1zIvk%{O8xJr2a!!dBz7`R&u+>2H8CYYxKlqY2dJ}@Fzof}eP zV%JX`q1Yh+Mz{LzUNBF=lU#V??j63>Hy$)6QH5@E4Ud3ZGQK>l$G5^f4<%q9=LO}C?7*wW zmvU~4z@x^l@+K_Edw+G)reE-6C z3SWGmkgpEDQLs0M*$!V%e1q@}$2S&V1mA=Bmf%~CZw`wHK2e82Sx zo$pKR8}il1*Aib(eEsoVf^Rs!G599on~v|o5mRpoD@I>< Self { + let tau_rise = clamp_tau_rise(tau_rise, tau_decay); let dt = 1.0 / fs; let d = (-dt / tau_decay).exp(); let r = (-dt / tau_rise).exp(); @@ -37,20 +40,13 @@ impl BandedAR2 { /// Recompute coefficients after parameter change. pub(crate) fn update(&mut self, tau_rise: f64, tau_decay: f64, fs: f64) { - let dt = 1.0 / fs; - let d = (-dt / tau_decay).exp(); - let r = (-dt / tau_rise).exp(); - self.g1 = d + r; - self.g2 = -(d * r); - self.impulse_peak = compute_impulse_peak(self.g1, self.g2, tau_decay, fs); - self.lipschitz = - compute_banded_lipschitz(self.g1, self.g2) / (self.impulse_peak * self.impulse_peak); + *self = Self::new(tau_rise, tau_decay, fs); } /// Forward convolution: s -> normalized AR2 output, O(T). /// - /// Runs the raw AR2 recursion then divides by the impulse peak so that - /// a single spike produces a peak of 1.0 at all sampling rates. + /// Pre-scales input by 1/peak so the AR2 recursion directly produces + /// a peak-normalized output — no second normalization pass needed. pub(crate) fn convolve_forward(&self, source: &[f32], output: &mut [f32]) { let n = source.len(); if n == 0 { @@ -61,23 +57,19 @@ impl BandedAR2 { let g2 = self.g2 as f32; let inv_peak = (1.0 / self.impulse_peak) as f32; - output[0] = source[0]; + output[0] = source[0] * inv_peak; if n > 1 { - output[1] = g1 * output[0] + source[1]; + output[1] = g1 * output[0] + source[1] * inv_peak; } for t in 2..n { - output[t] = g1 * output[t - 1] + g2 * output[t - 2] + source[t]; - } - - // Normalize by impulse peak - for v in &mut output[..n] { - *v *= inv_peak; + output[t] = g1 * output[t - 1] + g2 * output[t - 2] + source[t] * inv_peak; } } /// Adjoint convolution: normalized adjoint, O(T). /// - /// Adjoint of (K / peak) = K^T / peak. + /// Pre-scales input by 1/peak so the backward AR2 recursion directly + /// produces a peak-normalized output — no second normalization pass needed. pub(crate) fn convolve_adjoint(&self, source: &[f32], output: &mut [f32]) { let n = source.len(); if n == 0 { @@ -88,17 +80,12 @@ impl BandedAR2 { let g2 = self.g2 as f32; let inv_peak = (1.0 / self.impulse_peak) as f32; - output[n - 1] = source[n - 1]; + output[n - 1] = source[n - 1] * inv_peak; if n > 1 { - output[n - 2] = source[n - 2] + g1 * output[n - 1]; + output[n - 2] = source[n - 2] * inv_peak + g1 * output[n - 1]; } for t in (0..n.saturating_sub(2)).rev() { - output[t] = source[t] + g1 * output[t + 1] + g2 * output[t + 2]; - } - - // Normalize by impulse peak - for v in &mut output[..n] { - *v *= inv_peak; + output[t] = source[t] * inv_peak + g1 * output[t + 1] + g2 * output[t + 2]; } } diff --git a/crates/solver/src/baseline.rs b/crates/solver/src/baseline.rs index a61bd759..035344ed 100644 --- a/crates/solver/src/baseline.rs +++ b/crates/solver/src/baseline.rs @@ -45,12 +45,19 @@ impl Ord for OrderedF32 { /// Used for O(log M) k-th element queries via binary lifting. struct FenwickTree { tree: Vec, + msb: usize, // highest power of 2 <= (tree.len() - 1) } impl FenwickTree { fn new(size: usize) -> Self { + let mut msb = 1; + while msb <= size { + msb <<= 1; + } + msb >>= 1; Self { - tree: vec![0; size + 1], // 1-indexed + tree: vec![0; size + 1], + msb, } } @@ -66,14 +73,9 @@ impl FenwickTree { /// Find the 0-indexed position of the k-th element (1-based k). /// Uses binary lifting: O(log M) time. fn kth(&self, mut k: i32) -> usize { - let n = self.tree.len() - 1; // max 0-indexed position + 1 + let n = self.tree.len() - 1; let mut pos = 0; - // Find highest power of 2 <= n - let mut bit = 1; - while bit <= n { - bit <<= 1; - } - bit >>= 1; + let mut bit = self.msb; while bit > 0 { let next = pos + bit; @@ -83,7 +85,7 @@ impl FenwickTree { } bit >>= 1; } - pos // 0-indexed coordinate + pos } } diff --git a/crates/solver/src/biexp_fit.rs b/crates/solver/src/biexp_fit.rs index eec3e66c..1da14b19 100644 --- a/crates/solver/src/biexp_fit.rs +++ b/crates/solver/src/biexp_fit.rs @@ -82,6 +82,7 @@ /// the ceiling logic in eval_two_component can be simplified back to a /// plain `bs >= bf` gate. +#[derive(Clone)] #[cfg_attr(feature = "jsbindings", derive(serde::Serialize))] pub struct BiexpResult { pub tau_rise: f64, @@ -93,6 +94,25 @@ pub struct BiexpResult { pub beta_fast: f64, } +impl BiexpResult { + fn sentinel() -> Self { + BiexpResult { + tau_rise: 0.02, + tau_decay: 0.4, + beta: 0.0, + residual: f64::INFINITY, + tau_rise_fast: 0.0, + tau_decay_fast: 0.0, + beta_fast: 0.0, + } + } + + /// Returns true if the fit includes a fast component. + pub fn has_fast_component(&self) -> bool { + self.tau_rise_fast > 0.0 && self.tau_decay_fast > self.tau_rise_fast + } +} + /// Fit a two-component bi-exponential model to a free-form kernel. /// /// Uses a 20×20×(5×8+1) grid search over (tau_r, tau_d, tau_r_fast, tau_d_fast) @@ -116,15 +136,7 @@ pub fn fit_biexponential( let n = h_free.len(); let skip = skip.min(n.saturating_sub(1)); if n == 0 { - return BiexpResult { - tau_rise: 0.02, - tau_decay: 0.4, - beta: 0.0, - residual: f64::INFINITY, - tau_rise_fast: 0.0, - tau_decay_fast: 0.0, - beta_fast: 0.0, - }; + return BiexpResult::sentinel(); } let dt = 1.0 / fs; @@ -145,15 +157,7 @@ pub fn fit_biexponential( // additional candidate. This gives faster convergence when the kernel // is evolving smoothly between iterations. if let Some(warm) = warm_start { - let mut warm_candidate = BiexpResult { - tau_rise: warm.tau_rise, - tau_decay: warm.tau_decay, - beta: warm.beta, - residual: warm.residual, - tau_rise_fast: warm.tau_rise_fast, - tau_decay_fast: warm.tau_decay_fast, - beta_fast: warm.beta_fast, - }; + let mut warm_candidate = warm.clone(); // Re-evaluate on the CURRENT h_free (warm residual was from previous h_free) let (bs, bf, res) = eval_two_component( h_free, @@ -172,7 +176,7 @@ pub fn fit_biexponential( refine_candidate(h_free, &mut warm_candidate, dt, 40, skip); } // Warm candidate competes with the appropriate track - if warm_candidate.tau_rise_fast > 0.0 && warm_candidate.tau_decay_fast > warm_candidate.tau_rise_fast { + if warm_candidate.has_fast_component() { if warm_candidate.residual < best_two.residual { best_two = warm_candidate; } @@ -218,8 +222,15 @@ fn refine_candidate( ) { let (refined_tr, refined_td, refined_trf, refined_tdf) = golden_section_refine(h_free, candidate, dt, max_steps, skip); - let (beta_s, beta_f, residual) = - eval_two_component(h_free, refined_tr, refined_td, refined_trf, refined_tdf, dt, skip); + let (beta_s, beta_f, residual) = eval_two_component( + h_free, + refined_tr, + refined_td, + refined_trf, + refined_tdf, + dt, + skip, + ); if residual < candidate.residual { *candidate = BiexpResult { tau_rise: refined_tr, @@ -285,24 +296,8 @@ fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResu let tdf_lo = 0.5 * dt; let tdf_abs_hi = 8.0 * dt; - let mut best_slow = BiexpResult { - tau_rise: 0.02, - tau_decay: 0.4, - beta: 0.0, - residual: f64::INFINITY, - tau_rise_fast: 0.0, - tau_decay_fast: 0.0, - beta_fast: 0.0, - }; - let mut best_two = BiexpResult { - tau_rise: 0.02, - tau_decay: 0.4, - beta: 0.0, - residual: f64::INFINITY, - tau_rise_fast: 0.0, - tau_decay_fast: 0.0, - beta_fast: 0.0, - }; + let mut best_slow = BiexpResult::sentinel(); + let mut best_two = BiexpResult::sentinel(); for i in 0..grid_n { let log_tr = log_tr_lo + (log_tr_hi - log_tr_lo) * i as f64 / (grid_n - 1) as f64; @@ -334,7 +329,7 @@ fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResu // Inner grid: scan independent (tau_r_fast, tau_d_fast) // Upper bound for tau_d_fast is the tighter of the absolute cap - // (8×dt) and a relative cap (tau_d × 0.2) to prevent degeneracy. + // (8×dt) and a relative cap (tau_d × 0.15) to prevent degeneracy. let tdf_hi = tdf_abs_hi.min(tau_d * 0.15); if tdf_hi <= tdf_lo { continue; // tau_d too small for a distinct fast component @@ -343,8 +338,7 @@ fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResu let log_tdf_hi = tdf_hi.ln(); for ki in 0..trf_grid_n { - let tau_r_fast = - trf_lo + (trf_hi - trf_lo) * ki as f64 / (trf_grid_n - 1) as f64; + let tau_r_fast = trf_lo + (trf_hi - trf_lo) * ki as f64 / (trf_grid_n - 1) as f64; for kj in 0..tdf_grid_n { let log_tdf = log_tdf_lo @@ -356,15 +350,8 @@ fn cold_grid_search(h_free: &[f32], fs: f64, dt: f64, skip: usize) -> (BiexpResu continue; } - let (beta_s, beta_f, residual) = eval_two_component( - h_free, - tau_r, - tau_d, - tau_r_fast, - tau_d_fast, - dt, - skip, - ); + let (beta_s, beta_f, residual) = + eval_two_component(h_free, tau_r, tau_d, tau_r_fast, tau_d_fast, dt, skip); if residual < best_two.residual { best_two = BiexpResult { tau_rise: tau_r, @@ -527,6 +514,23 @@ fn eval_two_component( (best_bs, best_bf, best_res) } +/// Run one golden-section narrowing pass on a 1D interval [lo, hi]. +/// `cost` takes a candidate value and returns the residual. +/// Returns the midpoint of the narrowed interval. +fn golden_bracket(mut lo: f64, mut hi: f64, cost: impl Fn(f64) -> f64) -> f64 { + const PHI: f64 = 0.6180339887498949; // (sqrt(5) - 1) / 2 + for _ in 0..10 { + let x1 = hi - PHI * (hi - lo); + let x2 = lo + PHI * (hi - lo); + if cost(x1) < cost(x2) { + hi = x2; + } else { + lo = x1; + } + } + (lo + hi) / 2.0 +} + /// Golden-section refinement around the best grid point. /// Cycles through refining tau_r, tau_d, tau_r_fast, and tau_d_fast for `max_steps` total. fn golden_section_refine( @@ -536,109 +540,56 @@ fn golden_section_refine( max_steps: usize, skip: usize, ) -> (f64, f64, f64, f64) { - let phi = (5.0_f64.sqrt() - 1.0) / 2.0; // golden ratio conjugate - let mut tau_r = best.tau_rise; let mut tau_d = best.tau_decay; let mut tau_r_fast = best.tau_rise_fast; let mut tau_d_fast = best.tau_decay_fast; - // If fast component is zero, skip fast parameter refinement - let has_fast = tau_r_fast > 0.0 && tau_d_fast > tau_r_fast; + let has_fast = best.has_fast_component(); let n_phases = if has_fast { 4 } else { 2 }; for step in 0..max_steps { - let phase = step % n_phases; - - if phase == 0 { - // Refine tau_r - let mut lo = (tau_r * 0.5).max(dt); - let mut hi = (tau_r * 2.0).min(tau_d * 0.99); - if lo >= hi { - continue; - } - - for _ in 0..10 { - let x1 = hi - phi * (hi - lo); - let x2 = lo + phi * (hi - lo); - let (_, _, r1) = - eval_two_component(h_free, x1, tau_d, tau_r_fast, tau_d_fast, dt, skip); - let (_, _, r2) = - eval_two_component(h_free, x2, tau_d, tau_r_fast, tau_d_fast, dt, skip); - if r1 < r2 { - hi = x2; - } else { - lo = x1; + match step % n_phases { + 0 => { + // Refine tau_r + let lo = (tau_r * 0.5).max(dt); + let hi = (tau_r * 2.0).min(tau_d * 0.99); + if lo < hi { + tau_r = golden_bracket(lo, hi, |x| { + eval_two_component(h_free, x, tau_d, tau_r_fast, tau_d_fast, dt, skip).2 + }); } } - tau_r = (lo + hi) / 2.0; - } else if phase == 1 { - // Refine tau_d - let lo = (tau_d * 0.5).max(tau_r * 1.01); - let mut hi = tau_d * 2.0; - if lo >= hi { - continue; - } - - let mut lo = lo; - for _ in 0..10 { - let x1 = hi - phi * (hi - lo); - let x2 = lo + phi * (hi - lo); - let (_, _, r1) = - eval_two_component(h_free, tau_r, x1, tau_r_fast, tau_d_fast, dt, skip); - let (_, _, r2) = - eval_two_component(h_free, tau_r, x2, tau_r_fast, tau_d_fast, dt, skip); - if r1 < r2 { - hi = x2; - } else { - lo = x1; + 1 => { + // Refine tau_d + let lo = (tau_d * 0.5).max(tau_r * 1.01); + let hi = tau_d * 2.0; + if lo < hi { + tau_d = golden_bracket(lo, hi, |x| { + eval_two_component(h_free, tau_r, x, tau_r_fast, tau_d_fast, dt, skip).2 + }); } } - tau_d = (lo + hi) / 2.0; - } else if phase == 2 { - // Refine tau_r_fast - let mut lo = (tau_r_fast * 0.5).max(0.1 * dt); - let mut hi = (tau_r_fast * 2.0).min((2.0 * dt).min(tau_d_fast * 0.99)); - if lo >= hi { - continue; - } - - for _ in 0..10 { - let x1 = hi - phi * (hi - lo); - let x2 = lo + phi * (hi - lo); - let (_, _, r1) = - eval_two_component(h_free, tau_r, tau_d, x1, tau_d_fast, dt, skip); - let (_, _, r2) = - eval_two_component(h_free, tau_r, tau_d, x2, tau_d_fast, dt, skip); - if r1 < r2 { - hi = x2; - } else { - lo = x1; + 2 => { + // Refine tau_r_fast + let lo = (tau_r_fast * 0.5).max(0.1 * dt); + let hi = (tau_r_fast * 2.0).min((2.0 * dt).min(tau_d_fast * 0.99)); + if lo < hi { + tau_r_fast = golden_bracket(lo, hi, |x| { + eval_two_component(h_free, tau_r, tau_d, x, tau_d_fast, dt, skip).2 + }); } } - tau_r_fast = (lo + hi) / 2.0; - } else { - // Refine tau_d_fast — cap at min(8×dt, tau_d × 0.15) to prevent degeneracy - let mut lo = (tau_d_fast * 0.5).max(tau_r_fast * 1.01); - let mut hi = (tau_d_fast * 2.0).min((8.0 * dt).min(tau_d * 0.15)); - if lo >= hi { - continue; - } - - for _ in 0..10 { - let x1 = hi - phi * (hi - lo); - let x2 = lo + phi * (hi - lo); - let (_, _, r1) = - eval_two_component(h_free, tau_r, tau_d, tau_r_fast, x1, dt, skip); - let (_, _, r2) = - eval_two_component(h_free, tau_r, tau_d, tau_r_fast, x2, dt, skip); - if r1 < r2 { - hi = x2; - } else { - lo = x1; + _ => { + // Refine tau_d_fast + let lo = (tau_d_fast * 0.5).max(tau_r_fast * 1.01); + let hi = (tau_d_fast * 2.0).min((8.0 * dt).min(tau_d * 0.15)); + if lo < hi { + tau_d_fast = golden_bracket(lo, hi, |x| { + eval_two_component(h_free, tau_r, tau_d, tau_r_fast, x, dt, skip).2 + }); } } - tau_d_fast = (lo + hi) / 2.0; } } @@ -921,11 +872,7 @@ mod tests { #[test] fn fast_tau_in_valid_range() { // For various inputs, verify fast component time constants are in expected ranges - let test_cases = [ - (0.08, 0.5, 30.0), - (0.05, 0.3, 100.0), - (0.1, 2.0, 10.0), - ]; + let test_cases = [(0.08, 0.5, 30.0), (0.05, 0.3, 100.0), (0.1, 2.0, 10.0)]; for (tau_r, tau_d, fs) in test_cases { let h = make_biexp(tau_r, tau_d, 2.0, fs, 60); @@ -966,8 +913,7 @@ mod tests { // Case 1: Pure slow component — should yield beta_s > 0, beta_f ≈ 0 let h_slow = make_biexp(0.05, 0.5, 2.0, fs, n); - let (bs, bf, _) = - eval_two_component(&h_slow, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); + let (bs, bf, _) = eval_two_component(&h_slow, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); assert!(bs > 0.0, "beta_s should be positive for slow-only input"); assert!( bf < 0.1 * bs, @@ -985,22 +931,19 @@ mod tests { (3.0 * ((-t / tau_d_fast).exp() - (-t / tau_r_fast).exp())) as f32 }) .collect(); - let (bs, _bf, _) = - eval_two_component(&h_fast, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); + let (bs, _bf, _) = eval_two_component(&h_fast, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); // With no fast-only active set, the slow template absorbs what it can assert!(bs >= 0.0, "beta_s should be non-negative for any input"); // Case 3: Both components present let h_both = make_two_component(0.05, 0.5, 2.0, tau_r_fast, tau_d_fast, 1.5, fs, n); - let (bs, bf, _) = - eval_two_component(&h_both, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); + let (bs, bf, _) = eval_two_component(&h_both, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); assert!(bs > 0.0, "beta_s should be positive for mixed input"); assert!(bf > 0.0, "beta_f should be positive for mixed input"); // Case 4: Zero signal — both should be zero let h_zero = vec![0.0_f32; n]; - let (bs, bf, res) = - eval_two_component(&h_zero, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); + let (bs, bf, res) = eval_two_component(&h_zero, 0.05, 0.5, tau_r_fast, tau_d_fast, dt, 0); assert_eq!(bs, 0.0, "beta_s should be zero for zero input"); assert_eq!(bf, 0.0, "beta_f should be zero for zero input"); assert!(res < 1e-20, "residual should be ~0 for zero input"); diff --git a/crates/solver/src/fft.rs b/crates/solver/src/fft.rs index 2525f0c1..25513eb6 100644 --- a/crates/solver/src/fft.rs +++ b/crates/solver/src/fft.rs @@ -150,41 +150,7 @@ impl FftConvolver { signal_len: usize, output: &mut [f32], ) { - let padded_len = self.fft_len; - let spectrum_len = padded_len / 2 + 1; - - // Zero-pad source into fft_input - self.fft_input[..signal_len].copy_from_slice(&source[..signal_len]); - self.fft_input[signal_len..padded_len].fill(0.0); - - // Forward FFT of source - let fwd = self.plan_fwd.as_ref().expect("plans not initialized"); - fwd.process_with_scratch( - &mut self.fft_input[..padded_len], - &mut self.fft_spectrum[..spectrum_len], - &mut self.fft_scratch_fwd, - ) - .unwrap(); - - // Pointwise multiply with kernel FFT - for i in 0..spectrum_len { - self.fft_spectrum[i] *= self.kernel_fft[i]; - } - - // Inverse FFT - let inv = self.plan_inv.as_ref().expect("plans not initialized"); - inv.process_with_scratch( - &mut self.fft_spectrum[..spectrum_len], - &mut self.fft_output[..padded_len], - &mut self.fft_scratch_inv, - ) - .unwrap(); - - // Normalize and copy first signal_len samples to output - let scale = 1.0 / padded_len as f32; - for i in 0..signal_len { - output[i] = self.fft_output[i] * scale; - } + self.convolve_impl(source, signal_len, output, false); } /// FFT-based adjoint convolution (correlation): output[..signal_len] = (K^T * source)[..signal_len]. @@ -193,6 +159,18 @@ impl FftConvolver { source: &[f32], signal_len: usize, output: &mut [f32], + ) { + self.convolve_impl(source, signal_len, output, true); + } + + /// Shared FFT convolution implementation. + /// `use_conjugate` selects `kernel_conj_fft` (adjoint) or `kernel_fft` (forward). + fn convolve_impl( + &mut self, + source: &[f32], + signal_len: usize, + output: &mut [f32], + use_conjugate: bool, ) { let padded_len = self.fft_len; let spectrum_len = padded_len / 2 + 1; @@ -210,9 +188,15 @@ impl FftConvolver { ) .unwrap(); - // Pointwise multiply with conjugate kernel FFT - for i in 0..spectrum_len { - self.fft_spectrum[i] *= self.kernel_conj_fft[i]; + // Pointwise multiply with kernel spectrum + if use_conjugate { + for i in 0..spectrum_len { + self.fft_spectrum[i] *= self.kernel_conj_fft[i]; + } + } else { + for i in 0..spectrum_len { + self.fft_spectrum[i] *= self.kernel_fft[i]; + } } // Inverse FFT diff --git a/crates/solver/src/filter.rs b/crates/solver/src/filter.rs index 52d0fcf5..1173a4ac 100644 --- a/crates/solver/src/filter.rs +++ b/crates/solver/src/filter.rs @@ -206,17 +206,8 @@ impl BandpassFilter { } } - /// Apply bandpass filter in-place. Caches power spectrum. Returns false if skipped. - pub fn apply(&mut self, trace: &mut [f32]) -> bool { - if !self.is_enabled() || !self.valid || trace.len() < 8 { - return false; - } - - // Mode-specific validity: HP+LP requires f_hp < f_lp - if self.hp_enabled && self.lp_enabled && self.f_hp >= self.f_lp { - return false; - } - + /// Perform forward FFT and cache power spectrum. Used by both `apply` and `compute_spectrum_only`. + fn forward_fft_and_cache_power(&mut self, trace: &[f32]) { let n = trace.len(); self.ensure_buffers(n); let spectrum_len = n / 2 + 1; @@ -224,7 +215,7 @@ impl BandpassFilter { // Copy trace into fft_input self.fft_input[..n].copy_from_slice(trace); - // Forward FFT (use cached plan — no hash-map lookup) + // Forward FFT let fwd = self.plan_fwd.as_ref().expect("plans not initialized"); fwd.process_with_scratch( &mut self.fft_input[..n], @@ -234,14 +225,35 @@ impl BandpassFilter { .unwrap(); // Cache pre-filter power spectrum - for i in 0..spectrum_len { - let c = self.spectrum[i]; - self.power_spectrum[i] = c.re * c.re + c.im * c.im; + for (ps, c) in self.power_spectrum[..spectrum_len] + .iter_mut() + .zip(&self.spectrum[..spectrum_len]) + { + *ps = c.re * c.re + c.im * c.im; } + } + + /// Apply bandpass filter in-place. Caches power spectrum. Returns false if skipped. + pub fn apply(&mut self, trace: &mut [f32]) -> bool { + if !self.is_enabled() || !self.valid || trace.len() < 8 { + return false; + } + + // Mode-specific validity: HP+LP requires f_hp < f_lp + if self.hp_enabled && self.lp_enabled && self.f_hp >= self.f_lp { + return false; + } + + let n = trace.len(); + self.forward_fft_and_cache_power(trace); + let spectrum_len = n / 2 + 1; // Apply gain curve - for i in 0..spectrum_len { - self.spectrum[i] *= self.gain_curve[i]; + for (s, &g) in self.spectrum[..spectrum_len] + .iter_mut() + .zip(&self.gain_curve[..spectrum_len]) + { + *s *= g; } // Inverse FFT (use cached plan — no hash-map lookup) @@ -255,8 +267,8 @@ impl BandpassFilter { // Normalize (realfft doesn't normalize) let scale = 1.0 / n as f32; - for i in 0..n { - trace[i] = self.fft_input[i] * scale; + for (t, &f) in trace.iter_mut().zip(&self.fft_input[..n]) { + *t = f * scale; } true @@ -267,25 +279,7 @@ impl BandpassFilter { if trace.len() < 8 { return; } - - let n = trace.len(); - self.ensure_buffers(n); - let spectrum_len = n / 2 + 1; - - self.fft_input[..n].copy_from_slice(trace); - - let fwd = self.plan_fwd.as_ref().expect("plans not initialized"); - fwd.process_with_scratch( - &mut self.fft_input[..n], - &mut self.spectrum[..spectrum_len], - &mut self.scratch_fwd, - ) - .unwrap(); - - for i in 0..spectrum_len { - let c = self.spectrum[i]; - self.power_spectrum[i] = c.re * c.re + c.im * c.im; - } + self.forward_fft_and_cache_power(trace); } /// Get power spectrum (N/2+1 bins of |FFT|²). diff --git a/crates/solver/src/fista.rs b/crates/solver/src/fista.rs index 185b1a65..6669e024 100644 --- a/crates/solver/src/fista.rs +++ b/crates/solver/src/fista.rs @@ -55,20 +55,9 @@ impl Solver { // mathematically cancels in the gradient (residual = mean-centered signals). // Computing it anyway would produce pure momentum-oscillation noise. if !self.filtered { - let mut sum = 0.0_f64; - for i in 0..n { - sum += (self.trace[i] - self.reconvolution[i]) as f64; - } - let raw_baseline = sum / n as f64; - self.baseline = raw_baseline; - - // Per-iteration EMA smoothing for display baseline - if !self.baseline_ema_init { - self.baseline_ema = raw_baseline; - self.baseline_ema_init = true; - } else { - self.baseline_ema = 0.3 * raw_baseline + 0.7 * self.baseline_ema; - } + let raw = + crate::compute_raw_baseline(&self.trace[..n], &self.reconvolution[..n], n); + self.update_baseline_ema(raw); } // 2. Compute residual = K * y_k + b - trace diff --git a/crates/solver/src/indeca.rs b/crates/solver/src/indeca.rs index fc857f27..503c0620 100644 --- a/crates/solver/src/indeca.rs +++ b/crates/solver/src/indeca.rs @@ -138,6 +138,19 @@ fn solve_upsampled( (solution, filtered, iterations, converged) } +/// Return the interior slice of `s` excluding `pad` samples from each end. +/// Falls back to the full slice when the interior is empty. +fn interior_slice(s: &[f32], pad: usize) -> &[f32] { + let n = s.len(); + let lo = pad.min(n); + let hi = n.saturating_sub(pad).max(lo); + if hi > lo { + &s[lo..hi] + } else { + s + } +} + /// Estimate alpha from the interior of the trace (excluding boundary padding). /// /// Uses peak-to-trough of the inner region to avoid edge artifacts that occur @@ -145,17 +158,7 @@ fn solve_upsampled( /// Since the kernel is peak-normalized, peak-to-trough >= alpha, making this /// a safe overestimate. Returns 1.0 for flat traces. fn estimate_alpha_interior(trace: &[f32], pad: usize) -> f64 { - let n = trace.len(); - let lo_idx = pad.min(n); - let hi_idx = n.saturating_sub(pad).max(lo_idx); - let inner = &trace[lo_idx..hi_idx]; - if inner.is_empty() { - // Trace too short for padding — fall back to full trace - let lo = trace.iter().copied().fold(f32::INFINITY, f32::min); - let hi = trace.iter().copied().fold(f32::NEG_INFINITY, f32::max); - let ptp = (hi - lo) as f64; - return if ptp < 1e-10 { 1.0 } else { ptp }; - } + let inner = interior_slice(trace, pad); let lo = inner.iter().copied().fold(f32::INFINITY, f32::min); let hi = inner.iter().copied().fold(f32::NEG_INFINITY, f32::max); let ptp = (hi - lo) as f64; @@ -170,11 +173,10 @@ fn estimate_alpha_interior(trace: &[f32], pad: usize) -> f64 { /// /// Falls back to the full slice when the interior is empty. fn interior_peak(s: &[f32], pad: usize) -> f32 { - let n = s.len(); - let lo = pad.min(n); - let hi = n.saturating_sub(pad).max(lo); - let region = if hi > lo { &s[lo..hi] } else { s }; - region.iter().copied().fold(0.0_f32, f32::max) + interior_slice(s, pad) + .iter() + .copied() + .fold(0.0_f32, f32::max) } /// Full InDeCa trace processing pipeline with scale iteration. @@ -214,25 +216,17 @@ pub fn solve_trace( let mut solver = Solver::new(); // ── Step 1: Apply optional bandpass filter + rolling baseline subtraction ── - // Run a throwaway FISTA just to get the filtered trace (if HP/LP), then + // Apply bandpass filter directly (if HP/LP enabled), then // subtract the rolling-percentile baseline so the floor is ~0. let mut working_trace = if hp_enabled || lp_enabled { - let (_, filtered_up, _, _) = solve_upsampled( - &mut solver, - &upsampled, - tau_r, - tau_d, - fs_up, - 1, // only 1 iteration — we just need the filtered trace - tol, - None, - hp_enabled, - lp_enabled, - Constraint::Box01, - false, - 0.0, // no sparsity for filter pass - ); - filtered_up.unwrap() + // Apply bandpass filter directly — no need for a full FISTA solve + solver.set_conv_mode(ConvMode::BandedAR2); + solver.set_params(tau_r, tau_d, 0.0, fs_up); + solver.set_trace(&upsampled); + solver.set_hp_filter_enabled(hp_enabled); + solver.set_lp_filter_enabled(lp_enabled); + solver.apply_filter(); + solver.get_trace() } else { upsampled }; @@ -412,7 +406,9 @@ mod tests { #[test] fn outputs_in_range() { let trace = make_trace(0.02, 0.4, 30.0, 300, &[20, 80, 150, 220]); - let result = solve_trace(&trace, 0.02, 0.4, 30.0, 1, 500, 1e-4, None, false, false, 0.0); + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 1, 500, 1e-4, None, false, false, 0.0, + ); // Spike counts should be non-negative for (i, &v) in result.s_counts.iter().enumerate() { @@ -438,7 +434,9 @@ mod tests { } } } - let result = solve_trace(&trace, 0.02, 0.4, 30.0, 1, 1000, 1e-4, None, false, false, 0.0); + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 1, 1000, 1e-4, None, false, false, 0.0, + ); // Check that spikes are detected near the true positions let mut detected = 0; @@ -491,7 +489,9 @@ mod tests { #[test] fn upsampled_output_length() { let trace = make_trace(0.02, 0.4, 30.0, 100, &[20, 50]); - let result = solve_trace(&trace, 0.02, 0.4, 30.0, 10, 200, 1e-3, None, false, false, 0.0); + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 10, 200, 1e-3, None, false, false, 0.0, + ); // Output should be same length as input regardless of upsample factor assert_eq!( @@ -504,7 +504,9 @@ mod tests { #[test] fn zero_trace() { let trace = vec![0.0_f32; 100]; - let result = solve_trace(&trace, 0.02, 0.4, 30.0, 1, 100, 1e-4, None, false, false, 0.0); + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 1, 100, 1e-4, None, false, false, 0.0, + ); let total_spikes: f32 = result.s_counts.iter().sum(); assert!( total_spikes < 1e-6, @@ -538,7 +540,9 @@ mod tests { } } - let result = solve_trace(&trace, tau_r, tau_d, fs, 10, 500, 1e-4, None, false, false, 0.0); + let result = solve_trace( + &trace, tau_r, tau_d, fs, 10, 500, 1e-4, None, false, false, 0.0, + ); let total_counts: f32 = result.s_counts.iter().sum(); @@ -602,7 +606,9 @@ mod tests { let subset_end = 400; let subset = &full_trace[subset_start..subset_end]; - let result = solve_trace(subset, tau_r, tau_d, fs, 1, 1000, 1e-4, None, false, false, 0.0); + let result = solve_trace( + subset, tau_r, tau_d, fs, 1, 1000, 1e-4, None, false, false, 0.0, + ); let total_spikes: f32 = result.s_counts.iter().sum(); // Should detect interior spikes, not just the edge artifact @@ -638,7 +644,9 @@ mod tests { } } - let result = solve_trace(&trace, tau_r, tau_d, fs, 1, 1000, 1e-4, None, false, false, 0.0); + let result = solve_trace( + &trace, tau_r, tau_d, fs, 1, 1000, 1e-4, None, false, false, 0.0, + ); let total_spikes: f32 = result.s_counts.iter().sum(); assert!( @@ -647,4 +655,120 @@ mod tests { total_spikes, result.alpha, result.threshold, result.pve ); } + + /// HP+LP filter path should produce valid results and return a filtered trace. + #[test] + fn filter_path_hp_lp() { + let spike_positions = [30, 100, 200]; + let alpha_true = 10.0_f32; + let baseline_true = 2.0_f32; + let kernel = build_kernel(0.02, 0.4, 30.0); + let n = 300; + let mut trace = vec![baseline_true; n]; + for &pos in &spike_positions { + for (k, &kv) in kernel.iter().enumerate() { + if pos + k < n { + trace[pos + k] += alpha_true * kv; + } + } + } + + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 1, 1000, 1e-4, None, true, true, 0.0, + ); + + // Output length should match input + assert_eq!(result.s_counts.len(), trace.len()); + + // Spike counts should be non-negative + for (i, &v) in result.s_counts.iter().enumerate() { + assert!(v >= 0.0, "Negative spike count at {}: {}", i, v); + } + + // Filtered trace should be returned and have the correct length + let filtered = result + .filtered_trace + .as_ref() + .expect("filtered_trace should be Some when filters are enabled"); + assert_eq!(filtered.len(), trace.len()); + + // Should still detect spikes through the filter + let total_spikes: f32 = result.s_counts.iter().sum(); + assert!( + total_spikes >= 1.0, + "Should detect at least 1 spike with HP+LP filter, got {} (pve={:.4})", + total_spikes, + result.pve + ); + } + + /// HP-only filter path should remove DC and still detect spikes. + #[test] + fn filter_path_hp_only() { + let spike_positions = [30, 100, 200]; + let alpha_true = 10.0_f32; + let baseline_true = 50.0_f32; // high DC offset + let kernel = build_kernel(0.02, 0.4, 30.0); + let n = 300; + let mut trace = vec![baseline_true; n]; + for &pos in &spike_positions { + for (k, &kv) in kernel.iter().enumerate() { + if pos + k < n { + trace[pos + k] += alpha_true * kv; + } + } + } + + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 1, 1000, 1e-4, None, true, false, 0.0, + ); + + assert_eq!(result.s_counts.len(), trace.len()); + + // Filtered trace should be returned + assert!(result.filtered_trace.is_some()); + + // Should still detect spikes + let total_spikes: f32 = result.s_counts.iter().sum(); + assert!( + total_spikes >= 1.0, + "Should detect at least 1 spike with HP-only filter, got {} (pve={:.4})", + total_spikes, + result.pve + ); + } + + /// LP-only filter path should preserve DC and detect spikes. + #[test] + fn filter_path_lp_only() { + let spike_positions = [30, 100, 200]; + let alpha_true = 10.0_f32; + let baseline_true = 2.0_f32; + let kernel = build_kernel(0.02, 0.4, 30.0); + let n = 300; + let mut trace = vec![baseline_true; n]; + for &pos in &spike_positions { + for (k, &kv) in kernel.iter().enumerate() { + if pos + k < n { + trace[pos + k] += alpha_true * kv; + } + } + } + + let result = solve_trace( + &trace, 0.02, 0.4, 30.0, 1, 1000, 1e-4, None, false, true, 0.0, + ); + + assert_eq!(result.s_counts.len(), trace.len()); + assert!(result.filtered_trace.is_some()); + + // Should still detect spikes + let total_spikes: f32 = result.s_counts.iter().sum(); + assert!( + total_spikes >= 1.0, + "Should detect at least 1 spike with LP-only filter, got {} (pve={:.4})", + total_spikes, + result.pve + ); + } } diff --git a/crates/solver/src/js_indeca.rs b/crates/solver/src/js_indeca.rs index b33d0cb9..9f25bdd3 100644 --- a/crates/solver/src/js_indeca.rs +++ b/crates/solver/src/js_indeca.rs @@ -128,8 +128,7 @@ pub fn indeca_fit_biexponential( } else { None }; - let result = - biexp_fit::fit_biexponential(h_free, fs, refine, skip, warm_start.as_ref()); + let result = biexp_fit::fit_biexponential(h_free, fs, refine, skip, warm_start.as_ref()); serde_wasm_bindgen::to_value(&result).unwrap_or(JsValue::NULL) } diff --git a/crates/solver/src/kernel.rs b/crates/solver/src/kernel.rs index 1561bbb2..36675085 100644 --- a/crates/solver/src/kernel.rs +++ b/crates/solver/src/kernel.rs @@ -1,15 +1,20 @@ +/// Clamp tau_rise away from tau_decay to prevent degenerate zero kernels. +/// When tau_rise ≈ tau_decay, the biexponential exp(-t/τ_d) - exp(-t/τ_r) collapses to zero. +pub(crate) fn clamp_tau_rise(tau_rise: f64, tau_decay: f64) -> f64 { + if (tau_rise - tau_decay).abs() < 1e-6 * tau_decay.max(tau_rise).max(1e-12) { + tau_decay * 0.5 + } else { + tau_rise + } +} + /// Build a double-exponential calcium kernel normalized to peak = 1.0. /// /// h(t) = exp(-t/tau_decay) - exp(-t/tau_rise), normalized so max(h) = 1.0. /// Kernel length extends until the decay envelope drops below 1e-6 of peak. /// Computed in f64 for precision, returned as Vec. pub fn build_kernel(tau_rise: f64, tau_decay: f64, fs: f64) -> Vec { - // Guard: tau_rise too close to tau_decay produces a degenerate zero kernel - let tau_rise = if (tau_rise - tau_decay).abs() < 1e-6 * tau_decay.max(tau_rise).max(1e-12) { - tau_decay * 0.5 - } else { - tau_rise - }; + let tau_rise = clamp_tau_rise(tau_rise, tau_decay); let dt = 1.0 / fs; @@ -49,12 +54,7 @@ pub fn build_kernel(tau_rise: f64, tau_decay: f64, fs: f64) -> Vec { /// Used by BandedAR2 tests and the TypeScript port in src/lib/ar2.ts. #[allow(dead_code)] pub fn tau_to_ar2(tau_rise: f64, tau_decay: f64, fs: f64) -> (f64, f64) { - // Guard: tau_rise too close to tau_decay produces a degenerate zero kernel - let tau_rise = if (tau_rise - tau_decay).abs() < 1e-6 * tau_decay.max(tau_rise).max(1e-12) { - tau_decay * 0.5 - } else { - tau_rise - }; + let tau_rise = clamp_tau_rise(tau_rise, tau_decay); let dt = 1.0 / fs; let d = (-dt / tau_decay).exp(); // decay eigenvalue @@ -84,16 +84,19 @@ pub fn compute_lipschitz(kernel: &[f32]) -> f64 { let inv = 2.0 * std::f64::consts::PI / (fft_len as f64); let mut max_power = 0.0_f64; - for w in 0..fft_len { + // Pre-cast kernel to f64 once instead of per-frequency + let kernel_f64: Vec = kernel.iter().map(|&k| k as f64).collect(); + + // Real kernel has symmetric spectrum: only need 0..=fft_len/2 + for w in 0..=fft_len / 2 { let freq = inv * (w as f64); let mut re = 0.0_f64; let mut im = 0.0_f64; - for (k, &hk) in kernel.iter().enumerate() { - let hk64 = hk as f64; + for (k, &hk) in kernel_f64.iter().enumerate() { let angle = freq * (k as f64); let (s, c) = angle.sin_cos(); - re += hk64 * c; - im -= hk64 * s; + re += hk * c; + im -= hk * s; } let power = re * re + im * im; if power > max_power { diff --git a/crates/solver/src/kernel_est.rs b/crates/solver/src/kernel_est.rs index 89e3943c..6aa3a98e 100644 --- a/crates/solver/src/kernel_est.rs +++ b/crates/solver/src/kernel_est.rs @@ -169,25 +169,18 @@ pub fn estimate_free_kernel( let mut stv = vec![0.0_f64; kernel_length]; // S^T S v let mut eigenvalue = 1.0_f64; + let mut v_f32 = vec![0.0_f32; kernel_length]; + for _ in 0..20 { - // S*v: convolve spikes with v (cast to f32) - let v_f32: Vec = v.iter().map(|&x| x as f32).collect(); + // Cast v to f32 into pre-allocated buffer + for (dst, &src) in v_f32.iter_mut().zip(v.iter()) { + *dst = src as f32; + } + // S*v: convolve spikes with v convolve_spikes_kernel(spike_trains, trace_lengths, &v_f32, &mut sv); // S^T (S*v) - stv.fill(0.0); - let mut off = 0; - for i in 0..n_traces { - let len = trace_lengths[i]; - for t in 0..len { - let val = sv[off + t] as f64; - let k_max = kernel_length.min(t + 1); - for k in 0..k_max { - stv[k] += val * spike_trains[off + t - k] as f64; - } - } - off += len; - } + adjoint_spikes_kernel(&sv, spike_trains, trace_lengths, kernel_length, &mut stv); // eigenvalue estimate = ||S^T S v|| eigenvalue = stv.iter().map(|&x| x * x).sum::().sqrt(); @@ -219,30 +212,28 @@ pub fn estimate_free_kernel( // Working buffer for S*h (convolution result) let mut sh = vec![0.0_f32; total_len]; + let mut z = vec![0.0_f64; kernel_length]; + for iter in 0..max_iters { // Forward: S*h (convolve each trace's spikes with h) convolve_spikes_kernel(spike_trains, trace_lengths, &h_prev, &mut sh); - // Residual: r = S*h - y_adj - // Gradient: S^T * r - gradient.fill(0.0); - offset = 0; - for i in 0..n_traces { - let len = trace_lengths[i]; - for t in 0..len { - let r = sh[offset + t] as f64 - y_adj[offset + t] as f64; - // S^T contribution: h[k] gets r * s[t-k] - let k_max = kernel_length.min(t + 1); - for k in 0..k_max { - gradient[k] += r * spike_trains[offset + t - k] as f64; - } - } - offset += len; + // Residual: r = S*h - y_adj (compute in-place in sh) + for i in 0..total_len { + sh[i] -= y_adj[i]; } + // Gradient: S^T * r + adjoint_spikes_kernel( + &sh, + spike_trains, + trace_lengths, + kernel_length, + &mut gradient, + ); + // Proximal gradient step: gradient descent on data-fidelity, then // TV proximal operator, then non-negativity projection. - let mut z = vec![0.0_f64; kernel_length]; for k in 0..kernel_length { z[k] = h_prev[k] as f64 - step_size * gradient[k]; } @@ -284,6 +275,29 @@ pub fn estimate_free_kernel( h } +/// Adjoint of spike convolution: output[k] += sum_t input[t] * s[t-k]. +/// This is S^T * input, the transpose of convolve_spikes_kernel. +fn adjoint_spikes_kernel( + input: &[f32], + spikes: &[f32], + trace_lengths: &[usize], + kernel_length: usize, + output: &mut [f64], +) { + output[..kernel_length].fill(0.0); + let mut offset = 0; + for &len in trace_lengths { + for t in 0..len { + let val = input[offset + t] as f64; + let k_max = kernel_length.min(t + 1); + for k in 0..k_max { + output[k] += val * spikes[offset + t - k] as f64; + } + } + offset += len; + } +} + /// Convolve spike trains with kernel h: output[t] = sum_k h[k] * s[t-k]. fn convolve_spikes_kernel(spikes: &[f32], trace_lengths: &[usize], h: &[f32], output: &mut [f32]) { let k_len = h.len(); diff --git a/crates/solver/src/lib.rs b/crates/solver/src/lib.rs index 57f7dd9c..3d27d24d 100644 --- a/crates/solver/src/lib.rs +++ b/crates/solver/src/lib.rs @@ -162,8 +162,15 @@ impl Solver { self.kernel_dc_gain = self.kernel.iter().map(|&k| k as f64).sum(); self.bandpass.update_cutoffs(tau_rise, tau_decay, fs); - // Update both convolution engines - self.banded.update(tau_rise, tau_decay, fs); + // Update convolution engines (only the active one + compute Lipschitz) + match self.conv_mode { + ConvMode::BandedAR2 => { + self.banded.update(tau_rise, tau_decay, fs); + } + ConvMode::Fft => { + // banded will be updated lazily if conv_mode switches + } + } self.lipschitz_constant = self.current_lipschitz(); // Update kernel FFT if buffers are already set up and large enough. @@ -304,11 +311,19 @@ impl Solver { /// Does NOT reset solution/iteration state — warm-start is preserved. pub fn set_conv_mode(&mut self, mode: ConvMode) { self.conv_mode = mode; - self.lipschitz_constant = self.current_lipschitz(); - // Ensure FFT buffers exist if switching to FFT mode with an active trace - if mode == ConvMode::Fft && self.active_len > 0 { - self.fft.ensure_buffers(self.active_len, &self.kernel); + match mode { + ConvMode::BandedAR2 => { + // Ensure banded coefficients are current (may have been skipped in set_params) + self.banded.update(self.tau_rise, self.tau_decay, self.fs); + } + ConvMode::Fft => { + // Ensure FFT buffers exist if switching to FFT mode with an active trace + if self.active_len > 0 { + self.fft.ensure_buffers(self.active_len, &self.kernel); + } + } } + self.lipschitz_constant = self.current_lipschitz(); } /// Set the constraint type (NonNegative or Box01). @@ -333,8 +348,7 @@ impl Solver { /// Format: [active_len (u32)] [t_fista (f64)] [iteration (u32)] [baseline (f64)] [solution f32...] [solution_prev f32...] pub fn export_state(&self) -> Vec { let n = self.active_len; - // 4 bytes active_len + 8 bytes t_fista + 4 bytes iteration + 8 bytes baseline + 2*n*4 bytes solutions (f32) - let mut buf = Vec::with_capacity(4 + 8 + 4 + 8 + 2 * n * 4); + let mut buf = Vec::with_capacity(state_byte_len(n)); buf.extend_from_slice(&(n as u32).to_le_bytes()); buf.extend_from_slice(&self.t_fista.to_le_bytes()); @@ -386,26 +400,24 @@ impl Solver { // Recompute baseline at current solution for display alignment. // In step_batch, baseline is skipped when filtered (cancels in gradient), // but the display path always needs it to align fit with trace. - { - let mut sum = 0.0_f64; - for i in 0..n { - sum += (self.trace[i] - self.reconvolution[i]) as f64; - } - let raw_baseline = sum / n as f64; - self.baseline = raw_baseline; - - // EMA smoothing for display (damps momentum-induced oscillation) - if !self.baseline_ema_init { - self.baseline_ema = raw_baseline; - self.baseline_ema_init = true; - } else { - self.baseline_ema = 0.3 * raw_baseline + 0.7 * self.baseline_ema; - } - } + let raw = compute_raw_baseline(&self.trace[..n], &self.reconvolution[..n], n); + self.update_baseline_ema(raw); self.reconvolution_stale = false; } + /// Update the baseline EMA from a raw baseline estimate. + /// Called by both `step_batch` (per-iteration) and `compute_reconvolution` (lazy display path). + fn update_baseline_ema(&mut self, raw_baseline: f64) { + self.baseline = raw_baseline; + if !self.baseline_ema_init { + self.baseline_ema = raw_baseline; + self.baseline_ema_init = true; + } else { + self.baseline_ema = 0.3 * raw_baseline + 0.7 * self.baseline_ema; + } + } + // --- Bandpass filter methods --- /// Convenience: set both HP and LP together (used by CaTune's single toggle). @@ -495,7 +507,7 @@ impl Solver { let mut cur = Cursor::new(state); let saved_len = read_u32_le(&mut cur) as usize; - let expected_size = 4 + 8 + 4 + 8 + 2 * saved_len * 4; + let expected_size = state_byte_len(saved_len); if state.len() != expected_size || saved_len != self.active_len { return; // size mismatch, cold start @@ -516,6 +528,20 @@ impl Solver { } } +/// Compute the mean residual (trace - reconvolution) as the raw baseline estimate. +pub(crate) fn compute_raw_baseline(trace: &[f32], reconvolution: &[f32], n: usize) -> f64 { + let mut sum = 0.0_f64; + for i in 0..n { + sum += (trace[i] - reconvolution[i]) as f64; + } + sum / n as f64 +} + +/// Byte length of serialized solver state for a trace of length `n`. +fn state_byte_len(n: usize) -> usize { + 4 + 8 + 4 + 8 + 2 * n * 4 // u32 + f64 + u32 + f64 + 2×n×f32 +} + // --- Little-endian cursor read helpers --- // These wrap the repetitive read_exact + from_le_bytes pattern used by load_state. // Each panics on short reads, which cannot occur when the caller has already diff --git a/crates/solver/src/peak_seed.rs b/crates/solver/src/peak_seed.rs index 8317eb7b..583911da 100644 --- a/crates/solver/src/peak_seed.rs +++ b/crates/solver/src/peak_seed.rs @@ -8,7 +8,6 @@ /// 5. Feed into estimate_free_kernel() → fit_biexponential() (already exist) /// /// The result provides initial tau_rise, tau_decay for the normal iterative pipeline. - use crate::biexp_fit::{fit_biexponential, BiexpResult}; use crate::kernel_est::estimate_free_kernel; @@ -28,7 +27,7 @@ pub struct SeedTraceResult { /// into the kernel estimation step as a replacement for FISTA trace inference. pub fn seed_trace(trace: &[f32], fs: f64) -> SeedTraceResult { let n = trace.len(); - let bl = median(trace); + let (bl, _) = median_and_mad(trace); let onsets = find_seed_spikes(trace, fs, 5.0); let mut s_counts = vec![0.0_f32; n]; @@ -58,7 +57,8 @@ pub struct SeedKernelResult { } /// Median of a slice (copies + sorts). Returns 0.0 for empty input. -pub(crate) fn median(data: &[f32]) -> f32 { +#[cfg(test)] +fn median(data: &[f32]) -> f32 { if data.is_empty() { return 0.0; } @@ -73,7 +73,8 @@ pub(crate) fn median(data: &[f32]) -> f32 { } /// Median absolute deviation of a slice. -pub(crate) fn mad(data: &[f32], median_val: f32) -> f32 { +#[cfg(test)] +fn mad(data: &[f32], median_val: f32) -> f32 { if data.is_empty() { return 0.0; } @@ -81,6 +82,32 @@ pub(crate) fn mad(data: &[f32], median_val: f32) -> f32 { median(&deviations) } +/// Compute median and MAD in one sort pass. +pub(crate) fn median_and_mad(data: &[f32]) -> (f32, f32) { + if data.is_empty() { + return (0.0, 0.0); + } + let mut sorted: Vec = data.to_vec(); + sorted.sort_unstable_by(|a, b| a.total_cmp(b)); + let n = sorted.len(); + let med = if n % 2 == 0 { + (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0 + } else { + sorted[n / 2] + }; + // Compute deviations in-place (reuse sorted buffer) + for v in sorted.iter_mut() { + *v = (*v - med).abs(); + } + sorted.sort_unstable_by(|a, b| a.total_cmp(b)); + let mad_val = if n % 2 == 0 { + (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0 + } else { + sorted[n / 2] + }; + (med, mad_val) +} + /// Find seed spike onset locations from a single trace. /// /// 1. Compute baseline = median(trace) @@ -96,8 +123,7 @@ pub fn find_seed_spikes(trace: &[f32], fs: f64, min_peak_distance_s: f64) -> Vec return Vec::new(); } - let baseline = median(trace); - let mad_val = mad(trace, baseline); + let (baseline, mad_val) = median_and_mad(trace); if mad_val < 1e-10 { return Vec::new(); @@ -191,7 +217,7 @@ pub fn seed_kernel_estimate( let mut offset = 0; for &len in trace_lengths { let trace = &traces_flat[offset..offset + len]; - let bl = median(trace); + let (bl, _) = median_and_mad(trace); let onsets = find_seed_spikes(trace, fs, min_peak_distance_s); for &onset in &onsets { diff --git a/crates/solver/src/threshold.rs b/crates/solver/src/threshold.rs index 44595b35..eda60526 100644 --- a/crates/solver/src/threshold.rs +++ b/crates/solver/src/threshold.rs @@ -51,7 +51,7 @@ pub fn threshold_search( // Collect sorted unique non-zero values for threshold candidates let mut vals: Vec = s_relaxed.iter().copied().filter(|&v| v > 1e-10).collect(); - vals.sort_by(|a, b| a.partial_cmp(b).unwrap()); + vals.sort_unstable_by(|a, b| a.total_cmp(b)); vals.dedup_by(|a, b| (*a - *b).abs() < 1e-10); if vals.is_empty() { @@ -71,7 +71,7 @@ pub fn threshold_search( let mut conv_buf = vec![0.0_f32; n]; let mut best = ThresholdResult { - s_binary: vec![0.0; n], + s_binary: Vec::new(), alpha: 0.0, baseline: 0.0, threshold: 0.0, @@ -172,29 +172,28 @@ pub fn threshold_search( let (alpha, baseline) = lstsq_alpha_baseline(&conv_buf, y, pad, max_alpha); best.alpha = alpha; best.baseline = baseline; - best.s_binary = s_bin.clone(); + best.s_binary = s_bin; // Compute PVE (proportion of variance explained) let inner_range = pad..n.saturating_sub(pad); let inner_len = inner_range.len(); if inner_len > 0 { - let y_mean: f64 = inner_range.clone().map(|i| y[i] as f64).sum::() / inner_len as f64; - - let ss_tot: f64 = inner_range - .clone() - .map(|i| { - let d = y[i] as f64 - y_mean; - d * d - }) - .sum(); - - let ss_res: f64 = inner_range - .map(|i| { - let pred = alpha * conv_buf[i] as f64 + baseline; - let d = y[i] as f64 - pred; - d * d - }) - .sum(); + let mut y_sum = 0.0_f64; + for i in inner_range.clone() { + y_sum += y[i] as f64; + } + let y_mean = y_sum / inner_len as f64; + + let mut ss_tot = 0.0_f64; + let mut ss_res = 0.0_f64; + for i in inner_range { + let yi = y[i] as f64; + let d = yi - y_mean; + ss_tot += d * d; + let pred = alpha * conv_buf[i] as f64 + baseline; + let r = yi - pred; + ss_res += r * r; + } best.pve = if ss_tot > 1e-20 { 1.0 - ss_res / ss_tot