From 8ce0f94fa89ff6a20b98cd7b11b4a6e63f6fa6a5 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 28 Aug 2026 18:08:33 +0800 Subject: [PATCH 01/12] Add to_cached_chunk_arrays method to ModelData and FieldSet Provides an opt-in optimization that wraps dask-backed field data in chunk-level LRU caches via the chunk_cached_array package. Repeated vectorized .isel() calls hit an in-memory cache instead of recomputing dask task graphs, giving large speedups for particle simulations. --- src/parcels/_core/fieldset.py | 24 ++++++++++++++++++++++++ src/parcels/_core/model.py | 24 ++++++++++++++++++++++++ 2 files changed, 48 insertions(+) diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 9f8ecb8b2..6b3d53371 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -172,6 +172,30 @@ def to_windowed_arrays(self, *, max_levels: int | None = None): model.to_windowed_arrays(max_levels=max_levels) return self + def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000): + """Wrap dask-backed field data in chunk-level LRU caches. + + Opt-in optimization that replaces each dask-backed data variable's + internal storage with a :class:`~chunk_cached_array.ChunkCachedArray`. + Delegates to each underlying model; repeated vectorized ``.isel()`` + calls then hit an in-memory LRU cache instead of recomputing dask task + graphs. NumPy-backed (eager) fields are left unchanged, and re-invoking + is idempotent. + + Parameters + ---------- + max_cache_bytes : int, optional + Maximum cache size in bytes, per variable. Defaults to 600 MB. + + Returns + ------- + FieldSet + ``self``, to allow chaining. + """ + for model in self.models: + model.to_cached_chunk_arrays(max_cache_bytes=max_cache_bytes) + return self + def add_constant_field(self, name: str, value, mesh: ptyping.TMesh = "spherical"): """Wrapper function to add a Field that is constant in space, useful e.g. when using constant horizontal diffusivity diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index 848bf2fd2..e188b996b 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -9,6 +9,7 @@ import uxarray as ux import xarray as xr import zarr +from chunk_cached_array import wrap_dataset from dask import is_dask_collection import parcels._sgrid as sgrid @@ -112,6 +113,29 @@ def to_windowed_arrays(self, *, max_levels: int | None = None) -> Self: windowed[name] = maybe_windowed(current, max_levels=max_levels) return self + def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000) -> Self: + """Wrap dask-backed field data in chunk-level LRU caches. + + Opt-in optimization that replaces each dask-backed data variable's + internal storage with a :class:`~chunk_cached_array.ChunkCachedArray`. + Repeated vectorized ``.isel()`` calls then hit an in-memory LRU cache + keyed by chunk coordinates instead of recomputing dask task graphs. + + Coordinate variables are loaded eagerly into memory (they are small 1D + arrays) to avoid dask task-graph construction overhead on every + ``.isel()`` call. + + Idempotent: re-invoking is safe — ``wrap_dataset`` skips variables + whose storage is already a ``ChunkCachedArray``. + + Parameters + ---------- + max_cache_bytes : int, optional + Maximum cache size in bytes, per variable. Defaults to 600 MB. + """ + self.data = wrap_dataset(self.data, max_cache_bytes=max_cache_bytes) + return self + @property def time_interval(self) -> TimeInterval | None: try: From b962e8519a88bbb92bc6f701ad38e4ae0dda60af Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 31 Aug 2026 15:35:05 +0800 Subject: [PATCH 02/12] Update pixi.toml --- pixi.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pixi.toml b/pixi.toml index 4b699435b..79a4b96b2 100644 --- a/pixi.toml +++ b/pixi.toml @@ -41,6 +41,7 @@ tabulate = ">=0.10.0" [dependencies] +chunk_cached_array = { path = "../xarray-interpolation/chunk-cached-array", package.build.backend.name = "pixi-build-python" } parcels = { path = "." } [feature.rattler-build.dependencies] From 56ee45a93c54672a1d0adcf9040347052c4145ce Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Fri, 28 Aug 2026 18:38:45 +0800 Subject: [PATCH 03/12] Add benchmark script for dask vs windowed vs cached chunk arrays Results on ds_2d_left_agrid.zarr with 10k particles: - Plain dask: 85.0s (1x) - Windowed arrays: 16.8s (5.1x) - Cached chunk arrays: 4.9s (17.3x) --- benchmark_chunk_cache.py | 102 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 benchmark_chunk_cache.py diff --git a/benchmark_chunk_cache.py b/benchmark_chunk_cache.py new file mode 100644 index 000000000..41cdc6ec8 --- /dev/null +++ b/benchmark_chunk_cache.py @@ -0,0 +1,102 @@ +"""Benchmark: plain dask vs windowed arrays vs cached chunk arrays. + +Uses the ds_2d_left_agrid.zarr dataset with 10,000 particles. +""" + +import time as time_mod + +import numpy as np +import xarray as xr + +import parcels +import parcels._sgrid as sgrid + + +def make_fieldset(ds: xr.Dataset) -> parcels.FieldSet: + """Build a FieldSet from the 2D left A-grid zarr dataset.""" + ds = ds.copy() + ds["lon"].attrs["units"] = "m" + ds["lat"].attrs["units"] = "m" + ds = ds.pipe( + sgrid._attach_sgrid_metadata, + sgrid.SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + sgrid.FaceNodePadding("XC", "XG", sgrid.Padding.LOW), + sgrid.FaceNodePadding("YC", "YG", sgrid.Padding.LOW), + ), + vertical_dimensions=(sgrid.FaceNodePadding("ZC", "ZG", sgrid.Padding.LOW),), + ), + ) + return parcels.FieldSet.from_sgrid_conventions( + ds, + vector_fields={"UV": ("U_A_grid", "V_A_grid")}, + skip_field_data_validation=True, + ) + + +def delete_on_boundary(particles, fieldset): + """Delete particles that hit the boundary instead of erroring.""" + particles.state = np.where( + particles.state == parcels.StatusCode.ErrorOutOfBounds, + parcels.StatusCode.Delete, + particles.state, + ) + + +def run_simulation(fieldset, ds, npart, label): + """Run a simulation and return elapsed time.""" + np.random.seed(42) + pset = parcels.ParticleSet( + fieldset=fieldset, + pclass=parcels.Particle, + t=np.full(npart, ds.time.values[0]), + z=np.full(npart, 1), + y=np.random.uniform(1.0, 5.0, npart), + x=np.random.uniform(1.0, 5.0, npart), + ) + + t0 = time_mod.perf_counter() + pset.execute( + [parcels.kernels.AdvectionRK2, delete_on_boundary], + runtime=np.timedelta64(100, "ms"), + dt=np.timedelta64(10, "ms"), + ) + elapsed = time_mod.perf_counter() - t0 + alive = np.sum(pset.state != parcels.StatusCode.Delete) + print(f" {label}: {elapsed:.3f}s ({alive}/{npart} particles alive)") + return elapsed + + +def main(): + zarr_path = "../xarray-interpolation/datasets/ds_2d_left_agrid.zarr" + npart = 10_000 + + print(f"Loading dataset from {zarr_path}") + ds = xr.open_zarr(zarr_path, consolidated=False) + print(f" shape: {dict(ds.dims)}") + print(f" chunks: U_A_grid {ds['U_A_grid'].encoding.get('chunks', 'N/A')}") + + # --- 1. Plain dask --- + print("\n1. Plain dask") + fieldset_dask = make_fieldset(ds) + run_simulation(fieldset_dask, ds, npart, "plain dask") + + # --- 2. Windowed arrays --- + print("\n2. Windowed arrays") + fieldset_windowed = make_fieldset(ds) + fieldset_windowed.to_windowed_arrays() + run_simulation(fieldset_windowed, ds, npart, "windowed") + + # --- 3. Cached chunk arrays --- + print("\n3. Cached chunk arrays") + fieldset_cached = make_fieldset(ds) + fieldset_cached.to_cached_chunk_arrays() + run_simulation(fieldset_cached, ds, npart, "cached chunks") + + +if __name__ == "__main__": + main() From df22b0838054542b328fba6922098c5ddd341b09 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 31 Aug 2026 13:11:16 +0800 Subject: [PATCH 04/12] Modify benchmark script and add plot --- benchmark_chunk_cache.png | Bin 0 -> 64076 bytes benchmark_chunk_cache.py | 85 +++++++++++++++++++++++++++----------- 2 files changed, 61 insertions(+), 24 deletions(-) create mode 100644 benchmark_chunk_cache.png diff --git a/benchmark_chunk_cache.png b/benchmark_chunk_cache.png new file mode 100644 index 0000000000000000000000000000000000000000..0f9ee9edce245a991984882ee8191a20ba4b6efc GIT binary patch literal 64076 zcmdqJcQoAF-#;o+MT8tl1VJhh(OVELT67^u)Fh&g&R}$skOUFE_uiw+j2b~OIy2gk zL^q5uMjPes;Uv%Zcc163b?;iwUw76zj$>v%ZSVblzuG?UHPjTS&N7@OBO{|yd?c$$ zMs`x3jO_TEGpE2`Y)Y}y;9sK7a!;ML?97~9UpShQslITwx3+V(wtRWb#njQs($1Ej zimSQEqOVzn|c;b2R6696i1TjzVGo=$R848Qnjme~y`_^n-J4*C@(L zKXFT3!kl)yI5g3=YJ8Xd1o;~-RSnV0zw?JrUZ4?kd9w5pG3X^uD!H5dXgjc zO#0JO|I{;Q?xboSS3Hx+MwFO)6M|@G9qFoSpC!VQVJ=RyNaRR4zZl{m@t--zF9zh` z0>I0c^_Ilw$)n%Nd>`H9A>B84rM6%H|M{ob)|agM(w2 z=unHwN=hH#F3b5Qq2fAr4SPaR$~O1qO7xOQLR8eH<#f_Rm3l(@2F+bAE)+0Q$bl7-{KY0XUGCc4_?}0q%6yTOLwRYhBxa<7T14$nz{-*>4{gxIJ7Yh zIz)!h^S9PZj#b$NxjbfITKovNJ9l6G7#UeF%?bHq*WTz=+ZO$74#Y;NOFWL@GPE8A zSC#$Y!`Vgp!*iS_o#Vx)u&}fNbxDlGWSMD@PO*VyiX8?%^l{m3t(^_h=zGlNhX3H; z6T*C-%v5X0fbgWJ7FL3|zF^Y!M!*k6z@o7qIEpnU>%Ep=-DLGz=+81T+)43VQ)F@( zeqn_!PKpsWd4xk#)QWrPv{24TYlRVfCJuk8KTx{c&qhpyUbn#Uuak+X@D3{B7 zXMrurYhz?_q&$b|_0V7l`zCMO)ePYYHyJrOxfhe>hP35RQhxksp)($oKF##%^|+;_ z2?NDXAAA*x{cA!`Vwe}f|%pzdff=QGA(GFOydN!QdkzPo}LN#3HPKo8OaY8L5$xM}U@jAZHt~ij$;n107>|{@IVB~f{PKID5<3=B(8*6; zDDBpq6z|buY=2e|WW2_q5{J*t_~>})@@04c731Yzh(iiTAd9zKQHf=T)ux&Q-7DleAB0%M=aB*@Lbz?DmFstG|`$kA2vTITSB84%PFq?avlTi8tm?{oHLaXkZ;d_v&$YGV zg6`|{qM;C65Lx_X@(ePX?8bR2;c*SZZ3vhKLd#(7l^IV@_#Ai<PDM555*5}#bz9JD@UAOJPfx$w6V0h_fI+xc!{UYQaPZ9uk7erl zU%y_WQ+)OpcBUiLm+1NsEA?ZL$yv&2F2vHyY6vJ&dus@nWj$R*#rK>Du~bNXW_P?` zCuD4t{2xuOcQG;p#sAthkn19ei_I zZtnL==LJplY8%T6!Q?tGZyyv}R1F6U&oq?SOf6ZnfpIlPNC*lJ)t?SK8yvrBHs708 z+{WyIj5Yx^awCxsy4rZ0`pMokk*Q#QD{JCnK?%xT#V6Fx`w`-`_-tFkwu5!|?Q)eQ z(aO6#Cf#U32x{=z3YcB+z0y~Un89V6YK+;z@QVa?)cY}o~LXl>oLh%vsd{DTDpdU>!K72;wSN>dNkOw)V8PBwLKR^ zVX;&(pe7VLCwS$fhiZicHG@0yd5=R*GF@L`gvLpEM)%!O@(eJ)9EPAii+MStn|Ely z%mSroxYWcDKi5-T(jFnkv$!AZ{_8!P0ho>Wo%+RM6IzudUXuXM7!kX&@^V8krm!Jb zZY_mLpMyHIOG(4N`_P=1jh{E3X=J*+qhsxE3Pm`4(*xA{F|>^Y2hyPn0s~9N&;$3y zTnc(Nb(Ms4UYrSzKbT0+A!iM>#8`UwZa44A^?J6KaPnByBfVBz=p?Fz{d6PRPCTgJ z+t6rM1+a7bQ&UU@$dY^)JlSmJl<9*ywE8%@D-j?KPBwL{Nkp7Cz_f@Zx@~Tq3 zdxleMo!liCEY59iFKoxFXjBN{8gpTjUW6;FXhpW9SDPl}_+%BH8ZY98OA>now@kpf zyuRbDkht05{t>9 zpX%xomoWeQwNzc*3|X{@eEJ->nm1ZuX@aeScPELB@m2=kR*wG8zA<9XA|{cTSFkv> zda^UP&$Lz7DP_v+!58oqA6izG&ploDU(Q}7%tY(kuEu$TXEWAj?6W7BN`Yi>Og`bX zCyV#8c(fO{i31p>s;6t%uD+E9<2vaD{dPfN5H=piZ$9b{`F_RU%vyMADpYbqDMie+ zn28qGugt9tYPq1#zUR!$jCPSQ?HJ6K=ebGhBYiUwo%-!vJqr4{D{ zc%x41Q}5e{cG`5{u*@c{K0?U3QV+*wdkh_@<23muLeu!cyg4}m?5)U$4!a7>ux}p; zdU|^KphKxsi?4s$J9kFVXGIr)IHx)yA9~U6H1qI}vcZj{;b01@R%&Rd$0DRw-*u!Q zv~+2|_Dq}WBeD=WzsF=}{I{vG2G`eqY zn0y{-eS3ZZj7p>Hr__$BMKvk_Kyy$pd zzwVJ7VfHml__eVnmmVryHhQ~30PMjhfc%WPXI@##-dbq4o6 zS>oX+No9MKG2_(};rkS{T)3(D_7voopZhuf+xHj*EIuXCH1FmL`NgP9?hSHW8}4ni zkXh`&KPbWCORJ!yK`6|5E`fbrblNAjEPk`;#5*I{%%u|OtS9vJkM8?Z!Jj7&2V#kr z(fy;?cI`Wzq&J#(KKkqt&||5U+x1(M2m88c11D7t)3Y1irmTx+~N<~cC@-YFghfDJQpDsE?J%%znh(<5OkSyi5kW)Qx^Bd8C+fI z!GaH1V9TM~Qg2LB>h+kk%@>;{AN;h2{=$k+I`xP{L9v8vERAfm6>T?7Da)9@Dh`=D zSCd>+osYM^J_kXK!xsnn&oDiz6PiBkOUs*S*9mj<0(>dDm3a`u&!Dwa6K@Vjb~1f) zs&;ffeH;^tLxsKJm!AtPD1Oj0lH@X)hgyapPWjiVJ$$C#`!%hj_`5Hbx{l8jh6r> z>2ay%zIU&1biYe@vhd*bSg7r`pL5CDZj(haf(}13cn?0c=efTV;(iKM8|d`k?Ze(>9O@5Kwfs5q{dxXq$8N!H^( zSGF&IwDtLUU!9yYM{DA=56q}1%v1eAoNo{!C`&AL&byHF$DFqSXSM3)>LbdkVg|>g zH)Y8W3C~!%omw#rgB4`>rMTyMehdJj6!P)$-EcR4*E%PCtE98l-bnW4f!a`-Okh5W z2VIcUt?Pa4h~ngO-Da7W<-F;QSLqwBNI zYZnHxxAN0`(LL_{Jz?&)tjcC4+c}2 z!ZW-;3m0Ace$$5IY!B*0z24w|WP$aa_Z!6CUU!EANbU2U=vO~~(D(8~U%~x7 zPLZq9JsPkYs@2VoDiHgD?p>98T~r{szmaQ#<)h^4BXUg7>e zaiwagkbZJxuxQ)7ac^_71QZa#S0$Hx@fUy~amL}mf`?Q00Rrv8hlywCaAp*k>Y~qX z$!|dVv|+~Ld~NQ>ap#$7(hcOdiJuy*p7drMjQMT9#2-0cl(4w-!h%(eB9h{AtQ6w$ z*qe)HU&#md2I+cOq^VwXGxB#o^_mOuNQ(QBqg!J1A$q4fBEU$psh+y87vX9*c6$V! zDUpnnb^JMfemu@0rpr79#vO2u{v|@zW0rR=mp#YhT-T#h@~Vbkdhl?GlrCzby&+p# z{!uXw45(lboC=OljU@J+Xh&X|2&*A2`~UM9Qia9GTj{IEm%dZ&ojeEoJriSms+NTf zQ!+mJUR8HzlsuN;#U>sf9`h3!JKMvJ0m*i?BRpi*^mTmS>8vo@5GM4wI zukP{GQW9{y@GH5Z(OG{@kTIKHk}H;be>N!9Bf7gif=ySIY^KKr#~p2k8hNh2yVbd8 z!WL=b@zx}ziH|dNJmS$-eBD?uiy0g$xbDg44dycgR#iC|&T3(QS{oU9Z<9JpDzSJfzt0dUpzM(k~wNK5k@KJTH`O;0F7+ zB zKf#nl+%v27GME_ZFv@X%yg>1zpKWI7e22yNt{ditup)X3_YyYGAlou9M@nMibN^Ubz;M|1lFe( ztXl%6^y2Q$L$9p&?fdsSyh0Lu&Lw%SwY$_%;+YdRoq=w$QMt`;8Nc9dr`I3vgskPR zRE1BQ)%zqoN8~Bs#m8#{&AzhIk_pS$)f`-(PadX_qWwuuc4BE1QX)9HX2_{ubQIa@VBN88(Y0m`{NTIsa6H5KZOHkR3ZrcWP25^UzMYXol=uzD(IMRw*hY?6bx z>uXx!bVhDXK+&BmujSOMus3vQxwHWzNj6?x^&R^|y2f8LwX`5bK|nR@8p=?hCXUy3 zFSZ8LPdOp}C^u$ZjwIzq@01f)Pm)-*G@<@^lYjo+^z`&#d`XokiDKDVp#_=VLuK@ zC;ZWCjs|p=q}O0a>oYUILB}21MhLS8#QBU+gZ;C!yX$j;^n1k3>TPE5j70HCu-UlE z=$kS}6gP_;mNDA>FzITNX6aOu(|C^M~;H3`bdezJT|0EZiO zkXNOoRb2&)i^7pU zn66<>?uOkmrV|ZXIchFq10UL*+4(v34#2}0#;NY}95e}IYKs9`T`upDaPi6Z$Xxhp zqjxWsD{Z!ckbNn%;cIEnBNF|!fiP+vWS{(yr#nF2ip;f|aNUg6@#7ndyR6BUXz#zw zWwrRVqKc1`3Hj|9cip;XQu){TQlg4s{AHhHTILg@`DX3mx8VT4*#zBmQcGCn=n^qu zb>T6jxj8}*8OCQ`N2`Bck*SC+ZqFXi<|t&Zh9!wDS44;#t1{KJrShj7`}s@cd^W5E zxviCCj{W}VBllu^Z8b@Duo0Pkrb|-B{ZCRZL20@H^*R*XJ`;X$YjO49yTx%9nJbii zTO>vH!BaH}?sRRopK#aHoxxgmtdx9HFXK^PX6sg6uj+zsgw!}On@E!B^8xVWjm8kL zD|2<*=j_$>%=KQQ1Q)_#A>6gU{k`_~b_nGuwAm6}x74z47hAU2)ia8i&a;q(T`*72 zxFC8N2>c8bcxp!ejVy}rSllET+5HfDy(V!>4A8HJfHrACM0}K4mova;RC}RlLvzH{c-wd4?(3gvonqq}RndR;QL_Rl(}Ke(?Qx?=|2@kdeh-1HCHP!+??h z<0PPa0?zZj@1iNmd`EY;=Tcl$4AcQ*&7AO95x+|%_4<~OH9t77fEiHi8}(Mpb?qr>>Jns&ttDygfZcL><#G9SHQ(O8NP z(b3U62D(N@wlM)`oX6~EsAmPgF|eWKGCMDSoO)}O5>|^N#{^)x_*N7%zkP+$GLmuK!gxRau zn%MWDZz(zDnUEMw=ha^=f+$mS^Rc$V!b0Co#`&tB9V|%5J3+PLphq#BSKBCuByPem zC*f3R!yh%cuf}v~!zB?V2&px(ZO5VtP^)p~mi8^8PQ{f}}@F!^{=a$yNX4$rCqm6q~{Jag5zP@E~O2)<5qX zD(_hYkO{RkQr>K&o}mzm>Qi80Dqm@Pcj4yC9m3&~ z@@>7%J9<)T3g0@+4hp1UiuX<)AYM+Q+f7=~{d@h_Wp&-RW~NLNocpI?60g&fkvm!`dul27#xBwd zj)C?*76rVZjXth_r|!ddjdvC+fA>5v+!%d$DMousTc33@tAum#tE6L|;n6d>|Bv6} zL2NR0^fuE)DF60K?7s4#)h1B;OF6~OyW(e}>TNNi3AWKdW$l2F>Z#36i15VIm1$*? z#-^GgY?oYRJ(5nc8Px8*V%6XqOPJcyeEc{~*9ay@@)zYn?V{wa?ZstC013R5>&x^% zhyz}e)_Guh^n*U~bpwt8K1Mc>6;uJUt&CtJIV~g~`H7Jcb3gj@;dRnj2JI3LPyhWw0G`B7Tm9^<`7&`sC}IuZ5_ty(V-J+&)w zuc;+IAM3EI-;)sU0gSf2 zNv;$0YiOI;J(52qG`z6zmd>d6U~lqZ1RWO_w}8g-Zy%NF#5au=!Ti#FO|IS-E2mO+ z*M+GjHUbj@3RLlCRHGZIug`0in{`>8SWr>!n4W&M3`_>{9oC;$tW z%P(bH7FqhIgQ^hE^TOX|#FG5eNp-<85_XDfetpo>k2|V8qBRDZmeR{u3|6TMoy^9r z`T_EEs{E<^^NYOfoSdq}m4*hA44$PNGx{T~bEY5MLng2?)nA{YDl598N{9c*B$l6< zDGQ9WxTII@g?K=u>(&8_=MWUEyxVov6*iV{b}-3}D2N{5^w! zI_9!#WxFjSK;g^Jno;Ji%LKmJ7b{F&PI7*-xW^*U!b&k)7{5M^T3gyi0qH(aJNJG! zQ!nWn(UQcP7e$eY zYEin)S2Er`!(IM_Lq(XWbgL!^Er6qy=*^us@&+>DXuTKWyx7ub-;Xa^d)Lo>U?2fi zPEKL#VT~ar$Nk;y(JGq+8`uK9fHPu$2MfuJOU)d~c3qh;l71{qTPsk+A9L?ndI+O% zz!EfJH+$**`}aa;tXj_sj-6)|QVC>s&)pLf6B|2aK*w#ggvKHkvJ-ZTwGP*f@1tux zA7Zt|&&-dQjyo#~8+LmYQNEeGlAf$rVex`RSy_2`3bgg*3<}4cm;F=uf|kNfiW#8V zOH}h78(|iAcblPV!3s4WyCNbYQZB2BUcvmd zSK>-Vz0}>G9E|y_oK_Z`36x~?+@(FWl8KcktbQr3cZ@ANZBmVKDL4Gq_3dGHXAK2; zyzL|?9|)@cxZ)YJD%ikUAJesS6)m|9P4$tH7l{_F<25d%jQ`@_c47UEd zX8z3Os-Cs3GCwYN)gExEt9%Znpdbiek?>mKl0GTlGn}IwGDXz?xVFZ3Qp`AXVaJA* z+Z{v~KoY`ufg~uLU5_jNnSCkq*(R_!Mz0wI!j-QF(}_+iBQs;{ao$u}9vUd%$0CU@ z&E5->d^DmrsFD2~lgK#!x+-4em2?~z0fZiMs9fCnOcLU`*0N`E<3c9NuQt`%x>61S zk6|P_xk#_F%qIc;cnngwLQqYSsB6U*sNQF_JUb-8SjOU19E%;cMgXgU0c~#hpaicm zL-i`9w7VU8VWI3T)-7$=j)!RRy~@wFxsvXFaoul}T<*rgG$a(uWpCj6vn1$g!`xdA z072ZvWy?I(>T&)f>aq=x=*wHeA1#>5gzcik&khYo#>Ai*LeHFp_lqGj6(g)SEz0j1 z*VlHV{KI;yO#(KSelbn2p??&|Ee>3l>6^BturzZOF8|_D6?|2qw(RyYXcZR6b^WR? zQ@w#no~CBD)^|$==C~0j$6QxJ4x?>NTC0#*w8wzeKrjZ@N^kxIsw)hFC3FausaTVb zdGS}+FP4lU4eVMkcJ=h2A1g4KzRpQT3YmS)93R1=|1e8Di5_ccJ-TfUE98F z383-QnvI8*>TJdy2CjZn{|~_T#+qUI+s(`#s<5GK)r`a_iR-~da#R57tC*Q(C-t`7 z^E)BZY+0S8x7Ziqs#X5Yt>JZgXku-G3pPB@0^6G=GjuPb5O_sH zfIY4Nl+Wb;fq^o7Xd!S+p4L-9dfmRxcGKrGDdiLB%3~_T z-!n6h0&>7>tfg|vn{Bii(#*O zmc9eEPSe`N<<~Rm4~l(J0}aBA(3qV-F&>vCgj+Q#7y~VPSu+6y!w}tkakjHVO!ke( zpe}tvn8~hx-lLb>&(Zak7&rLBL%~$~X&Xi&!enGLXUH?o#8P-qI3UTa&L94XEadX* z{ePLHh2xpvBM?vaKmV$x{rJN2Ki}?wM9bmJSIZ0UR-sH&9vj7c^cVZoZ_r8QlyPiQ z?1Y-AdDUQT(H_4C-i#OvrrENXk`g6EgpuR#_nW_6n!PjC>=R$V&isD0f7m<=g*i}D zzFX+2q893Tolykd|1YXn`{l@P(k-YoTs0xKmkP;ZJ}qDHgg5YC+CcW>7M@vaEkCqw zKt*c3u^>Y`Jy!Pl#ec8qgJgY#1iVX(v#!W_oL|%PVTIDAh`+b0rtw)W-aY!s|gfVj-%>NW=wg{MIDPZVcS^7QEgHNs8 zIVhKy;NVVH-?aV5j{da*w@eXp$0mbX2H$gD7 z;^M55R;%(QG4~(1Pm;UR+5+YVefBqw_!AkKEaXb_CRdmv+AU|r=dBWg>AS~2X8-jN zwl>%;+XYrCiit{r$iSU78y4f1zvT02hRK0WKZ9&E(<{A5(nu5 z^~~uWhs83$h`n|eay|qF4FGSzw(SL!9|*7`t^iOjg@5`)(|a`AbIV-(1hA|tXBBNY znGOmgrM`3hlRNWN_wn^hwyWY0PfhWGyHOWU`0OqxWc5Fs#4d6KZu7jOE`E&vaEW8R z$7f%Q6eTjqhV7guX=%8csq-N$lJ$3OKu8o$KB5iUcZY3G7bm@L^lTx=RycS;-$ z`qutgQ4(<}8|3>eFJW>QUj7SscDysA{nR9dKjZNj)nwa~OUA2qFO_LV<}$T15iS z4f;AW6SXH6hfx=k1YV?Ac9uAwuwcngLPZ( z*X$M=ozHs446d!7Bsz4e*#(fUQ!LHPF~;7dMOQ;b7w-_OEoZhssBi($yaA9mn7hQ2 z2$&kxhkcMAsp6zlecrzepm-^fr-)CrK*u#^19%SOrzG1#7|r#cWM^fKtdqJbP#Iro z7by9}q4Lw*k^*2%zz2nf9NxB;?(=ofZu>6j8l(YMWChx&QAyW{7F7-BmAcI;=yI_! zo;fx;i#7F-_()0B`{?Yih+0VP+mXzZ zSm!O^Hi+Y{*+I3>kg*A@&#m-6#g_`I1}Qe5uv3M!yDX26T#mq{=JMMAkP_ee@%k7# z%TkpwA8qOXX~@|8GlOFy$dofW3pG6(c({%&X((PD9XE-LkA~;vLM9xWN5QMx0 zQp@6TJ3qaXnykJ%sBjDeacdE<_;Q(JTk%pICzMn<3o1YuGjX<@>EMB=e$bGP@ zsI-r0c_}r+V``Tg9^TOuB(ATFYCFN)6+E5&Ap5OF+eapwAQ}Gmp>ez4E4~X9oqH`3 zpie*r+Ri==@Nm4zuV06Nrqi2Oa2=+vb$jTC#N*cBp;MucYP)x0J)tu;-yM+m^H_O+ z3>VL<=&?Tk(~0F^)91jvjYY7?xnXxrSuE0*-eVBlmFwQ)4Eb#W1L;wdsm69@t*hH= zJ~XOyklU=2UF=%vu(-~{fJHuy8)D;GbJf|&TyA2KtP6=WC?5iMVByPU2!9fC!P5UF z;Uu|g!9&`Mf@5(%Vmv&mAtiNH`33GWp?x&C;&IkHmT&^GxgvbsdSLhMF-E~VM(rfO3Ic?c+HolH?PZ{=9Iq>@y|d1kQVP^Fn@P1C{=dm3Bu*2 z372D8-wA`Nme7L~uw4gj*(E{PC&lg8dj=85OdEPjJwplLI~f_yGr9){)Sn+P91DI8 z!tlC_AhJ^taW2kedS+(ieHaospWRs}cHgSCUx1bxNmS$0;KPq1}L6E z!m~y^8Fml?a_Wx&E2#IXGqZk7drhcFzJZ7X+D2VpoyX?q$V1#pbNoZKB7iQYauW** zl3bOGmy0xCoF)efOm;yz#VI2%!jOv)tcQMviniC8pT9VPEHf_66QtQ)9QY%yFnO)_ z!GS{JFDfsS?M-V-h!D&3<-{zv&aq*5AwTim+`Wz%Z9l zXc6~y3&u==rZEtE7;6y^_O%>>?XI?Q7{~`-dGR`x^*)5WhZOHys!7kTET(rtRn&sO zE9!(5;17;t{f7Y5mjTJ+!e9JJh3tPu0w!MZ&#iTIqTx`i+3^o%H?FI*lI?)1Y_+?!TJ_c(WKz2z z{)NX+o$Ceqgbn40{gA1QaSnxVyn~pk7v05+x#m5|qF{bh0&p?57f0suxh>CzaX@?@ z_?QDA2Wrcv$HyrfNX3{<#VBL}x=7VT7wiRVCEpA$XAmUeGsWSO;&fBgTZik)MOHE9 zW;XB3&D@twk+(tV_Tldc9!525t$aCXnOmK6(^#(9y+H(~^>)8j;3Lzw z+KQLjyP`b$1u8*~`cdcXI-b*@v{WOt1weE1cr*F!8Zz@xY5~c^o0KZNGs0EpkgcmR zIP5a-m2WX3Xw~vi6{&Q@^2pwF1Q%p3btKnmZ-}bjuDi>&^VaFFJ@N2on{MXY61N=a zm0KY-5j)Upe~6LPXCbaAGlPYs=)~(uzx^XG4xjk`?%dyZd6dp(vQI7FLK%H@m$e#* zfJ5*6*LxB*@}m{s#B9Eg^!rZb>Yo0)-%VgJT{!b3w|h<{%N-_fYUrZqg7r0kEP0;{|6FrF=R&+Oga~+K%fYf|+f@=YpKHh;94yCTYlFc8zR(q!tMw9wY zO4yel1Q9Gi=^Oq0;`N@GxVYi!>MFGAIldZdS+jN;KKXQZ8vf3Cz1Rn= zcNXhewY0Qeiwcn*f1^maPg3?^O9_mv3dH40!P)kp%nQaQ~er|w@4{znEVhS(7bsjtc3MhQycS2 ztQc#50=X&sTtI&4Qlvw*0lxaA{(pO31MMrcq-S&udAWhpYctqYZ%m@!qVmsHKv-nV}E!>Bb9%3fM_xaHlNg%{mS_HYu;NY+Z5XVVgkHM&m#8O8z*Y6cCJ}0< zlCLepP7yB3{C&$-+5D05K^sD1bUnQF0kS5YH{j2GJfAKOn~3{iNJr5H?$IV%nyKZ= z-^BvU6ChaUBe30Y1mkw^VgJ_w^Q*DQjtI!#ThcwY*UD@;dpS%S^Sr%dOoSco2Ek zsdj+6L4UDEze4HlnLoQweOcL;69S_uyz5W@RDg?#?z?;GM)058i2u08XSpb%*<`k|e!*2KX){%=26oAp}s+YpsKQW~D>Bj88Jcmqb;PaZv0{(Pecn7~r_($(x2f2j5<>{T5 zK3dT(nXO$rYWjovll`9uY*|$wtomq}{v4TR>~2R^b$7(svbyl<6_1b((Iw`}m)GHc z9>ujPFZYy+MW>YtLxntx&GHQG8}~_B@59%(;>lU1l!=ujR9TQSaw?c(^xNOdvhJAF zFfn-C!lq=B649@xQ2xr^_nQUXH682mn#Um0t!ZVIf9m|zT(GxL7V+rcQHk*|=v!v2 z8wkYQfNm>PD4jf=)?)E#e;2RLd9q7Mb6$h0aDlMUg{eL&-9j(0tC^j$yVuJppD zdfk8v{q zu%u7QG)a#9Wy$tVDUm{&wEc6oJ4yf2MM0qWKDd4BmIX+aH{t$O(ENt06y5O#?b@rZ zdIk??YSI^Ma+RNuuGeM6G={+vi!RZH!M>2Br6t_a)c8^I;q1aPKFbDMH$YwYcu~H@ zNT9TFezXz-!T|RG{812kK2l|)Uk5+{!FCg1y!j|yqdwBkh4ppkvD^EF<=tx#%`9v& zqFe2&0NA=lQ}?f6QU|1?gasuoHCgokqdbr{9H-OgZ;SRPJ4mSAsh#lVp zzI{=Dzh=r_F%s;1E-jXYJXR&0aGiHp;N>48yc%)}bpuhx% zPlqcn;j1Uy3*^Z*Nc&4oV+^wqxO~*B!)8H7_7gvJv;V`$grykEJ*FQ_W5?p31ocUE zrASUn>`$>Ad~WJnJTJbkL98KV^(uf}u?T#rfnsn75K*@Ro3OO8?i)*az^u)fn2NXE zsV61qPq~(tmqQ9WJ3F;xPQq8AwU4;0Jmb*nPgqIA(!ZTpb)ob$y2%c&2 zPV&T@44-VD)6WFCq?)<6+2Dr=HQI0sdjqus``N+RnP@~EcF*{G=!|1tOK?HtXPY`* z$MGg~-V9DbuF%sHeXhj5>T<|7ul(_-$jDbLogE#X#hZyC6G?PO^CT>4XR2Juh*0YK zGW2-X@simh-^nw#sH1PmIW}2AB|Q=yX4859Dusb6}^UX%VBYg8pP9B(H*Z`2m+eAe@VYK{H~@R$NrXi=9}2_cGrVn?jU)akbn=l zWq~lPN0!Y*?FhsRNFIW&I0Hlbrjw-2urPP1&-Ra*DGKl|`?i}Ulr6T%>R_sl#Q@T+q9QLu$JowS(?Fx?n$U}cF>dV&4G zLxs=Ed7?=nMDP-9z2vp>4?7xk82={PC!@>P?ovv>Uuevwu#eA56LrN||2qw2o2u27 zm9t?vP>k+L5B2F~soq%6Bc?>wUB01Ui3*vp%qz;8N{RYq0xx2@@uvLy3r8oX<(HAs z(H%36Bkd{qU?oF-^O9qP47Ojxs7hP^11#PGskds7a=UssbuSX3Al~-Kn=QpI3}EOV z$#n7H$J~SeZUj?_>lirH^G)72q?oxqL1c7kP*P5i+jC`}6ug;HTWbJ7{EIP3kjFCv znd!y8j6g%uKKR7ML=dq&d;r%R4Y&jkrcd3T)@Iao4#hti-J_YHcyo^Q6(II8NP5?L zNq_~;c&VqOBVW+R5W7hGznSBsoPRngezVAJQp_gD=Z+0hTuJ$n)t0x6cKx;LDMa?c z*AczJYdniQpbES2hKJX3yBuD&Z)}xqm;nCNN=2lK<&u+TW$GQ8dw<-fYdatCS?p6; zYcE&UO;RHH5V8eN-o5yXZ49|$1Z>hPOtV~f>SG?&Dae6$9scWG`~k<_gCVgZ{Bz6U zqQ1Kq;>rHhws*+{No5uTY*KW;CIhiHGWIH91r*o*wP1Yb?@Z{36k9wCbKU6+RW5?B z7GEp;3z$TwK6KYS9IYJ}^IDzA%;va-2+ZtN(NV0w@3H9QWz%R9lW3mP>r>Wn=r)OtJTY!TT)j6Z- z;^GqG34_6?iogGFR^pVo(toP%>CUcmT%0SdzmfOFF;KGbL#<7Cm7qmK+Fj!sFzRMM zj#^NveHtkw5HFxeO7(VvFk4ha1nu&TqtbOEr$U3CpIh4i=#CkH)osuSUJfoC)c_Vc8L<*)FuQ zev>(TbR#?Z=bgXPl%d=<5_K>g^<>bsu!6`w$aD3tf$Lwr#k7DhuMO|eTS>P$J)yMx zHSF&kzOU~cF@YEq>d!%Y?M+(_Ib#u+!n^u2HZECaKcD{ByCj~v zH=iUl|H^W6AXoc7!=I-jfaG3LoK)LT>m1b<{Uy!L^-=F_#%*~7+3)Ct1PEP5Ma2`Y z80C?X5y5whzQ>9A5F{2(;=e`-M~l=k``?{JGG>MZ)ef9WP7OZJtr-x?wSJf@qMj6N=12r_ z67r1m5pi*rbT`@A*+*NBR^9Kab*}~Nfm=@H!tD@Ptb<%%ZttJ;z<#^o%{u-jG0wC& z=dYBS<>#y}gq-~|UV|mIPXNk!@}f98NjwXtMmz5Ae8kv8?hcZDmu7eg-pJt)_TI*G3j{By+DrMuGl{#fF^qOF(2 z749Ne)g*Uq^9FE(3=D#QmdEeq2t=l2{e~LWwW0$3W$+4v_@B3OEw>f=>8L$&OkE;A zx<@+5U7-DQk&R=h%KouA-Jlz{f%Z{xyYg+U)Sp%5!@*r}Rkz}Kxoc5yjm%;#o8a5n zze~<1!6UR{mScrsfnR=<_|S@i_wAoGe79)#>CUe7%BaPNw&)HoU{c6G1)PlCr*3@~ z7H+^CmLa?FW-e<`XOQOQ`G4=-H@i6905HwGYhz+wTg&}wGC<1vuZAN4O_|^>+8GzK zK|k)Vr&B19l^XuT!Yw+y6F2~Ychw65)Rx)Y;PJN#K*Gr;{sXVDiy`#*T@()+0Q`Ti z0VK}^jnyj_WhZC2g*wO?I#6_oEnSbg^ki^JkL6qG%Z2~kv^<^MvX`3^lDa1sTq$9`G5SFx~js3 zGM8r|6>Wm+KgO|Bs}h;YqAfdDn9u`HsQ>+MgFZbXj=QY6n2!FOjV#>l^?yr4eYV=3C{gH;iY*NF zcYA3mwu@@#DLoQ2HL z@9Z}2>VF8q#N1PTM!gHchl|SbU z8i0RDOG@tV27k>MXE2kjV!(N!iCjWOd`jU1+;r#Ll2y@Z>Cq@~ZEE z?MdkQ{1rWVKB|DzXm`!@z#JP9;~#7FWTXD?F6XJ zcKti9+{L+aM(!HbNl(wau!pVq`m3h9i|@Qa?#+Zv-29yKG``*g`?F9CHtZyczo$!O zP21*GmnYfKKd!5_6gwWAq{5QWw9^^-q5QQQYU6CIvvMLUmgy>s@cMZD?cJbIF#ASZ2W>_Q^PTII{ z+{}%9;!Z&5j#>S2+UPe&pOHsH5cyfU>|?vQ{;e<`Hpe^YKZ*9})m_AIPKhcE2~s<0 z+Yz?TI;_$0vgRP7y)SDNGfz3srax_CcyI0bG=)embB%m*!^h-lSBYy~IqV4j3s8o`f@hbPxD1kLO$V_Ma zxA>-Jwj-6F+)swJ|H4K=X0pIWET$`KCXc;3$>m?XH_>Y+ZTp$~IBK2{+-IawA}w+A z^4G5l)1uTnXXE}?&?XZ+>6FTtBFS7z&ubzVyGKD>&2hVWAE|@qmnckc49tE%R1@6Wb5~4NqoC8&O(i zOC)2;(Kc!RfYH?6Vr98Ml>M}N+f6WfWvlkIWeL!cGV*kDQfjk@R|hm4+j-1ec^wv6A;eioe>vXX#YU; zAEGnc+QslF+}@p~B_A=zfwbtFSLj$NZ$z(%HG=EV%_GsAY4cpO+uH?`ojsLay~6df z3(3nKoq3)Z7{Hlx`S!IYMVH=T6TJe>Tyxgx=l>4(=Yfl<^9y@*q;i%9V4qZqICu6V zcV5H3-=P-IyMle^%l7V>3p5%eqL*ZZ=t^v)iMqVd2-cK1#Y{E{xk9qlEKyB0=R(I|SMC zq6;WX!Xb4~=4`?LiJ2vxCg>E1xee0n(~ce`?BBI=JA3=Vkr|Iezj3at`B!lA@wo8J zQ@_0#p47Dso2(?}N4>|yx*p7Zfw5f|u##K(w}TuR)F(y{OYenPL!#$`?fI+vhDz~& z>b(8BP3c?6?uJ|RoG3sG+D_+`UKIZKT6{@el+LV}mmA+oK*o-e&}9xKDIDljql4|n zqbtr8##v6$Jn(H{v~?6Jd2c3bvXZc0)9&x=*=zaGPRw(2O3j~OIk!&x+dEO!&Gbs+ znf{zsiMj$z4f$qqjlFTxz0XdCZLe47mBDF=2P0!Fk3$^gn-x~ZD?0G1RSPN|+<9d- zzUgw5ImAVtAD+f9-4)D9Q`A&=^(U7d_e6q!*>+_nK&!$yRj-cv9M`*~s);tKSS4wM zgqxKWX4gBVeU0?xgAnge?6T`#30e91ndXHSbM0#mhQ>QDO@c7ZAP5Ds`CLc-bSKBY zly`f_H*|8;BfhPh>rRDc=x?=bJOQJ6>d?6s1t&W*?9$Tp=(tI}tgbt&`1Ug=Co%<7 zHWR1ShFdya`ta0@m(lz?7fxE9cQtAq?lUU!4_PSlEH?SYXTjRGPehIUPQ|i|kJfog zZ(?w&RN3f&+Kc zF+J~6q;| z_H3R7Izf~8SWQ;03&LsFB7O{?t3H%wtPurF-0ncAFH6s3dZ1|;DZGNcM!V99>twltb0a3r50IT}!<9`kxmD(97>7m3k4L=qV&A!h za_5bN`EIRgyv{4BI~+;nkJN8TMGET_D>R?XX=bkDTjAZ~Y^?bFCCSXt|L_5FsXqm+PME;U_G$D1c{i9{=cb^-i~?B*E8Gb-w6ooPP8vELI{ zXCw}drNP&!gueE$pz5Xlp{*rX*<2t?Y_YziCI)50_)ucEyc$>Hhb*o5?C(tz*l1E` zt(F-_x(IDl%n{rkd$UXrd-I$$P4 zTPhD5hX)>8u=;f4?{VdxLhxsG{BKjlKOSHbg;S^N!N{sec5`tlL&ca`Ku+Xo?zG^A z(U`1-Q-qMvP%#(&9eIRk^TOW-s?gz}!ym4p8E*@HY{Z@^Fli$aR-S%kSdUa1QFEc3 z=pp)_GA^#2Zeby$09in(T^2MguWAT{d>;ig+L*&m_oUJP!5k@;!Lpai_{OlV_nc&= zWSwIDg~hsj#?)J@Hd{8lw*q0W1G2uc<3BbAgqgL`6>ZoC_vf2M9UgB+zap6@#a^?w z3FhFCx3jw?KL~a0$KcOny>Q`z=_)g{G$5TKf*enWWDkbKMD(_R;E{^fL;tP*>CGw{ z+beQN=B##FJ2Uo2;lys_6v+!WYs+9mKB!1GhYoHjq>H~dJQH@I8n;_zQNq%ijZ2Vi;9I94;XZSVnytuSXxFUQ%lM%x~OhqAiva$y%_m%3NTJ zGMcaS-+GHjg3{s4FbI~A+&N=qJx;03$NFGPdHY89X~XAGL)|GDa=W zIhM_RHg?&pR#jO#ddMs^EIDwGz}UB2=kgcp%|b*ussH69U&*%pai!vy|srq9S2K--Gg{Yq@i`?Cg07^<-6euq~$Kr zdXMbvYLCKCHO`jFV0&7@WuePHdF>_|Yup0HIHP)^UsU;%?4gZfDMo!7*c`)R5`MaL zn|x337kJWMAT1^BC8_bOc;|t>A1^UXW}JqWYo~<3Knf-_O@Za=m3BBNjH?oDqw;J# zyrtn{QZ9Wp**qV4y00YPzb`XfF__}zCykxWsRO?tG5YyTX|tjxS58{Iw{#t;6BljD zzZ?)|Dt}ALp{zT?d}iNtm#UuKc=AJ+u5?b3!rU$j?@!Md&eO87vNo?R&yj%Wj)7)3 z7)OJF`OyMW9NWn*6~rNi_)Hit{n6bfuQZB^X>Ro$k3rv&S%%F%5pD9Y_^aFdN5Lr^ zT87n!2^{R1g*idTe%VyAo545oc_m*>Tw1ZU_>*(S>3SLv{QzKr$6@v{BK7G4O?<3_ zCLZ%vN$!Rif|W&D?~*h}x0ehk$eP)Mty~x7=p57W<3F2NLiIvJ)kZQyqj_DL&@LaI zddwcklvNzObK5y)VLVjdRWVFFYeHlwNNXMX;S>GYnuTrO0x$A{ zi{5sm@fIjB6~x3|Ag&z*qYwAm4wN%Pk2Hf*T#r1@KmDDD4wuV^li58WYlp*W^`I{xH8jc2;)TAXIq6)`#z8OaHYN&nF9GeQ1fp zV`;E@GR(yDyh3O&^Us>VsHP1J%i`HuRx=p|x5&y-1O^`tcdQFai-g|o{=aCp_Z*Z5 z&=aw)-;XOAv@>Z~rp$B#3Zx)<4LRaoX19)|obV{g&+Q$MpqfZqU{n00BS$~9I|*#5 zVvVbuyIKJz#(HURyV(2w=6_QNQdK(ma${aCE>_kR-s~0D!baSrU3WUR>d;tnlZQZt zwSjR*9R*L3#_0Bs-!|WH!^WpNrb6^jS&>~scE?SlFrS6C*g)o)ZZ~PV^$@cPWw9DNlBLg<*OQy>zL6PhcS*9tmMlU@)_|RG2Wf%BJf!G zt9owUdbB%*m+et1Uso6^0S>rUM4AS&lP)}vc9lV|*AD*X4Db$T=vKHE$gG3iAqS`) zoJ6thRqyvPfJ4SdXHIVXLkZ8X!(X3^SRa^p`j&~-m`Q2&RV0tmck@uamxo7so|7B8 z2k%a**xV>_nP|;aQ|SWt81SrC)v$iFT(|GvobjX#IrbNR!)vO=K=@F*-S;2HelhSX z_|!H?GqxVzh9rmRcTDX#GXcD0dW@)%yu@T5WAKB^KUyj=mX`X~W)1AK2-g#Kr-qCQ zZa8f9pqroUv%ST~64oi1ymsqqtf8XA$4*WADqt2Vd{f}{z6z?{3NeRJ#K|FD_XWivyM z#TK*|DQ@busGIhwW<~n6cT0toFT3QiXpt^ve^p*UUwUrf*aZ6gpO@mp+&%V6dx|VvT3?YGwy}{d)SGm?cV}7ne#_S)6iWyia8BxLzpRx%y~GblofRprujx| zC`EUrTP22ZWRR z+umiQ2?-hfjXlGw2hm34P zqpA(gJ%O>bq~~1?AYDuS^yzji=+G`~Tc%MMvfOZ)W@g7!Myw3f#2B&SrPsge53a1W zo1Qtej;RqhV`a~@j~?)RdhUZs7<^)>t`o0?I4AX!=kIi0BQ8Gk4?oE#i-zNhM-qK_ z>?YLn0n>pqv$nS8wH;Fyg6$^?th7_GPdqt6vvoh$WsM%f*2CCc+8hvbNq_QK(W4Nx zkvbX(>gwxR;&qkXICJYEqtAcy7-jW)-u3zFrfJPbMt}O<%x(&!y$AVwvK$abnnqd{mfrvB#K z_kUMI_&$!M6oGTnvY2P&%R9Q_EXWu)3%(W{iY?sONo%)aTEw;uLJ zFtJA9WmydB>}guRfPnTJ!!b!nG_VP@AYOT>42y!fZ?GwB()znpR=-iD2kK8DJCD2 ze0iQLElHpt*`BqrgNt(8$@s)ygI{X=BnY5fqwuTgxAxsRiR;ZU@@YfD#})2WIAqD~ z%#HTEl{hq!#-4qr8%-l>vs3%4e}@QKQp=jV76QdZ8h<<($)#LIPONkMHmyUuD+pYPH>JHJgQ6`HjDw0HaTe`rb z2h3-VO!lYm{V1P6v6D6>8rJZPMl-`bPB5ICc{7l-cetzPcG+gLfUG87?%?J(()f{B zbV>(B$i2G8TY~wQXb)8Gs4EgL8s^VWeeE-PeFb(gzE!7~20u+7IQLR)N2oHb*Sc(1 zPY`$<6DVHVwVc;RUlEQopGB4K+-snkDub+$3`^~t9Gz`^%R~3vLYi@p_;w`lSg~t8 zQKUthzLx2bD^8hq=;c&v8Sh7VRQ+M!6DLoy^7A9Jj*Ry+yAmaT88*-*Ex9bWlnh2N4K4qY7fz;+-2ij0f+M_4C85s zfP28I$Ogt5j1_d7Hv{7(gE{W$Z3I!k{^7@S49&LLW@$7$Uz0ReXo7!2Fo0b(Qw+@w zah780RIiGQU_n~ZN8%#ZA;5t#`~F1M=ac>leb3{^g>)d*I@ed>ffO+BI29-ygWK?f z*a^Yd!6BjI>|6@?FDwv}l(P}LJffI)0b{+ubm^Sw_9^t^rpQKMreG>mBTal)HPV(% ztVqAwGkKvv+j}!59(XPqMJ)DyTGRL5412|Y&AAj?w*F|aNYPkGbegf0FCWpD9&f%>GTd6U0 zwKKs-Cl2>3*IQsNGGX*}@URh4sw6Wj>+ZTiTyoskzPy@#1wOfo55?^rd~7#veO@U5 z%Fubdi6-_Eqvx6}3R4VeR4lJa7$n>&U@{&>oRA=fw}XjPXK1xQJb4&=&^M~`y@oG% z0R+E`{uLdzlBXUP{hv3mNQ3LPzwOG(Fi9tH1ht16FRm#DwUlT_`FCfR*%?Krq!s$=Or;Ax(zWZfWG(ppY?+Oo;lty_U({?K|2l@V zZLV8vlNtP+xnUkgi(9v^Z+Tdz*)3tuX``v*UyyG31a!KivpfA}Utchofg|LL+=F?D zDDncpq?yJWV8m-YFraMVO}7d{uz}X+1ppjVV54%cA~-*Ec{*jxwwAe~_U#cHxzRUmNc&x9zYo_N6vSA}I^KO02LEYVI@jkzzH}me!gP?lW#j$}YO<$8 z`P$Jsx{4O@T5dCk6Ubx);0DKCkAh`Gj2hdKwg$f$Z@B{(P|T*`Dy{ zfT2kfgg5q``X^UEn08)^f2PoEy9Trbe8cMEG%kgSo4c7^Tv2hNml}1~{D`w#mDB84 z*F4d=%<^No8|&hSpeWLQ?WbS3%zTbX8aV1I9a^mYhW)9Ry(YlAA3&0K^lo=A>r zkw@!1+!N8ai%=0iY>IQ9lIwo-JVtDg1-IYyJFt!~yS`SSPRoT-n&Y;F3u4TgwGb56*?S? zHfs!L&Glee)W@+mHbUj<2pulj^E%^~>ML9qbY3rT<;-#kyH?LtG|LKeW6BBBMY64S zGdd1Et2ry(S-38@K>^+0@_XbAs~1Unmk+{QZ8=-YBv{UY;t8Ghc!3=Z-6h*Cv_4Hg z95&gvt3g)tpy|Q-?5~5h94x`sMIjlFDy#~W{Lf3T1gN23v+7V4mdY)rUUNR~&UBiV zmmR@6+jSX5-LOtMX8|{LTgYhvpDhjy#xr0cytfp%-%Tar!J{IXZhNw0;|I<=?xSFx zQ^Bfb+h4ITDJL{v<;AgN7fkk6&*hiD`|u@4$@#f`-?@)9xu;j3^tnx}zjyx4S91NC zpeR}>V#WSJLD!t)?mV}BmcLb06-_tmG37=Q4+^gJmf@nRN&BgtCNfuMHv8OgZ*jhy5UB-o zV%O2ifG^+jrFnK0S>fEZ*Rzkwm>j%1CE3thOb}*DHBA@VuW^xcX>p-dlfSLd<)R-W z!@Y?c%0CZ0Ep4g#LZ+~J_`52(-tU3jhG<4c$Gz5>=*}sOi>460CO7B!sW;heW9=88 zXep~dpwle4oOR+0sh6R{Wo<0mwM>(|x|8lXQ}&gI&blbfS?{`=OlLE=$cO!h`Ti{1 zvo@4|B8tyx^WL~nF8(A>#WCnwb5N_J1p2l^Af=t%bW2-~MP$_D>7V6~4GlMiR_$sQ z66)-K>1aRi+lS&=sk|IT&b%SXvpGrPe>!qM^NUq{_fkL0`-FEL!T~ap^M|iiJ;}U_ z3jIF(N<-RN)oZ$?stPZfE)=n9Pv6Go-ORE=%!|B>tg80!{!$0Mph6+(!;Xq4B@W4IOmp4US3^s zQYD?cRX@qt^^1$#a8U0x->{n8b!GXkryHYPr5cH9Z<_HV^M{^x(Hm;KMWLPfQZ`BLe7Wq91%RwcqE3MqpH;m8Z-rc_#n|-!c(baM_1>fR*fJS8$NN%P;w5}0Di(E%{D#X5XsEk}UrP~g?*7nMjq9?{ z&uzJySzT@)LMeOy9p%!=^-#lb(TIi^>IPIsbrNMl8;w?PppK9(S}{>E=y&F9R09uc zzp{c+5v8f4omk;6*#OHMh14M!C7oyIqODC7iV6fAKV7i6A~mD)k$@))N2f4l=#8J? zOY6*Vn!Bsxe?_`N@C%w4Vfztc15Io89E?~&V>H}%K%ZpO@;sdr_s9Glir}*e?Gi)VxX3bW|W*dCc0|A zSt@vVFHbp+!nx;WitZbUmy5~L6vl1tf|LGr>{jjCdvMlsdPNy~v#SRzrhbg35JH0s z^zq-bcxkmxeIXt0a+K=}3)j7|OH0}q?WbBSyyxrrTLYgMxLxyPW(2D&a@c4mdv|{* zy+%;Eu^diG%Wr8wOP1%E#~41!GI6*S>_2k*PjJ%#nxGl2ir-U{n3w4fWkbMe6U zfba*8dGEJ@)Z)2$D5l+M?;#g5b`R4~+Cv+kedUqjKpbshDpM&>Kki2wg}1P4S6o-N z&-7bal-R^-X1Dxw*y_?P0}?DwPF{FWc=+BU(Zu#?!CHSC|e zUs>f^!5uFM;9m6OLCP2m~gW)tC_P7<E=!pE zWbEcdt|p6a?H_Vz6|Ob4S8Eu4Ut7wTH{i*&G{%ZScYo0sau72LG$TZ{W>lJFwU6ml z7HifWbI$GsiH1#69ch!hv4`q@%*}?4xyBym*vhpetktIr#~arK8<*{$TRp!E^d9((ZZo;k~U@e!YPWhxDkWQ&yc5 zcG(B_Z{1N=_6s{#%Sfxj@UF;aAx&#lDl12F;m@aNsd%~X6bn8b?UXGGZKsVDL_17j zej2@s!90Fu?fG$mO_|j_;%I2)2!erxwo0GtPoTt6P5gpEbhc&+x z;oa{I^a4lQ%CbM)>+Yd&_h&u&^Ce|wbcXd5b8-`RM;E%p;e~N& zw>I7RzuT$FE;q#myWPL8yo}i$XR!98s*)@(b07a?)|eoSxx;(rGIVSA-n$-Bbg|>q z7>YN_i&s}hsY;@R8rAr8PZ9!(2GV0{P|z`i>WYj6ZiSX{yu&`q)6}ECE*yen0ELmJ z?dMIS>1 zLMW~eoeb^7req&wt81a#2nvp<#2Kp1XK*1IkS@<|cLs7B(s|Iy!MKH{tuG@aE%MlO#j&@cK^sJcy_? zsZ-9o^a9(Od8)1UyQvOd@A9q&se(4c;f`UcE2F$jM3yXlBHQBJq!^a(x~A4p?#J$` zYwzXFf{y<(A?LUEEgFqoSo{z(B*d-V^Gtj?zyf>4@QTBDJo33;#2?^+%f7(H&(BYM z_43Z2i_kr%U%9t9Zt$AI&gjbY4--4b7yY~1il65(Eg4ru;)L5W)!yafwayd{bk1o1 z4pu6Cd6Dy3$5^nn8SB-8J4YgJq8V0rPB+$QIGDx_5?wuBH&9+!&jimo*{P7;2cubc@#)-Icn;sLCw8PAzeE= z?utF{iDk83w(KA>$x=M_y6UPEOnyIdA2pHinWz>Z-4GBmyOE~K zXvxwsu4pE-B6cNdPwQSN#FnzU05D|2cYb3CQ4+X!_3)&`%SGL6@Yf&1g(N?P5+_~g zfe$716;AK*6YEi=z%a&hoA?9h!-+YlFUb(ZgL1&R1gTewzO$PSN@uQ`OsC)=EV5oV zex4W*Uyup~QmMhD8X2w4B%W2PD;BbZ;*o+18xRTk_!~UwXVO(eJyQBovPBjUy?m}K zwG2Si{f;HA(Ii}a1rDZS8nK5fx;pEy^$T=L$n(ev{4pN&wP*E6L>>9e(D!IdG<}ir>?vKD_T_aun^?mfDwM2* zdI=GMos4xC;#{)gh8lhqZLtPXm>*el0Cfe6MVE|>gmmnVGcfolnu!@bX7z~=v|j>?()oUVhCRaD=T1f!$!WL%7p$BUx>vT;p& zxFN_&tzBKE%hGF>Dn7Zbcc0=UDk{tCez~%$szLadFLz6Lc(yC%Y1r|#!uwKmeV0{2 z+O#S*2W!vRyrpiKHbcp3s#nvJ27#WV474AYLrmL0@TbXk!><`jm6LvoTGos$mQUrh z>zP|HTT>iRDqp^S6)V|h`tTe!IjLMP5TBPZRLs7L80FA*{^B*L%D}UyXaPH47yx$F zotY~2LsnM?Q%_lr@Nbf_*3UU_f*PST&5MM&4mZknPv|-T7 z2!^>iB6DQ>VMd4bC2wS5dA_Q%Wkv6&bFi>i@CE~pSzTf$>f>5W;D}XUQ$mtm)zoXsg(jH z7vnsdTgTq~!SZ(tEM-Z$t!9~fl}Yby1Pm9`Z<56_`R-MaoT~!ZmqHVFiRw|$01w0K z81czcTIeXMEXqV`Yh1Pw~ zrN^>R2_;9qXk9ZS`p6RQ^iaVg3U59&f^<~uSwj}(bEk#Sv=km;+uR6ZhU46f=dB?k3$mp;F=oXvAb>NmFQI*~YWHoffG0!8p2} ziOu~KHNI*2Y?gdeNg5+g*dwa6-lc7moAXu z&IvhDj<(j$J3(v;kXrW-%gsx7KzU&+M-UK^Vz^pM|#z_30e#`W7SW+v;y*>bFlS#IGtaBSkX;b`|&| z8gcb$Q7$&EA!lw_=HgVyu~7jk6xSKr$(HnBz<-}@SwQ|3iYF&$PEG2dZhBGucLBHg zBoNFbDO^K|UfPp{AXT-u{*jl%nQ+y5feIrdg6%nS6}woAX71CBdr;Sn{sQ?PiRED} zq~iNnLQk}46}1-;y-JdkLon2oDES?;1&^}Sv#5q5)hjY7D^mh)Et^t%{8L?bRX+Qb zWZ8(!6jjr<4c#8}icrN%rU14R@b`=S(~zQ1Y<{(J*xQcwvi6@VupvKr2j}kG+#K)2 zPh?Q}RRSW~PzMW~ExF_O;sER^QT%l4(pxNknh=%t$wY^>Xz>r2j>ZJg-ls;LWrq|V z&_9=)(8WPi;6q$ zGKy#SQ93=Fuii8ZO=>0z4XR0T4qIQDSFvmb@~mdUPiNNVuCUAugePa!_UEBc!#+=a znn1a*`TH=PO99(DmD)oOfJ&wC{b@79kl7 zpzA^krF5>UcB(dz?J9McD@FTny<R&9!HJ|T z$~|;gs9H=EXJSw|hUz!<QAyIUSG4JX1)5S)&MzHc4Y!E;lwMP*FWfOk2ec@E9LJ zJ(5)ZduH5KF$$djNg^&bv}%QiF!?#CcTK{655L7Ek5~dP$Spvj=R+WCSWUz0K3#wHzJsm^)~Lv zQ=_8XVfDa8S+CF5aTr7Q1nI#eNj(BRBY-1o1mzc^(qjeZ0BkoI_9eiN7h3!C`Si`F z?+_scwBV$Wk*h-u@pQ@#Tbt0)*}4+qfCntFeX*rU<2Pad!qh_Cj+Y5%NBXFwZx;g) z;Z*zP#-d0!P4A%9C^EDL1`L{&!wAps=I+x-P+P8$-tsC91*S@U5FVioV^S(AT0zu;B`0inU^Xu2sgQ$ep z$f|q_+AXlO85%y}ed7CmqNIXWBL?)!3hC+28DEjc3A+bep??a*i8{XrUxMR z1Plx}Kn6V}bKeRxyt0@Fx1D)tG%ke`Bytw;Q7X;28+sDG74jg!Y-0wq!H z-1yWjC~}|^kI|AQRNa7W{q4$kPkETbs{~-{b4o42r%sZPkr;qdPd&wL@K^GmY7-mp z;sc&q-gdCfCegiFxL%e!e)E|E%3W!9l3X69iD7jqz-{f$%na`^whLm8$=@K-9S8-u zBMzD1SQso+q#0cwb3+M~7wnQfW=SPq*p}!#H%<9&Md!Z^PHR zukyxQXaDjkY*Samk1KS!{yAsG$e%j43)Mv7D-20sF-+q-Q@Re(kTJ86cL?Y%S39ya zxZ}qGwp=#LXQAT|j0|9yB+C1oiJ|Ho9xM?=M|tbxMC zDl029QRBy=q>`r9YU2!9u>KTbABK6e*#BaN7&L@I-tlqb)pXQPn)^ zDb(K=0z-m=9?|*u__*Y7u&~UYy0o=&?uyorK}Nc{0;dsQ84`C192x+>>|GH~0rQ^g zB|Il!C&~)|?E|T31CeG&VTF5+nrxj!udv+7_Dqsu!c}|DJ}I>s*pwn%OS^`0?~ib< z@rb-yvjmbwU};dRWsiS_jBPx0EA>2LkTSj&(+&EqZfJPe#%kcQqvIsI1XRy{}FXI9auO_rF+EF)IL3rXb68=i$5SFn3kU~y~xZ^F2i z_W&-LJTzwj>3F)deQM-a{C?8o=U8Ogi+TC70Ti@9!`KsS=e5Y$7|>4K3f(GA$@}z( zu$0uhkuU_?Y)8rA`mhCoS2=Si_o^GqL>S&Pk%VhYjLQ}mn&WE$<3T1eOSArFiWtTg zw8ty$6&w2D-rfhkR&!IcBG1mVH8&D`_|h+Kz@rls7Z_7Ozut{HZU3M}sCgAs{>Fqn z2PmUqo(FxDymt>T!@1uZIoCv?=G1@Hg(OL^rKka9MDhjvvAw7-qTA)fnM_>yZ0G4_ zSn&J}nS)z6Hu{Ax!?-OxbXxO2di0Wuwb;Sj*#0dX#`*owBl~~PevNVLWz1f=%VwyZ?imy{$U1dHy8CdMM9#pW zn(FDdFMq|gzboGVj|=;6vw9hM<}%n6?RF=D)UtSWs$>#IdP0r>*kB_MEEDQFRE`B_ z0w2ykM#N=caR@R8b#_KYH{R%Y>n!rvMI}SSSr8mHpPQARO3u<_v;Iv-8p@lo6k#aT z!;?@@u!Ohea4nx8`(>3 z|LLG{pTm7#Dr>I))qbUL(DyVf$Rebok0`9_g)Nx$ujzUFpBL(yp*#EY@yjEeK7laz za`M;KePgc1vz6JFE05?g^p#7fK%EM`thpcMF6o5^u}pVCHlhayu?+H$kH~?QpV?ko z_okL!D66nk_hN|-fvW5|^0ou0d22WP@-Ev1vFMr%K~&r+>+nieQDgO%@%fV{{fN6z zQ6H8tCx5_9CZ%{`mx+tWFQpgHT9S0Ma_hg7wt-2{uD%dIY> zC&;!VO8M7-{85uD)v9fIH(j=<6WxV#b97hoX?Q77J3vAAAI0@jv4_Nggd!)swTCOxigNT%j*>Oox=JD+=RUV5 z2}HSx+ZnOnx7q;Kl6s4NWozx}CSucOAaKN!9z#C;7mS#2H2re;5Y#cLmS8(yvoq@V zK@xCG*E8f9l<+{QGHm--9Uc|f94aDaG`K)LZMC#TABW@{D8s$tFd~Sn*L8K0vTvhJ z)1nm+(#`5HLMoASTWss>hDsL?tv)m^?tVd3Mh>|GULj{POsIgS2M| z^3C<74Cd_CeA2Di9se~(Ev2K#9KWPREbID)2EmriPO+zBIHZrX*`!D^&G3>LFK;bS)k5dnzBSHf_3A3z)I!$X$sF^*4hMjV29Vw8NkR*^70c$^a^@7(U)DwH-EPoAB2Ox!Ct9w)d) zJsMNpX?q@Zr0uV`-gsKT!Y+niuz{|0-1%M%k9jW7yucN{Fr7UpPhLChR>9}#;Pkw~ z>zrt%LJ72OP`w*^Ou(~mH|j3g-~TT;J|tzKb762gt@})wTv)xsZRv;B>L}5H2irgx zuZ8rvxeXS=M&3YrnxO4Qa!K+-eoo{nEIYR0iis*4?4#MbGA4$tLg8&qg0|(|jq-SS zSG(UYA48V?-6YZ+mj8Rp@}zA8>9pDnTVX+JwmqoIquVQ$J__kf(FGc`a@%fv5ybvH z9tlgr>*)nbL{x7AwG?Fb+l!(bbOY}`AI{rqL5Loq9(bBsk6Xr172ywuDhG8(px zc`|V~Ql{=GAncVP)C~o=q`u$6biZy*O%3|Plea*{Pe77W7#IB(v!#E(d-(nfG`s8V zm0vj5l=^CEz3DgNJ>Oa!P@*lbZC5B{kF=19L5!OZnE}F$Zu<@6K>bf#RZrMl8x}#v zG($>uITk!Sa{ODN9)rX!e{`AC3&j({TKt!zx^#-?wT)T`WMaEe3N;{|fa*>b4#K-) z0WcH|ufZ9e4r7YPpjXiX)1zX|w=e1?<-z1?2^P!WJ;lR^&#-C>>5d6AAdj0~AeA`X zSJu6PAnM0Ww6tPK`FK1&<2%G(*``FubVeZSSL5M2qY(MIo%wu#1dM66tp&;&@A8e^ zWoYBcE zYXg9pE|?+D(f!X+CA@?YU4PaWO6#$hOiA;k_1~!(__;T4@Q{NRTHfSz90-&^>$YcU z=(!;E6A}sP*jX%kxnh)VN|P6R5kxQ4567MIc7?w7HiT$Lv;3{~t79YicA*ya(?^VY zdQQp}68wh`s>aIsmx$sX_U!eyG<4G?O9GTkiD+%f@-)H7dH#H4*0Da|*;2_Yl4NOS!C&Ph;yQHR(SGW{@^@${BUuj+Ca1U%}m0)dV{HiW$D zPf7p-GADt(KNWa5N`T{QT7uP>bAb#r2T`8V&X*a=5MsLU>I658Je0ok@=z0$b+RE( z)kP+UsW9>JwXrl~Wlp9>Z*3z#lm0K_-aD+Rc1s(N6;V-qt$-8@AfTXBrCUIHFVaPN zSLq$>B8W)uNN)k@y@(1g9YP2lR0IqyAWcHaw|1gu=6rKz=9>9^*Y(RE@6qGV&fd?n zp0)1P?sY$Sd9@cajd$)$D-Q37O;6Su)d9*VmQVhwe z$jQVWKu`{LrNxd@N;NQc0{6A`FM%Cuj!}25VB!D#eukovdq&5UUAAZ5h77(sqgowY zF+Y)=DF_wrArnwlxYtXQQ4qr`br{P7Y&IYqh)nj~BipVmXoR2}+7(eI0(L!^uqxr~ zvwR8jLxLA{DN4wp)-=5%cHFgYBVLp#7{I^Q@eA=Hu6aO^pE;G_*&%$eh;Qps+3+{h zqX<5$H&dVT9S(2`I+n-C#nLGkHG!}hc>BCy8s_$UvriNJ%GS3s* z8p=9Ahc)OvCWqsa?WZ>zd}aCA*AvOcM}2Di265~f3e}L1;;P2yy$v<<#4vG{?j+kN9H&5j()DfGor~ z-yCzFV&tN=c1a&AvIt&rYmK;s{*pocDCRCw9YXyCS;`lNu_xv8mEZgn#EJFUmAm$> zYcbyF=0%~-zMb>$21Ti>p01z@ox8{P30zlkJV?(k1JmwZB*@9-Cm*ouzS>IhJZMj! z@9mhVp(wC9lI2QiJX~DrK#bZ8SY{rpj^GtP&x{TdK^OmgZ8|}|BRfu5f}5;QzYQ>* zs)5%0UZnq8Kx(QpF(2HK&XyAbIw`dWF~%RC>GaGXo5jLun-wRVIAbyG)n|-ZH=p)u zLhPeo=3URbAa5PW$kyDCcWS-ddp%Ni%OPL6HGy5{NeoLL`|3Nr`GO&>z!fpMelmmM zqd!1nIqQ4QmEgEnpyR@z#vhZVo}slzBSi{Ae$J=BZPj9iO~-N|F5360eapI4EkaYS zFZe+tWek}&E$zusVsA1?>bnd>%gT%DcpjtpePzWSKXAFaV`E_JWD>>Ka)GdV@r#f`H3(n^;JB7>Tw}Y8*1s*Hb zb0bwkSVi(V*}W)ItJ#=iL8J+YaJgz&A>rA_fMj_?qOVM=gM>m%vk{iMX}H3%#Fwt%GyzIe^TAkX?78W=GHtfjLzeF)u+cGWighEPUcU;-u#d71h$2eq$<@mq z!w;J9;KLO?dPmWt$WO621+2*u40$g34#78ym?K|{3Lqx}r35uEqDu-@iOB1qZa567 z;RX1deGEn%#UJ;I7#tJw!X@)A=hG@W1*!X3ZmkhP}#Lm4#L z?{QI&Z*MWAIT)db*k#T8)_vt}$6xma z_4%8r(XYVdDM*7_M|R~__u?Hw-cq#K4+uJ){r&LW?F#WmKOZfBU`{PeuAHJWG5}Q0 zff!?yAcaMh#e?y}o6$O1UE;$LM_F`g5mP8J?l_{NQL&UV2-?r;Q|vO=xMmV&F4tGh zcug9=v!nyu0jn2@uwETPc<%qd!^pOIcXMJ@ofRQ`bz*MGM_u>Fnir;p#OXZh9l*hp z<(QW_#>#o>L99l|`E_?5Ci_#at7}rJV@tD~49)+HsptU9!tm}Qw9n2IZax6fiqX7j z(^JrL_@)%Gi1*V9&zebe=a2RbC}tL(J)^nKK9awyM#(`75R9Z7jQZe3*t?;eH9 z%cqopHuJnKJyb_5U7@0K=TJ9qA* za3DC@DniI3zwlz$=LjN&hIO6AUV^PTL^&NS!mO!s7F#!kiNh*okdAxuTcW%10eH2I zl$#^KSwVXifeMH)?aOLos)ql)b>$Bm9xC}aS?<@*>SyeXPE2dTvAx_xpMqRI+&&)S zvc-aq%_t6M<}z+QsU~ZmsRye0k8sW(crvfUg^o4ZbF;*!R^xq-aQLz5QUFRXC=`g2 ze}A4Ix@bp~(G}`rI5<)#Kvwxjeuqe==7y$e{PJRRa5m;#_9Sydxbv6aLnvn>>(q56 z@;Ogsi`ir=eHJKvz0A!yfB%u32#u;P;R;5AYM@saUTqIW2kfhD5W_s{%+6>e5XBdS z)rMrnq-95L%ZW+Lp*=77uG!CM=;lOxaNvW5S_;|?5l_`#Q$O+8 znRD)%{Iqy7?-*?l_Rd&Q2JR- z$LxMV!MDFO2j}lF{t(sR@Kz2HX5{(I-VN1@_2HDY?9NQqAKW&W!2YTjo$cDSfmVMR zggz%Y%CNbt7)0Qead+Q?%L;?JtAha5X?K%X)6eqDqKrsgegb_R$Uz*#9mXpKYR`A{ zs*ktxucu-k?{_Ch^Xw^H(l4-*KjxX!oE@uUS)(N_jVIV{`&;&P#?n|ndS^a&Gqtv*im(!CRj*rP8c)vEBBNwAVEP178oKlNku;oAGscd=$bmyG_1%(!V9 zSAFdeltKUV=g*kU)Bc7w+Ut{>he8_Ly}vpQ&AQdz?ZtmDKd?;SXz+!sGA>U$FO?ri zt&Sy)#!XJv?)xmplU^)I&0x-DymZ5+C-s*>>+$ideVGXXsZQZ8^tmvRp_n}-QLg8( zx`(v1v;Rm{EVeZLTf&rhpvcC?*?$$H!Xy3TpFe*d(2@ZC9xd#govskAw&Dl_3}lAD zZ1}Noa%q7u$@N!DW-&;pt99)E?AEd~I~nw-S^QqSWX2iJIni4a1GigW&U&;8InAs5kE0Gpu%+SQGO;H$m4n7M!!>+Z-V6Wr zHZAWamykStTK}P@3s<-eU()C|i1GaLK0*y{e8_?!AN~A)7LWro5gHO=%qV^1#*It_ zX8B2mMat<#YS>5cYY#4`X=ok0!`#R!Yp_>VtWEgcxwO|!vWqlNuY=Te9iyCp^y~|f zrI%LBAIBJs3vs_W^qorUw(TRV``w8=IPV?oGSrvc6hW!lQn&c}WT)cX`l)(~xB*tgijpNGGKDfM|GBV1K=q3s_9=Uu z>nreu23qn8KhbiJF7~!s!Tt(5kcN3wKm9`SP31vYk+F54B5%$+0h__Hll-cFMWV+? zY2JFpvfBO9=Nvc*I`n~B#(2$jm7@m*2Hw4Fb}2!o_#UCG-|WPqOqI|3!s|U+9U2W0 zJZW!-QCu$<^A1?!_6gkXJENu+qe!zC){&pANBzD~fjlhc9kxcp#&GtV#qYO#=vkG* zmC+ddwvr}P*L(qSM_DY%HXvzz)jkt;?a{l9vb~14S3^K6Rt3LND|9!5C%DT%Mi2Xj zYtf^q%;N{+xhCunR=>uHQRG3_6M9==pDEz@&?R3BXTePBXY)Eo-SX31RUbTu33u+9 z7-Bqknrk{lRZz4PKR`wME0fjM=Jeekd4yerH0&y(inT0f=3B;UhJ?c-I`{z=qn=Rg zqReM^%(FdkMHvrY)rCDczGj_pWT6;Y&6v4Af9gW2@$XZnC3OFg^e-;prM2e00S({o zoI7(;q&z01Wy3#o3Eq}1D4;CsK2GvWbte_e3@)1?#t#$mSxrjkuJieHxshj@pa|5g z$CUrbgCCl_P`8o`cVQx=Ms1Njy59ZsEo`s@hcn~1_Pujgbf+&%22j#UXLk{7;Z4t8 zcaL@z-mJVfdiO^N52!QbVFO7Oq6yC)zo~UJH|QjuKTlPv zsw!1cyj2A{z5=hyV>_bp8q{D;}-WGRJ zrit&^5z8|?JoN^au@BG&6R}Gox9Ek|p7z)kjjun_%m^Ll?y@=M>TfbUFn01oW*5ouKC@LJURV=?Q(44oTwGkKb07n~_X1Wqeb(eY&uxDo z##EX-{mao;!`A&lQ5A>pFFVdTVY!Lu=h{>^;8IloGvpR9QurvF!uwig7=w8HcxOjP zUYXa;sY9$?HglX~abssN)@(DEwv6# ztEPumO}ti9`w43sWzF1KDfyv)sUN#^%5!xU0zeDD^c-ViH_rYwYO^MG4ApzTv_VGn~kK&KJ#UnYK$NRL{ zG9UKI9XEStkL)grfvlOk-}oNYdYD!2p_qQ+j2#Pd0XzAEkTk%ZMOrT< zxHUHCnNz)Is&{nN%sf{2)qjhfl*tFW5w^ob}LcDZH3w!iZkhbDFnTdNgE{bDjNRCJ;}UxH1kMSOVwBab4<;8pWS>8>I7?rEq*?}(xI7x6$M zC7EOv@qIK!_>dAhwVbr-9&yhC*Oilp1OMNy`Phn~1(gT_vlOefm}*u4`+Lg>G`@9+DQcX{a1n&K&9Hw04ea zDSWM=+j&N9XjDy&J|6)Ceo17rGp((9X2O0UZ_8?r1XtfK0rhdQ<3Nhi@5ZsQ&tWnn ze2;F$3CkR%62Y5tCVc0S-?~nQM{dXGR`!aaS$-5w6Api;%6N_u)=>6cLsNkdwD_`` z9CbR4I1_%r9l$w}o44Pubg=l=1M^E=n|2?v#%#}GS1%mTSSwEdkijI_&#+su;~-@N zDGtjE2v-#QR=Muq?4Y0Vo7d&g4AVMp<56oK7&egNvN!XZ?%JDsCgDccXvD*<7Le!C zC!5`nFJMc0Y@Ze9HhLtrrX`Dg$f`4&X<3o5$=zcuc=KN5zf-|qOX)D9&fz*qwIz~~ z@}9WpW8IT|%&{yb*9zx#ldR^@?Qf79v#22|w;sKJO~b;6nBG9xG~~Pl4~rJo6RAI$ z&X|sfJkI7tn0NWHNZL3eq8I8vBMOt&8o)@LX883kGwYId278>L!%By(^(!E6N-PC7h0c9hF(b$PE? z&dwz7v_N~DxeW4je%1)?d4Ch8!3~M|c;u?2H;i5uMt@`_qU) z@Y_e;kjcyok*m4rf)kczC?^8(%d(Arvd|YDkasgdOKUyBo5<$*nIqyMX>$BK2)j4K z*}2W@oLu{rnA2TL9wk1t8qtgxf%zzhLJB=(9}MRquguc_D-i^X3b5KkEjba+A_?#4 zr1cMj2RYEZyNVyusMv^grKFy~gDCwKdj|Q2pB?h5$)0?v!$DML{&{a6YRM=3oa}QD zyxAhWnS`=r=1Ae{`9i+=C-S8?7i*^tzPMi0CPp|MEAvcqVv)GIl2ME{q{6zIGIOvN z@1w|&UlJ00d2ta6XA%U$)xYHIde<^^C+E%$^_gG?)e1UHcW-o(-G*G2-J@(&PsD^nJmSbHn@S5FE0`b?L3kmWl|3S$S#DH061nm7)X;FtsO}p+*nHaDx;i6yeWj^U+x{6=drua5 z74gpxhwD1uqa3am2-0G~xaja)$|Y43hW#HAThE^tIePwL3%Y<#ZNejIa_X+#v#`O@ zPIH;tGK;dK8J$^ef{LAoI+l=iLNy?BnO>Y{TUeiZy^>(Bv2S4&t#(gtVNk=6b8fr9 z1%>shb8x&r9qa)8u!ejnI1_Ts`+u+63uKBAZMALh)=*l!K*opsHxJv)WvCu4?7pIA zemKi+f*v_4daf3fhn1{H7Hy|R`ZvPn*Cdu753?+l?|T$UNt-PcNuRKc1zp2Q)-@=W z3&@r9qX(@Bm6bOB^Ex(aw_?W}3#dn(B@P=N;79MGNCe6j0ahwQ+RE}(+U~pFv2hg_ z`v=ABG2Q?Sg5{A+zL@EsPn0!821BVGrUyi$Z*$uc^uUy2Qe8*6gjJ8zT_Hf=|PVSe_k{<#$B%JRTM=*Nv^y z1bBdpt?Ii^mYn8A?GrK$-03$or_0e^F0U&#LRycybDY_^cK=?Ip^icDTj|Z((d>sV zgwk4b3$tF57h64n zd?u??lNw4;HN>%Sus!QTECxVB#kA>#bgW0>k(ercsNlK$?(*3@@9?~(`>R~T@LQe2)O8GpZ$xFBB#Hd&v^$Yi zxYHt^_PJB8h!DBob%Q6gTT$3O06S$Xs71yeEKiViaC8^Nv>Z;04ow}&**cLWW1W`l zU#H&EKPJB@yL^rOduI?h6Js(Xgb5T$TUpvn3=WD^Dp&Iiu2oBjK!@LeUL543jhI8l zkTNOavY{28i_Hw{8$Bu!-l!N-=-C_^dQ*b@WM9rAPlmce^)pSGIa1sg>5(+b{D>kn z%8fklDI?#h9GS0e85r=8^^&g?htZ}YXD<6jk*0;!uqpn>Qd<8VB1OBMM8coXe~vd{ z0;|GG9iNhxURq)P9+SW6!^K*hDnVW?a?m45TqHtohbnsg6LP@uzo;RFo{BA9K1cpJ z<0qh6PTP_WS>e1fZ%>|AeoETD=2h+Vvj)K6GJ?=S{UFfk*XavkCq&uW6QM=$*IZ;_ z?r4C-7LQtI!7>BoA|zJ_rP-Q1VggP@&yO#lE zNIK~WX~+0HFW%B4O&xuk2C2-CdFsuzXEXFDVr6>+xoEINF6Nr7?uxP4o_nb#zGs4J z)UA8p{rr>rxV0Y?X0rVC*ZdninB)vrL|loJ^}Qw;ZuLuRa_is6%D zXTxWgjMz4093+<)Lerv7Wk5|W@;MDEUD9;PEv$UnAnl@?%wXMT-{|kQR(z5hmY@_6F?ntpqU}5UVu8 zSXP}v-GQ`FHd#7@I9tRUhOh=TJonk#-G){jH!+dQa>vuQTb&r$;vXhr@os-sKp*zS zZ$aBdl4u>nKeZ$4YM(>M&t>}nk(Ew0LgM>^zsR=&sbdKQAzmZMF)ObsC9Y+d*-RX+ z*7|&#q(3kxEIe|9eAyQ@(DeJ1&2OgRUF8@g)_lbUuhV;l^G}-}KN%+S`@z3Fga>~T zIUIyxUy~WDG4zcVw~!y7kKpxm^X2~gL9RRH7e(R3(=r7trLOyip4gP<MaQ^K`}3$J}Yf{KcwRWKcNyK9VVLpD&|9;aZB3=A*0= zxj&8AI1K2WL-~Jja=d%mbBSzIQO1Z=^}QS_x#2-NV>O~r+H-F_%|B8;#UWg+Fpg44 zhkAj6CogZRMaymdQexd&h>)dQF`E3lCoPl^Qh3<4l8ot5ERT>U_h4NIBPy3jx_L?N zz*Z#wexR#K8q>LWil>Wufs7P?q9j$Z-2Muau_n7NvKeXBD=TdE7E8FnA&A4DAgd3* zeE>uxEpyj3Vje^mn?%Fm%fh^PiGkGe%VwoAatRdb)>YW&xRNDLtHyXtF!P4k+{>4m z4-yab{=|Bw7Q0ka9EV=sOR?H}$two2 zgbOIo@Ds8JfEq+O{LuFeo}SWt{%ymoK3UGr=C8*o3bGXRLOs**;~3X`*~H;pY ztwr!qt+}J_&Ao9qR`SjKOkjuh)k3u=-c@F(Hfzc#Ap)~7lV;o}FPlYi&rtOb;K_C? zcLZbFWHBvkGaR(%Xns`)$ro$h_8?aXRJrWw{(%3)MEr$t`oxhcUyHiq^b(ZH%Bi+ zZUsL3^xy|&Dq>Ti2(Z=@y*+w$TES{~L`_TCgRD1EPau^g$q&Kwc2aEdf{~g|!Pwcf zx!0;67KXyRNg7R*;y={glgMS9JfXEdN-^Gb8t+?uqnGnA<;^WN5^@2HRhSFppOfqd z59roVqg;-K93KH(l`nM)YFt*|gD_ugQlFeqQ+nNlro7{OMt@ZnU0?^TD!>d50P6T8 zyNf3j41AjDpW}|Y4NH19ZPBrYEuE-*^8xGAyYemLch6FuRQR#%$Ydx28IWDOJi(llV9y{^*Wnz`VEK z3KXT%xyp*vF8PX|k)pq%O7DRVi<7 z_vf_x&A`Sq=xlRattJ9kkWF|~kN;q)PRB6%sJPr$+{PWm|NY&4|5C5quby-QNlp+T# zwy@a4+4qCND&TXh)TE+Gyd$?g`i{(9t^f1LPhc99-%uEUH*t=3ZyC-GuzER>VAC-n z8M0TdD3#HYABbfZxmCk>NJQNwZ?Cdq*BTP`{I^V1*ZCTDi2d+%i&)F0 z4{rN!H8jd`NrVA{cFAAWOFf$ve(yF(m44X+0~#rmMCSO#pylqJJ=5@Az>JO-GWq(m5*O$h@{>*2ED1ak;m-`@e`Ch_9j1YT zJ>NdRLgzOzy6U5w+BN5im+xJYKHMLNLP96mPmaX!~;evI$~LMjZ~Hw)KhK5 zBYlY4kCG|(=e{XyLNn&I6tvQtd&W8)!)yz@C%8T!i@m-PxgH~;^E){RlXXHF6%{;Zjb|Tp z7M1ya3`Y_ZjRyR~|4pHn0PMnNB2p+>7jO0`tP$8%auvhxaS3(O!5AFqOn;|s=>B>C zJ?M3f*Ls7^y%MnRYnBS;umgwtSEJ+F?>_=U{$9#?`eCS?N8-89gF=>+bKm8)TiWKr z#bt98aJ_8IT^t){3U1n&+&@k>w}>?4vhISvKR zrL~vu`y1q>M!U-9h-YU=27(<+EpU}X3tt_kOg=>MQ~7t(!#qTEzKF}OM9vxz&uqo` zC^qH~`yZwbtjWJh2oaswZapn1P4AnCtxfDf$_L>iNlSYQ{Z2Vi9M=bR*3~jf+>aB|QS7~&vG7KkwNiauVpdki4!Ihb8K>h1w;pq1s^s=es9hfscX?~z zedP7OPvqVZn{NoAsD}6dTom=n8{wMWJ7bns%-z?JL2tBQWHUGgDRHaAl-Ki?A5ICGY;Cw5 z--Z>n_(rEX!QMU2fBP-;ALmdXQhNRD5okCfho|-$o(bAc32S&r0r2FcWkYK% z8vJ@Wc_TrN6i@`~N8w1P&Kmv!-YDZo^%rL2rww)Qkle>AcI{b-l$D8Si%9}^8Nj+Y z>A8m1_d46>CaH_QZ#@ES%|B9fD=mIsV1xLdbON@pvCRH0pviV(Xjt7x-)OCP#PDu$ zhut%Y9oDCu`CEr~rra;fXq2z2UGwF>qKP2)pCTMr2Jn;+C5k2R$f>7Q>jwfFKgl(1 z$x8i~*DF2s}YPR~BdKA^?w+~8O=UrAT)aYkkQ0~oAe($k9Gp)WX$1=rBLOpsN z_+L~=mo`tx9|%WJHV@S3$z4!Owdt6lv-S#c< z5d$Y>FH!w|S<93g{G%Qk{;K+C7CJQW{I|dR?dAvy8xb38!-zCzLuVIo_)j~jLvvJU z6_lO)_mvv|w7J7ek9tlHFAYEc&s%Sh7Hl9B8=7W^y+inaK>goon(oa3t2C-2M6Mlx zBAKhE|1Zlf|1^7g|Bvg}VLe~SC{Q|;Ar0*ejyh$rDQ~CfRRs^`YxxdUG_+J$LBb`;EYsN?->R;w( zEMVfTG)#OLtn;nAv=Vw{`X01-;$9Jyp^s{Hb|3}36TcZ6PCgTs{8mbga1WO zSH+Y7!TG^{N=w(j{qTX*i)+>Cb!;Jn|HuqkM3qOq3?66%$7(ZIudKtAsB1J=3vIth z}!b)KAf7SOp@80fd5C@gDax-alq z7IeQWCWsCYt=CFiBA-`nwEa`7vHdo~LojtZe{@OT@TwH$?|_2uPdi+?bJloR>Z4d^ zyYEf=iB)mOf9u8y@i)L9;RO&>Qp;^r8>hYF?(mKyL7q(Y<-3w$L;N0`RQmhhw8g*v z#cCKm_gkQhCzmwz+CLi#C-X#A3)gaWxEb7EclA?N`Li%>T3@m=5ofW-nRPm;s`yE z{A9-Nz9o7z%6;+kNr8voWdPrEoI7_jBf%rP7I>#rglJw^0{LUCz-U`YTwE?RDoz7s zM2YyZ-FNGtLpR>*XEVL9Qzl?p>BwDzF~B2=#2OR_kCa{F*4-TfEqT_dySj0=m~`%UcZ>J53S7aoBEvARTO;&#U*tp36yzTPeP z*-~-$3!I#VI5T8!=plE6R@FSsHrDDV%wz3R_9d#@yt6|Eg~rC_0AIEh^ugj{n~g}s zzUS4Arg6~)7U>_$YzNdHQyng~4-7I+*NwVj&)e^3Y1Y^Ul!RvJG8^&pznYbOFqtd? z19!hc-}w8+H#(ioxzTpwwwb}vat`T&m+ri;Cu)#@NWk-1LE_EFBxW*vE%=9chOqa6KEH zm%FGA-67C19AzH>mW0XnSc?&M#l1hGu9$3cA1HFdtF6qB3{?`M>HLA{aeAz-4rYG) z|Jd?Ppt4K`)+5iohfe=^85FbdA@ zQe`!@zw+sl2SK_f>EI}^b93ie&!%Q%F#oKCA>2_B_cG-A7*|oEs32pa8iwb;M#Be+#UEHR0z(KhEojV=EKn* zP?)BxJ548G9mam<%sbq#b=aDRf7>`z;@LpME5sMLg3s#-O=>pbTQF`w6n? z(wO8vMtv)VirK3Rs1Ra>$1Olh&={`4*kfW%9=M7)@8!EZHOp@;fXGf$qHgnND~@!D zP%*`%=fhLqj>|Kcsrdd@UKHUw5b$C0`0;LUKGdErZEmmC5d*C&`)wWamRW5eDQMMOMu%_uTylk!Q3PQV+UpC*1hCX@$9|f0E z_qu^oz)zZ4GGq_tRA92rQ?bZDLq}o*=;QPA>b8lx6YbZE29ly*ejVEUY}#rET&%=| z74IWxHQh23( zx0e$2CK_Ke&Tw43nCs!S`0X=P=w4d`E>Bzq!q$^uKT@eP@D*xuoFtMXn&-WMPXI>8 z%@+6ord&S5CvDLdWgLFRq%BZ@nQt$Q2){KG*UTms24od_W>am^M^nf4?mIjwZn3?6 z`1IAgTbr9poy@!@#hICZ_f$IDEufLT2f_1zTq{69YQ2pKxks ztj;%OG*bF3uXML5n>yCCeRtt@D`MKbr|rd?FWVSmmvfA^S_a)mO~<`{eSKcdD=7wZ zh?AhtU!YvvSaHvuJ-B{81zcb=F0Z2bt%Fzz3FIorc?38z(JTF!o*n^jqSFi)E&AC( zx_jQ!qi_JEjm8<&%>#RfyV?S;^PQjxtudeXd`|yyFM}G(VAMvO4&frGiSE9rmHq9x z=^qbE=tI2P9qoujs|#^$#Cmg@TZuqiQOBu{bv@6q(-AF&)Pkb7 z2aAfTU*<>Y$U>;xymy>U+61>-P$!9PmlqIRXil{SvKA%e>x5!6GBafu1Z^#yBCJ!v z?Yix~J`suQgP6yLn8VKWXVJ+

%DNL0QqyjP0vyh0x~Nd;Ok8FPfxJ#4yB1YdyH^ zcuoKT0@+P%?U=yZONmVoqVBBBiV<%6;AEKepiSCBx?Lz+eDbm6(%ycBNFft!lfq5@ z4{DBed}~*iuD^>W&t%#bUt~(XsfEpo4d!9`UIj$ z8_|oF>H9at5ArkLjiNDASZMT0GTxapxobU$L1_8|74!EARLr_z&-udjx82s;5}g-~ zA3#D!uUXLCqy(Wq%+ARx>@aX_7f>A_?3VtT33(JI*D!DK>$ zYi2>2$I5he)An+Pgif@rITZ%Rv+GuM#|r$VxL47n=##%C4!agZ#m9WhX?}R%*@>MY zBJ+ntW-r$jQTJdl*@a)c%wVLWjrT8$LYxB2^5>IA#k`-|H?;WVcS$a6Cu%R%q(?u! z;y9@QNiro!kd@85vT-!-R)%ACEPTW;e*DnY`ZelzI$!#j$Hu9Yxi*1s;It(pZFE~u zN2bI6_t9q^MnA>Ob^(ay5E0R{h$G1Ak~jvk4Ru=yXGYM z8V%J)ZJ2C<@U(ibq+0fm|>IY=u+;l z9?tuD^wv?o#ip_ucAYDd3%?%m(V(ejW9{#Z_3!Q~eqhlUQHQSC^2FW%u(u!91?H(4IKe(Mpb(KUC9kmgJ3fikJw*1D};p6CZz%6~cBx=be zh>@gp4Dt#9T~4;OoK!pG%nttef$mC#?VjG#Tu~79;e)3li_^G0aDe(m-Wx11dG10O zwQJwMT*I^H%P?sP8YnJ9lqzgMq2aQ~JX4!3@h5P8_Hj}rt-qGocAEhS`SM}b$@_BA zLd_MJu29La@9Y@&&z=nBe8-N!&wwpEV!k;6BSSYq^3AW&n=vbRm_{@m=esdiW)r@v z%WLxEN`jcDE-7)TCO;PsjH`vU3HKYQ<>cg&TX#^6+Qb~@fU&W5-kqc+dZ<8qjqrG) zpO%5)v)O{2QNHs$%$Y7a@m=uM{a8lo9lLJ>kB*Y~x~YvE0dA&Saw4d}c}7byUbvUI z4Q!yvt^O=OgE$=PhG0xV@3{HF7|(=t4~#5(7Mz!55!&-$lZtloi@45X!PKeEB&_B8 zZpj;YEk0FKQ*#2~)IB9hyEUy{S6CbKpbW__G%23IhjFM|hwq?0ebUhF-O36eaSRdw zpBV-^V24W%#)VQq93m4b#L}p~G_sSH<)|s3v4UEm4;_$e$OQ00&ujW(?P@(mkg`#06VoNxc<_wPiG~<(F7*4UIfV0xQhX`Eo3W=3Kr>94i~+eckC3?UMYgo zfM4Um#Fczs^!C<s$NO|v#CuUcnZAd;p;PMyHs_6;@1V7@+X!0o~$ zge1C=!6CKP-{Ycwxe6HW9O62_!cuY5UlTmOI^l~}J9!0~jD37+1#oRO-MQDSJ0Q(q z1x^&sW7XW{NtV!=GyXnQ)5hZUOLq6nAH`4Bn%z4$`Pk_k7neN1j0f?Dzad<)Htj^< z$r`S_-H|9!XC0Au&Go~clU-9ZSNv6zvcQ;`@0?j-0#Z(&d{5xl5(Q0|t0@R{1xsM} zm+qmSUv~fbeQ#6zwNc24ZPzQ+BQ=2 zA%S^7xM%9tTEIB>a(K3IWVn?&j49RadQN-tWZ#$a8SuK}Z4jsquCjG}Ap~3!5S%LR zup7+(MfkDoa4B`>rI+7WBlW+n#SNsnj#j(lZVa=9nnNzsx;uGK{=WPB*p4=yrde-y z@KF!{tQ8uqYl+k5OPq{dG%p;0xw9Hvjx>^63NEXuvZ3^|$tJI_HMX#b-GTT~$>_W0 zuq!iS(N?un37&bt1t%Z5>-k+AS8m0zrxKGj>2=Nr+THjF0_|=4igE*Ww6SHLf_}v4 zgEQa*3iL}&ikSkx(`Ld3%%(t>%(V*-0;P0U$ZaNB*IhyP!Q;6@V3;uxT&)I9t*wMu zpraDVY51h0fXi&=pBt?a*&2l=xqn30D8y6{Qv3fc{HElUi4<)4{k;Iu2|N%;AtUp! zKvO4`Ov50E@ZMkmL}?ltGO5<(tk<02*3V4bnzfi7sdD`rVOV<40CiqK4TGE~@vUX` znT;ML5k>HcSM}xr&T4_^Bdoq@>FF(Y{@ozbswRzrFQe`nD3nIMVN1u{ef zMad-?I$H&^N4jR~EiVI&^`E;(?v~C_h`#4dg2b9NP_74nhskhi1hjl!Ez!DfkMud9 zpo`9{niuA};4v9(kqw!CaO4i;B>nIt)MS9igrAZiRlL6X?S=4^GKhznsVQw6WK$p! z28jW`){AGneDe^Th~=;*)!$6G!W zs%VrAKa@|&^90C^P%T5yr#t45L1^?|A7GW@;Nw&8GxDt&YkynmGG~o`V4c3e4kZ(* zp-O5%=bGIDzJVr;7VvkO8%zV6Go4%10ZLu1X6_2Z(P*N=7F(ByZU1|ZK1G2;r>`bm zbDXqX|2e_Q$*DHg5ibB#+0IZsqBIaQHxCnq;Ou9B@bJg%lw}h|**eXo@c0gbJ0z6;^L`r0;HTt`q002)TWk z5FMfy^C*U!XGre;smoxlI0C$gmY`E{W8vN((%&gDJ?SeCQ$YvJohbWr_VV&VFF^!$adxnz)xjIlOjA#eDrIVp zo(6QD&JL+!r43~Mg)%JcOc(?amWQ4D44hjw5~jfFEgfVOF2)If3AWYnC4nEm1d|!D z?n*tcft$WKdw8_Qb0%z&5eVAEzc2BpDVw*)OkDf+;@qbA)ZZcJ#6V3RE?0)O)#)uU zi0$0#MZzg}B!EzP$a&hAP_rg#Q0bKAX!q*GrR;4bFZyiuu09Wl`6VnbJJof1LX3{N zkOZPff1{Bs;!%NhW814`#G$e6bi{p!FvN3PZ}Cv6N(eXEs3$1J3*&)<90#GeH?s@QlH}3(8HL0ZZr%c;{%L3vp#$P{iW| zo`_!GMwEhNQ8W_4M9eqnUulS?>lT>@jRI}jOod|UcOd(mEgNB*=@AXq;p;K26lBv| zZe~x!f(r1Q*?{PrAp(?rc#O2Zfkn;Zsr}W^GHTbI{In$3zxkAaPBpmlJuvyf#(?4P zMtT3`W`HUtTifPM8;Gw4`aqMeitXrfWZh<3IDT5$+nd)4bAF~m;lwn+bA6j2)zC|W zZxH`&T7Y*TKHUJx$@AJNPG(s3Fq3{Y-CVt%l)nNKLY#UvwRLoofu*tF(9%X!UxVlr z7poTZ8nf_obBu(lU169wlvp2Pb@kOVASN2}Oq!_?U0(dt#Vb{;cwKdG9ka-j( z9p2!_CqR{&G6eF*h;Ktf#Sr5C&^`E8;T5wB5s{Ie;E0%Qy`J3spp~bOZ*FTV=Cyzn z$c$((KJhr0mfKEMHMMV6z!&Yj74JHV&&tfq6fT1e>;!r<%XpFTKt5{&A#gNi9xM&} z?<-agbg0`Z#n~7Zc{ee6=4t}wty|j=rCDj#)HW~+cLyoiN?41;>u70B;zw($R(?77 zZs;ikQDJ-?8fG3b7AB-uW#Y)04N z2PI5=I?s@}fM!Os*J+5phtsWt1n71eSw5?odh!w7ZpDy4-SHJib1+DD$5k*`! z{#?;wsLohVpY+FE=4l7UuJmQun{B$OYih0Aq-6iU}=vK7~Fu}xj>^xalRLd?hGM6X^xbx}Y_=jWoGt-8^h z_BbJbRX)3oFZG!N=l^CddhmvNXC9<`osA*MOfdY~{Zy!RIZaLYzLa>6zHkJZ|K1OkKe z{P?)9t;7RQc;tI>1+1|8Wkv6twGW%)QVxdY*9h=Qjnwxgz$}&c2aR+%Vhmf=LjYh} zDqW`LF}&t>ibX{=`GiN^>`Kep>P9w8{(d?et8n(P4NniC6%3DF6c|W3EV21sKoA}x zhTO31_sYO!zznzdHrr)>h-U45ocS!;#cLs5vNe(?GPKRf+4lp_Qu$XKBT&%jd>E^O zV-r@|&87IhTS#iycG~g&iPY$MyF#gUjk%k4a;CQvquN7Y(75T_82lua@`RA@Va17= zG08^*Y>#o5cj)+T_1=~u!8V>syV2%o^_F}LvGUf~<)n8ytILo{bHi4ud*|qyz zd*Kdq#m>jk0tymPcbZhvE}EgUw_X9C0Nl9f?KY@^gs7Q`x9`tpb}BpXZb@Qld3it4 zn)lZDxXsk%CAQ&S9G;pyb_e8GpYt{^=xoAflSjlmDFw(6CFC_Nl^ z92yn$GS>IjHLje7#W~9M@enwu6AO2>61;Z9XYjE$NvoPAv?@gms#We_&Y#JnJNdz6 zz}g;WxhXjIEUmGda^sE1gQB?NYB|4zegm_6La)wmFBg^zGN4azosl9s)m?oQ6Kx=WW3?bZ;X0va}y)3MrcCa-NG6<**qB+ABDsEG~JLv6c}kjAkij^}PC zRb#Hxb+a^oToO8}ssJh$vC*xiVi|%GS+#6yr5VoOk$mY`K#Q5rjB!(6U?m8h1!Me% zOMy^XJ(3PICyUM`(6Wyq11;+hduT_?NZaD@{y>+!uK3FFfu{PQS@{eD_qh*hj>{A7 zKt`$yc?VdvL?k;mC#SbMx?Mf+YRTI~AnQg@nzNc;lE0j+JXv5kVe(+k6MbpCM`TYt z@I(lSA2Y%@9RjJSde5;u?+likh`aRZ?ThsfNPb1l)ffjkWcW4gb4ZwSJl3kcvNvdR z*@e!wrjI=3DZ{M*kOY>^94y%J&Dq&mlNcHm+_QEco@&d)k=GWlaz%{zHa1%jbk@3X z4^RpjnuG*Z5kFFD6I)!sTyPmTEo~4gQ>RFXBdIlM7Zi#48#AXefX=g12s}" + for name in methods: + header += f" {name:>15s}" + print(header) + print("-" * len(header)) + for i, npart in enumerate(particle_counts): + row = f"{npart:>12,d}" + for name in methods: + row += f" {results[name][i]:>14.3f}s" + print(row) + + # --- Plot --- + fig, ax = plt.subplots(figsize=(9, 6)) + markers = ["o", "s", "^", "D"] + for j, (name, times) in enumerate(results.items()): + ax.loglog( + particle_counts, + times, + marker=markers[j % len(markers)], + linewidth=2, + markersize=7, + label=name, + ) + + ax.set_xlabel("Number of particles") + ax.set_ylabel("Wall-clock time (s)") + ax.set_title("Parcels simulation scaling: windowed vs cached chunk arrays") + ax.legend() + ax.grid(True, which="both", alpha=0.3) + fig.tight_layout() + fig.savefig("benchmark_chunk_cache.png", dpi=150) + print("\nPlot saved to benchmark_chunk_cache.png") + plt.show() if __name__ == "__main__": From 6c60fcae00cbc84ab978599937baf23c83263217 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 31 Aug 2026 13:09:55 +0800 Subject: [PATCH 05/12] Start working on test_backends.py --- tests/test_backends.py | 83 ++++++++++++++++++++++++++++++++++++++++++ tests/test_fieldset.py | 2 + 2 files changed, 85 insertions(+) create mode 100644 tests/test_backends.py diff --git a/tests/test_backends.py b/tests/test_backends.py new file mode 100644 index 000000000..30d307ac3 --- /dev/null +++ b/tests/test_backends.py @@ -0,0 +1,83 @@ +import io +from pathlib import Path +from typing import Literal + +import pytest +import xarray as xr + +import parcels +import parcels.tutorial + +BackendT = Literal["WindowedArray", "Dask", "Zarr", "NumPy", "CachedChunkArray"] +BACKENDS = {"WindowedArray", "Dask", "Zarr", "NumPy", "CachedChunkArray"} + + +@pytest.fixture(scope="module") +def nemo_dataset() -> xr.Dataset: + ds_u = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/U") + ds_v = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/V") + ds_w = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/W") + ds_coords = parcels.tutorial.open_dataset("NemoNorthSeaORCA025-N006_data/mesh_mask")[["glamf", "gphif"]] + + ds_fset = parcels.convert.nemo_to_sgrid( + fields={"U": ds_u["uo"], "V": ds_v["vo"], "W": ds_w["wo"]}, + coords=ds_coords, + ) + return ds_fset + + +@pytest.fixture(scope="module") +def nemo_results(tmp_parquet, nemo_dataset) -> tuple[xr.Dataset, Path]: + run_simulation(nemo_dataset, tmp_parquet, "NumPy") + return nemo_dataset, tmp_parquet + + +def assert_fieldset_backend(fset: parcels.FieldSet, backend: BackendT): + # a bit of a hacky way to check for the backend.... probably better for us to change how backends are stored + buf = io.StringIO() + fset.describe(buf) + return backend in buf.getvalue() + + +def run_simulation(ds: xr.Dataset, output_path: Path, backend: BackendT) -> Path: + if backend == "Zarr": + raise NotImplementedError("Doesn't work at this level of execution. Also will likely remove Zarr backend.") + + if backend == "NumPy": + ds.load() + + fset = parcels.FieldSet.from_sgrid_conventions(ds) + + if backend == "WindowedArray": + fset.to_windowed_arrays() + if backend == "CachedChunkArray": + fset.to_cached_chunk_arrays() + + assert_fieldset_backend(fset, backend) + + # TODO Create the particleset by seeding 1000 particles + ... + + # TODO Advect the particles using a RK4 advection kernel, saving the results to the Parquet file + ... + + return output_path + + +@pytest.mark.parametrize( + "backend", + BACKENDS + - { + "NumPY", # reference point + "Zarr", # not supported + }, +) +def test_nemo_identical_across_backends(nemo_results, tmp_parquet, backend): + ds = nemo_results[0] + ref_parquet = nemo_results[1] + + assert str(ref_parquet) != str(tmp_parquet) # just covering my bases with Pytest fixture usage + + run_simulation(ds, tmp_parquet, backend) + + # TODO: Compare the tmp_parquet with the reference parquet diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 75ca0868a..e575a425b 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -542,3 +542,5 @@ def test_fieldset_describe_backends(tmp_path): fieldset.describe(io) actual = io.getvalue() assert actual == expected + + # TODO: Add test for the ChunkedArray backend (can also refactor this test at the same time) From 1864cc6953801a7c01c876d4ef74641980fcf17f Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 31 Aug 2026 14:06:22 +0800 Subject: [PATCH 06/12] Write test_nemo_identical_across_backends --- tests/test_backends.py | 52 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 42 insertions(+), 10 deletions(-) diff --git a/tests/test_backends.py b/tests/test_backends.py index 30d307ac3..895aaadc3 100644 --- a/tests/test_backends.py +++ b/tests/test_backends.py @@ -2,11 +2,14 @@ from pathlib import Path from typing import Literal +import numpy as np +import pandas as pd import pytest import xarray as xr import parcels import parcels.tutorial +from parcels.kernels import AdvectionRK4 BackendT = Literal["WindowedArray", "Dask", "Zarr", "NumPy", "CachedChunkArray"] BACKENDS = {"WindowedArray", "Dask", "Zarr", "NumPy", "CachedChunkArray"} @@ -27,9 +30,10 @@ def nemo_dataset() -> xr.Dataset: @pytest.fixture(scope="module") -def nemo_results(tmp_parquet, nemo_dataset) -> tuple[xr.Dataset, Path]: - run_simulation(nemo_dataset, tmp_parquet, "NumPy") - return nemo_dataset, tmp_parquet +def nemo_results(tmp_path_factory, nemo_dataset) -> tuple[xr.Dataset, Path]: + ref_parquet = tmp_path_factory.mktemp("nemo_ref") / "ref.parquet" + run_simulation(nemo_dataset, ref_parquet, "NumPy") + return nemo_dataset, ref_parquet def assert_fieldset_backend(fset: parcels.FieldSet, backend: BackendT): @@ -55,11 +59,31 @@ def run_simulation(ds: xr.Dataset, output_path: Path, backend: BackendT) -> Path assert_fieldset_backend(fset, backend) - # TODO Create the particleset by seeding 1000 particles - ... - - # TODO Advect the particles using a RK4 advection kernel, saving the results to the Parquet file - ... + npart = 1000 + lons = np.linspace(1.9, 3.4, npart) + lats = np.linspace(51.6, 52.5, npart) + z = np.ones(npart) + pset = parcels.ParticleSet(fset, x=lons, y=lats, z=z) + + def delete_particle(particles, fieldset): + error_states = ( + parcels.StatusCode.ErrorOutOfBounds, + parcels.StatusCode.ErrorGridSearching, + ) + for error in error_states: + particles.state = np.where( + particles.state == error, + parcels.StatusCode.Delete, + particles.state, + ) + + pfile = parcels.ParticleFile(output_path, outputdt=np.timedelta64(6, "h")) + pset.execute( + [AdvectionRK4, delete_particle], + runtime=np.timedelta64(3, "D"), + dt=np.timedelta64(5, "m"), + output_file=pfile, + ) return output_path @@ -68,7 +92,7 @@ def run_simulation(ds: xr.Dataset, output_path: Path, backend: BackendT) -> Path "backend", BACKENDS - { - "NumPY", # reference point + "NumPy", # reference point "Zarr", # not supported }, ) @@ -80,4 +104,12 @@ def test_nemo_identical_across_backends(nemo_results, tmp_parquet, backend): run_simulation(ds, tmp_parquet, backend) - # TODO: Compare the tmp_parquet with the reference parquet + ref_df = pd.read_parquet(ref_parquet) + test_df = pd.read_parquet(tmp_parquet) + + ref_df = ref_df.sort_values(["particle_id", "t"]).reset_index(drop=True) + test_df = test_df.sort_values(["particle_id", "t"]).reset_index(drop=True) + + np.testing.assert_allclose(test_df["x"].values, ref_df["x"].values, atol=1e-5) + np.testing.assert_allclose(test_df["y"].values, ref_df["y"].values, atol=1e-5) + np.testing.assert_allclose(test_df["z"].values, ref_df["z"].values, atol=1e-5) From 414bf7f190916441d88e7ded174ed5b45f4a6b5c Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 31 Aug 2026 16:34:38 +0800 Subject: [PATCH 07/12] Migrate implementation of ChunkCachedArray Moves the code to a Parcels subpackage (rather than a separate package). --- pixi.toml | 1 - src/parcels/_chunk_cached_array/__init__.py | 4 + src/parcels/_chunk_cached_array/core.py | 188 +++++++++++++++++++ src/parcels/_chunk_cached_array/lru_cache.py | 48 +++++ src/parcels/_chunk_cached_array/py.typed | 0 src/parcels/_core/fieldset.py | 2 +- src/parcels/_core/model.py | 4 +- 7 files changed, 243 insertions(+), 4 deletions(-) create mode 100644 src/parcels/_chunk_cached_array/__init__.py create mode 100644 src/parcels/_chunk_cached_array/core.py create mode 100644 src/parcels/_chunk_cached_array/lru_cache.py create mode 100644 src/parcels/_chunk_cached_array/py.typed diff --git a/pixi.toml b/pixi.toml index 79a4b96b2..4b699435b 100644 --- a/pixi.toml +++ b/pixi.toml @@ -41,7 +41,6 @@ tabulate = ">=0.10.0" [dependencies] -chunk_cached_array = { path = "../xarray-interpolation/chunk-cached-array", package.build.backend.name = "pixi-build-python" } parcels = { path = "." } [feature.rattler-build.dependencies] diff --git a/src/parcels/_chunk_cached_array/__init__.py b/src/parcels/_chunk_cached_array/__init__.py new file mode 100644 index 000000000..89920eefe --- /dev/null +++ b/src/parcels/_chunk_cached_array/__init__.py @@ -0,0 +1,4 @@ +from .core import ChunkCachedArray, wrap_dataset +from .lru_cache import ByteBoundedLRUCache + +__all__ = ["ByteBoundedLRUCache", "ChunkCachedArray", "wrap_dataset"] diff --git a/src/parcels/_chunk_cached_array/core.py b/src/parcels/_chunk_cached_array/core.py new file mode 100644 index 000000000..902d396d4 --- /dev/null +++ b/src/parcels/_chunk_cached_array/core.py @@ -0,0 +1,188 @@ +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np +from xarray.core.indexing import BasicIndexer, ExplicitlyIndexedNDArrayMixin, OuterIndexer, VectorizedIndexer +from xarray.namedarray.pycompat import is_duck_array + +from .lru_cache import ByteBoundedLRUCache + +if TYPE_CHECKING: + import dask.array + import xarray as xr + + +def wrap_dataset(ds: xr.Dataset, max_cache_bytes: int) -> xr.Dataset: + """Replace all dask-backed data variables with ChunkCachedArray wrappers. + + Returns a shallow copy of the dataset. Each dask-backed data variable's + internal ``variable._data`` is swapped for a ``ChunkCachedArray`` that + caches chunks on vectorized indexing. Coordinate variables are loaded + eagerly into memory to avoid dask task-graph overhead on every + ``.isel()`` call. + + Parameters + ---------- + ds : xr.Dataset + Source dataset (not modified). + max_cache_bytes : int + Maximum cache size in bytes, per variable. + + Returns + ------- + xr.Dataset + Copy with dask arrays wrapped in ChunkCachedArray. + """ + from dask.base import is_dask_collection + + ds = ds.copy() + # Load coordinates eagerly — they are small 1D arrays and keeping them + # as dask arrays causes expensive task-graph construction on every .isel(). + for name in list(ds.coords): + ds[name].load() + for name in ds.data_vars: + var = ds[name].variable + if is_duck_array(var._data) and is_dask_collection(var._data): + var._data = ChunkCachedArray(var._data, max_cache_bytes) + return ds + + +class ChunkCachedArray(ExplicitlyIndexedNDArrayMixin): + """Chunk-level LRU cache on top of a dask array for vectorized indexing. + + Implements xarray's ExplicitlyIndexed protocol so it can be used as + a drop-in replacement for the dask array in ``da.data``. Xarray's + ``.isel()`` with vectorized indexers will route through ``_vindex_get``, + which uses the chunk cache. Other indexing modes delegate to the + underlying dask array. + + On each vectorized index: + 1. Maps global indices -> (chunk_coord, local_index) per dimension. + 2. Fetches missing chunks via dask_array.blocks[...].compute(). + 3. Assembles the result from cached numpy arrays. + """ + + def __init__(self, dask_array: dask.array.Array, max_cache_bytes: int) -> None: + self.array = dask_array + self.cache = ByteBoundedLRUCache(max_cache_bytes) + + # Precompute chunk boundaries per dimension. + # _boundaries[d] is a 1D array of cumulative chunk sizes, e.g., [0, 15, 30]. + self._boundaries: list[np.ndarray] = [] + for dim_chunks in dask_array.chunks: + self._boundaries.append(np.concatenate(([0], np.cumsum(dim_chunks)))) + + def get_duck_array(self): + return self.array.compute() + + def _raw_vindex(self, *indices: np.ndarray) -> np.ndarray: + """Vectorized indexing with chunk caching. + + Parameters + ---------- + *indices : np.ndarray + One 1D integer index array per dimension. All must have the same length N. + + Returns + ------- + np.ndarray + 1D array of length N with the selected values. + """ + ndim = len(self.array.chunks) + assert len(indices) == ndim + n_points = len(indices[0]) + + # Step 1: Map global indices to chunk coords and local indices. + # Normalize negative indices (e.g. -1 → last element) to positive, + # matching standard numpy fancy-indexing semantics. + indices = tuple(np.where(idx < 0, idx + self.array.shape[d], idx) for d, idx in enumerate(indices)) + chunk_ids = np.empty((ndim, n_points), dtype=np.intp) + local_indices = np.empty((ndim, n_points), dtype=np.intp) + for d in range(ndim): + cid = np.searchsorted(self._boundaries[d], indices[d], side="right") - 1 + chunk_ids[d] = cid + local_indices[d] = indices[d] - self._boundaries[d][cid] + + # Step 2: Group points by chunk using a structured array for vectorized grouping. + # Encode each point's chunk coords as a single int for fast grouping. + # Use np.ravel_multi_index on chunk_ids to get a flat chunk key per point. + numblocks = np.array(self.array.numblocks, dtype=np.intp) + flat_keys = np.ravel_multi_index(chunk_ids, numblocks) + + # Sort points by flat chunk key to group them. + sort_order = np.argsort(flat_keys, kind="mergesort") + sorted_flat_keys = flat_keys[sort_order] + + # Find group boundaries. + boundaries = np.concatenate(([0], np.flatnonzero(np.diff(sorted_flat_keys)) + 1, [n_points])) + + out = np.empty(n_points, dtype=self.array.dtype) + for g in range(len(boundaries) - 1): + grp_slice = slice(boundaries[g], boundaries[g + 1]) + grp_indices = sort_order[grp_slice] + + # Recover the chunk key tuple from any point in this group. + key = tuple(int(chunk_ids[d, grp_indices[0]]) for d in range(ndim)) + + chunk_data = self.cache.get(key) + if chunk_data is None: + chunk_data = self.array.blocks[key].compute() + self.cache.put(key, chunk_data) + + # Vectorized fancy-index: extract all points from this chunk at once. + local_idx = tuple(local_indices[d, grp_indices] for d in range(ndim)) + out[grp_indices] = chunk_data[local_idx] + + return out + + # --- ExplicitlyIndexed protocol --- + + def _vindex_get(self, indexer: VectorizedIndexer): + key = indexer.tuple + return self._raw_vindex(*key) + + def _oindex_get(self, indexer: OuterIndexer): + # Delegate to dask for orthogonal indexing + return self.array[indexer.tuple] + + def __getitem__(self, indexer): + if isinstance(indexer, VectorizedIndexer): + return self._vindex_get(indexer) + if isinstance(indexer, OuterIndexer): + return self._oindex_get(indexer) + if isinstance(indexer, BasicIndexer): + return self.array[indexer.tuple] + return self.array[indexer] + + # --- Convenience methods (direct use without xarray) --- + # Note: .vindex is inherited from ExplicitlyIndexedNDArrayMixin as a property + # returning IndexCallable(self._vindex_get). Do NOT shadow it with a method. + + def isel( + self, + indexers: dict[str, object], + dims: tuple[str, ...], + ) -> np.ndarray: + """Dimension-name-keyed vectorized indexing (xarray .isel() style). + + Parameters + ---------- + indexers : dict[str, array-like] + Mapping of dimension name to index array. Values can be np.ndarray + or xr.DataArray (in which case .values is extracted). + dims : tuple[str, ...] + Ordered dimension names of the underlying array. + + Returns + ------- + np.ndarray + 1D array of selected values. + """ + arrays = [] + for dim in dims: + idx = indexers[dim] + if hasattr(idx, "values"): + idx = idx.values + arrays.append(np.asarray(idx)) + return self._raw_vindex(*arrays) diff --git a/src/parcels/_chunk_cached_array/lru_cache.py b/src/parcels/_chunk_cached_array/lru_cache.py new file mode 100644 index 000000000..117593439 --- /dev/null +++ b/src/parcels/_chunk_cached_array/lru_cache.py @@ -0,0 +1,48 @@ +from __future__ import annotations + +from collections import OrderedDict +from collections.abc import Hashable + +import numpy as np + + +class ByteBoundedLRUCache: + """LRU cache bounded by total stored bytes. + + Keys are hashable (typically chunk coordinate tuples). + Values are numpy arrays whose .nbytes drives eviction. + """ + + def __init__(self, max_bytes: int) -> None: + assert max_bytes > 0 + self._max_bytes = max_bytes + self._cache: OrderedDict[Hashable, np.ndarray] = OrderedDict() + self._current_bytes = 0 + + @property + def current_bytes(self) -> int: + return self._current_bytes + + def get(self, key: Hashable) -> np.ndarray | None: + if key not in self._cache: + return None + self._cache.move_to_end(key) + return self._cache[key] + + def put(self, key: Hashable, value: np.ndarray) -> None: + nbytes = value.nbytes + if nbytes > self._max_bytes: + return + # Remove existing entry if present (will be re-inserted at end) + if key in self._cache: + self._current_bytes -= self._cache.pop(key).nbytes + # Evict LRU entries until there's room + while self._current_bytes + nbytes > self._max_bytes: + _, evicted = self._cache.popitem(last=False) + self._current_bytes -= evicted.nbytes + self._cache[key] = value + self._current_bytes += nbytes + + def clear(self) -> None: + self._cache.clear() + self._current_bytes = 0 diff --git a/src/parcels/_chunk_cached_array/py.typed b/src/parcels/_chunk_cached_array/py.typed new file mode 100644 index 000000000..e69de29bb diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index 6b3d53371..ccb945eb2 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -176,7 +176,7 @@ def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000): """Wrap dask-backed field data in chunk-level LRU caches. Opt-in optimization that replaces each dask-backed data variable's - internal storage with a :class:`~chunk_cached_array.ChunkCachedArray`. + internal storage with a :class:`~parcels._chunk_cached_array.ChunkCachedArray`. Delegates to each underlying model; repeated vectorized ``.isel()`` calls then hit an in-memory LRU cache instead of recomputing dask task graphs. NumPy-backed (eager) fields are left unchanged, and re-invoking diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index e188b996b..0e3338c99 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -9,11 +9,11 @@ import uxarray as ux import xarray as xr import zarr -from chunk_cached_array import wrap_dataset from dask import is_dask_collection import parcels._sgrid as sgrid import parcels._typing as ptyping +from parcels._chunk_cached_array import wrap_dataset from parcels._core._windowed_array import maybe_windowed from parcels._core.basegrid import BaseGrid from parcels._core.field import Field, VectorField @@ -117,7 +117,7 @@ def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000) -> Self: """Wrap dask-backed field data in chunk-level LRU caches. Opt-in optimization that replaces each dask-backed data variable's - internal storage with a :class:`~chunk_cached_array.ChunkCachedArray`. + internal storage with a :class:`~parcels._chunk_cached_array.ChunkCachedArray`. Repeated vectorized ``.isel()`` calls then hit an in-memory LRU cache keyed by chunk coordinates instead of recomputing dask task graphs. From 7a9b822d3b5313da424041acf6b9175988855049 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Mon, 31 Aug 2026 16:39:00 +0800 Subject: [PATCH 08/12] Remove isel() method This was a testing-only convenience method, now handled internally by Xarray --- src/parcels/_chunk_cached_array/core.py | 32 ------------------------- 1 file changed, 32 deletions(-) diff --git a/src/parcels/_chunk_cached_array/core.py b/src/parcels/_chunk_cached_array/core.py index 902d396d4..53775a082 100644 --- a/src/parcels/_chunk_cached_array/core.py +++ b/src/parcels/_chunk_cached_array/core.py @@ -154,35 +154,3 @@ def __getitem__(self, indexer): if isinstance(indexer, BasicIndexer): return self.array[indexer.tuple] return self.array[indexer] - - # --- Convenience methods (direct use without xarray) --- - # Note: .vindex is inherited from ExplicitlyIndexedNDArrayMixin as a property - # returning IndexCallable(self._vindex_get). Do NOT shadow it with a method. - - def isel( - self, - indexers: dict[str, object], - dims: tuple[str, ...], - ) -> np.ndarray: - """Dimension-name-keyed vectorized indexing (xarray .isel() style). - - Parameters - ---------- - indexers : dict[str, array-like] - Mapping of dimension name to index array. Values can be np.ndarray - or xr.DataArray (in which case .values is extracted). - dims : tuple[str, ...] - Ordered dimension names of the underlying array. - - Returns - ------- - np.ndarray - 1D array of selected values. - """ - arrays = [] - for dim in dims: - idx = indexers[dim] - if hasattr(idx, "values"): - idx = idx.values - arrays.append(np.asarray(idx)) - return self._raw_vindex(*arrays) From 9d9fbbb666ea6a4939c27767314a05cd3a8a70a2 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 1 Sep 2026 12:15:46 +0800 Subject: [PATCH 09/12] Consolidate naming to `chunk_cached` instead of `cached_chunk` --- benchmark_chunk_cache.py | 2 +- src/parcels/_core/fieldset.py | 4 ++-- src/parcels/_core/model.py | 2 +- tests/test_backends.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/benchmark_chunk_cache.py b/benchmark_chunk_cache.py index fa457079d..e0e13573b 100644 --- a/benchmark_chunk_cache.py +++ b/benchmark_chunk_cache.py @@ -83,7 +83,7 @@ def main(): # } methods = { "windowed": lambda ds: make_fieldset(ds).to_windowed_arrays(), - "cached chunks": lambda ds: make_fieldset(ds).to_cached_chunk_arrays(), + "cached chunks": lambda ds: make_fieldset(ds).to_chunk_cached_arrays(), } results = {name: [] for name in methods} diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py index ccb945eb2..7f6d65ca7 100644 --- a/src/parcels/_core/fieldset.py +++ b/src/parcels/_core/fieldset.py @@ -172,7 +172,7 @@ def to_windowed_arrays(self, *, max_levels: int | None = None): model.to_windowed_arrays(max_levels=max_levels) return self - def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000): + def to_chunk_cached_arrays(self, *, max_cache_bytes: int = 600_000_000): """Wrap dask-backed field data in chunk-level LRU caches. Opt-in optimization that replaces each dask-backed data variable's @@ -193,7 +193,7 @@ def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000): ``self``, to allow chaining. """ for model in self.models: - model.to_cached_chunk_arrays(max_cache_bytes=max_cache_bytes) + model.to_chunk_cached_arrays(max_cache_bytes=max_cache_bytes) return self def add_constant_field(self, name: str, value, mesh: ptyping.TMesh = "spherical"): diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index 0e3338c99..255e7b056 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -113,7 +113,7 @@ def to_windowed_arrays(self, *, max_levels: int | None = None) -> Self: windowed[name] = maybe_windowed(current, max_levels=max_levels) return self - def to_cached_chunk_arrays(self, *, max_cache_bytes: int = 600_000_000) -> Self: + def to_chunk_cached_arrays(self, *, max_cache_bytes: int = 600_000_000) -> Self: """Wrap dask-backed field data in chunk-level LRU caches. Opt-in optimization that replaces each dask-backed data variable's diff --git a/tests/test_backends.py b/tests/test_backends.py index 895aaadc3..31286cc7f 100644 --- a/tests/test_backends.py +++ b/tests/test_backends.py @@ -55,7 +55,7 @@ def run_simulation(ds: xr.Dataset, output_path: Path, backend: BackendT) -> Path if backend == "WindowedArray": fset.to_windowed_arrays() if backend == "CachedChunkArray": - fset.to_cached_chunk_arrays() + fset.to_chunk_cached_arrays() assert_fieldset_backend(fset, backend) From b581629407b8685ce1f9f338396501df193c957d Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 1 Sep 2026 12:27:00 +0800 Subject: [PATCH 10/12] Update benchmarking script and add data-generation --- benchmark_chunk_cache.py | 2 +- data-generation.py | 294 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 295 insertions(+), 1 deletion(-) create mode 100644 data-generation.py diff --git a/benchmark_chunk_cache.py b/benchmark_chunk_cache.py index e0e13573b..d21047e22 100644 --- a/benchmark_chunk_cache.py +++ b/benchmark_chunk_cache.py @@ -70,7 +70,7 @@ def run_simulation(fieldset, ds, npart): def main(): - zarr_path = "../xarray-interpolation/datasets/ds_2d_left_agrid.zarr" + zarr_path = "./datasets/ds_2d_left_agrid.zarr" particle_counts = [10, 100, 1_000, 10_000, 100_000, 1_000_000] print(f"Loading dataset from {zarr_path}") diff --git a/data-generation.py b/data-generation.py new file mode 100644 index 000000000..28c09a48e --- /dev/null +++ b/data-generation.py @@ -0,0 +1,294 @@ +import operator +from functools import reduce +from pathlib import Path + +import dask.array as da +import numpy as np +import xarray as xr +from zarr.storage import LocalStore + +X = 1000 +Y = 1000 +Z = 50 +T = 30 + +X_SMALL = 200 +Y_SMALL = 200 + +TIME = xr.date_range("2000", "2001", T) + + +def prod(seq): + return reduce(operator.mul, seq, 1) + + +def _rotated_curvilinear_grid(): + XG = np.arange(X) + YG = np.arange(Y) + LON, LAT = np.meshgrid(XG, YG) + + angle = -np.pi / 24 + rotation = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) + + # rotate the LON and LAT grids + LON, LAT = np.einsum("ji, mni -> jmn", rotation, np.dstack([LON, LAT])) + + return xr.Dataset( + { + "data_g": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "data_c": ( + ["time", "ZC", "YC", "XC"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "U_A_grid": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "V_A_grid": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "U_C_grid": ( + ["time", "ZG", "YC", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "V_C_grid": ( + ["time", "ZG", "YG", "XC"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + }, + coords={ + "XG": (["XG"], XG, {"axis": "X", "c_grid_axis_shift": -0.5}), + "YG": (["YG"], YG, {"axis": "Y", "c_grid_axis_shift": -0.5}), + "XC": (["XC"], XG + 0.5, {"axis": "X"}), + "YC": (["YC"], YG + 0.5, {"axis": "Y"}), + "ZG": ( + ["ZG"], + np.arange(Z), + {"axis": "Z", "c_grid_axis_shift": -0.5}, + ), + "ZC": ( + ["ZC"], + np.arange(Z) + 0.5, + {"axis": "Z"}, + ), + "depth": (["ZG"], np.arange(Z), {"axis": "Z"}), + "time": (["time"], TIME, {"axis": "T"}), + "lon": ( + ["YG", "XG"], + LON, + {"axis": "X", "c_grid_axis_shift": -0.5}, # ? Needed? + ), + "lat": ( + ["YG", "XG"], + LAT, + {"axis": "Y", "c_grid_axis_shift": -0.5}, # ? Needed? + ), + }, + ) + + +def random_dask_array(shape, scaling=1): + return da.random.uniform(scaling, size=prod(shape)).reshape(shape) + + +def _cartesion_to_polar(x, y): + r = np.sqrt(x**2 + y**2) + theta = np.arctan2(y, x) + return r, theta + + +def _polar_to_cartesian(r, theta): + x = r * np.cos(theta) + y = r * np.sin(theta) + return x, y + + +def _unrolled_cone_curvilinear_grid(): + # Not a great unrolled cone, but this is good enough for testing + # you can use matplotlib pcolormesh to plot + XG = np.arange(X) + YG = np.arange(Y) * 0.25 + + pivot = -10, 0 + LON, LAT = np.meshgrid(XG, YG) + + new_lon_lat = [] + + min_lon = np.min(XG) + for lon, lat in zip(LON.flatten(), LAT.flatten(), strict=True): + r, _ = _cartesion_to_polar(lon - pivot[0], lat - pivot[1]) + _, theta = _cartesion_to_polar(min_lon - pivot[0], lat - pivot[1]) + theta *= 1.2 + r *= 1.2 + lon, lat = _polar_to_cartesian(r, theta) + new_lon_lat.append((lon + pivot[0], lat + pivot[1])) + + new_lon, new_lat = zip(*new_lon_lat, strict=True) + LON, LAT = ( + np.array(new_lon).reshape(LON.shape), + np.array(new_lat).reshape(LAT.shape), + ) + + return xr.Dataset( + { + "data_g": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "data_c": ( + ["time", "ZC", "YC", "XC"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "U_A_grid": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "V_A_grid": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "U_C_grid": ( + ["time", "ZG", "YC", "XG"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + "V_C_grid": ( + ["time", "ZG", "YG", "XC"], + random_dask_array(scaling=5, shape=(T, Z, Y, X)), + ), + }, + coords={ + "XG": (["XG"], XG, {"axis": "X", "c_grid_axis_shift": -0.5}), + "YG": (["YG"], YG, {"axis": "Y", "c_grid_axis_shift": -0.5}), + "XC": (["XC"], XG + 0.5, {"axis": "X"}), + "YC": (["YC"], YG + 0.5, {"axis": "Y"}), + "ZG": ( + ["ZG"], + np.arange(Z), + {"axis": "Z", "c_grid_axis_shift": -0.5}, + ), + "ZC": ( + ["ZC"], + np.arange(Z) + 0.5, + {"axis": "Z"}, + ), + "depth": (["ZG"], np.arange(Z), {"axis": "Z"}), + "time": (["time"], TIME, {"axis": "T"}), + "lon": ( + ["YG", "XG"], + LON, + {"axis": "X", "c_grid_axis_shift": -0.5}, # ? Needed? + ), + "lat": ( + ["YG", "XG"], + LAT, + {"axis": "Y", "c_grid_axis_shift": -0.5}, # ? Needed? + ), + }, + ) + + +def _ds_2d_left(x, y, z, t, time): + """MITgcm indexing style dataset.""" + return xr.Dataset( + { + "data_g": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(t, z, y, x)), + ), + "data_c": ( + ["time", "ZC", "YC", "XC"], + random_dask_array(scaling=5, shape=(t, z, y, x)), + ), + "U_A_grid": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(t, z, y, x)), + ), + "V_A_grid": ( + ["time", "ZG", "YG", "XG"], + random_dask_array(scaling=5, shape=(t, z, y, x)), + ), + "U_C_grid": ( + ["time", "ZG", "YC", "XG"], + random_dask_array(scaling=5, shape=(t, z, y, x)), + ), + "V_C_grid": ( + ["time", "ZG", "YG", "XC"], + random_dask_array(scaling=5, shape=(t, z, y, x)), + ), + }, + coords={ + "XG": ( + ["XG"], + 2 * np.pi / x * np.arange(0, x), + {"axis": "X", "c_grid_axis_shift": -0.5}, + ), + "XC": (["XC"], 2 * np.pi / x * (np.arange(0, x) + 0.5), {"axis": "X"}), + "YG": ( + ["YG"], + 2 * np.pi / y * np.arange(0, y), + {"axis": "Y", "c_grid_axis_shift": -0.5}, + ), + "YC": ( + ["YC"], + 2 * np.pi / y * (np.arange(0, y) + 0.5), + {"axis": "Y"}, + ), + "ZG": ( + ["ZG"], + np.arange(z), + {"axis": "Z", "c_grid_axis_shift": -0.5}, + ), + "ZC": ( + ["ZC"], + np.arange(z) + 0.5, + {"axis": "Z"}, + ), + "lon": (["XG"], 2 * np.pi / x * np.arange(0, x)), + "lat": (["YG"], 2 * np.pi / y * np.arange(0, y)), + "depth": (["ZG"], np.arange(z)), + "time": (["time"], time, {"axis": "T"}), + }, + ) + + +datasets = { + "2d_left_rotated": _rotated_curvilinear_grid(), + "ds_2d_left": _ds_2d_left(X, Y, Z, T, TIME), + "2d_left_unrolled_cone": _unrolled_cone_curvilinear_grid(), +} + + +def save(ds: xr.Dataset, path: str, chunks: dict) -> None: + """Save dataset to zarr with specified chunking.""" + store = LocalStore(path) + ds.chunk(chunks).to_zarr(store, mode="w", encoding=None, consolidated=False) + size_mb = sum(ds[v].nbytes for v in ds.data_vars) / 1e6 + print(f" {path} dims={dict(ds.sizes)} ~{size_mb:.0f} MB uncompressed") + + +if __name__ == "__main__": + dataset_path = "datasets/ds_2d_left_agrid.zarr" + print("Generating ds_2d_left...") + if Path(dataset_path).exists(): + print(f"Dataset {dataset_path} already exists") + else: + save( + datasets["ds_2d_left"][["U_A_grid", "V_A_grid"]], + dataset_path, + {"time": 15, "XG": 40, "YG": 40, "ZG": 8}, + ) + + dataset_path_small = "datasets/ds_2d_left_agrid_small.zarr" + print("Generating ds_2d_left_small...") + if Path(dataset_path_small).exists(): + print(f"Dataset {dataset_path_small} already exists") + else: + save( + _ds_2d_left(X_SMALL, Y_SMALL, Z, T, TIME)[["U_A_grid", "V_A_grid"]], + dataset_path_small, + {"time": 15, "XG": 40, "YG": 40, "ZG": 8}, + ) From aa3c98a6b13c1128a0733aa398984b644421cf17 Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 1 Sep 2026 12:35:55 +0800 Subject: [PATCH 11/12] remove py.typed marker --- src/parcels/_chunk_cached_array/py.typed | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 src/parcels/_chunk_cached_array/py.typed diff --git a/src/parcels/_chunk_cached_array/py.typed b/src/parcels/_chunk_cached_array/py.typed deleted file mode 100644 index e69de29bb..000000000 From 35b5324ef5d9214f47fc8c3669a56410ede001ba Mon Sep 17 00:00:00 2001 From: Vecko <36369090+VeckoTheGecko@users.noreply.github.com> Date: Tue, 1 Sep 2026 13:27:02 +0800 Subject: [PATCH 12/12] Fix typing --- src/parcels/_chunk_cached_array/core.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parcels/_chunk_cached_array/core.py b/src/parcels/_chunk_cached_array/core.py index 53775a082..9d545c1ad 100644 --- a/src/parcels/_chunk_cached_array/core.py +++ b/src/parcels/_chunk_cached_array/core.py @@ -44,7 +44,7 @@ def wrap_dataset(ds: xr.Dataset, max_cache_bytes: int) -> xr.Dataset: for name in ds.data_vars: var = ds[name].variable if is_duck_array(var._data) and is_dask_collection(var._data): - var._data = ChunkCachedArray(var._data, max_cache_bytes) + var._data = ChunkCachedArray(var._data, max_cache_bytes) # type: ignore[assignment, arg-type] return ds @@ -112,7 +112,7 @@ def _raw_vindex(self, *indices: np.ndarray) -> np.ndarray: # Sort points by flat chunk key to group them. sort_order = np.argsort(flat_keys, kind="mergesort") - sorted_flat_keys = flat_keys[sort_order] + sorted_flat_keys = flat_keys[sort_order] # type: ignore[index] # Find group boundaries. boundaries = np.concatenate(([0], np.flatnonzero(np.diff(sorted_flat_keys)) + 1, [n_points]))