From 52adba9c12aed607b417f917744e6b95917c5c05 Mon Sep 17 00:00:00 2001 From: unalmis Date: Mon, 16 Sep 2024 01:19:46 -0400 Subject: [PATCH] Improve quadrature over velocity coordiante for effective ripple --- desc/backend.py | 56 +++++++++++------------ desc/compute/_neoclassical.py | 46 +++++++++---------- desc/compute/utils.py | 24 ++++++++-- desc/integrals/bounce_utils.py | 46 ++++++++++++++++--- desc/integrals/quad_utils.py | 29 +++++++++++- desc/objectives/_neoclassical.py | 22 ++++----- tests/baseline/test_effective_ripple.png | Bin 13812 -> 13483 bytes tests/test_objective_funs.py | 17 +++++-- tests/test_quad_utils.py | 20 ++++++++ 9 files changed, 181 insertions(+), 79 deletions(-) diff --git a/desc/backend.py b/desc/backend.py index 5aab199a2a..ecbc915cbb 100644 --- a/desc/backend.py +++ b/desc/backend.py @@ -75,14 +75,7 @@ from jax.numpy import bincount, flatnonzero, repeat, take from jax.numpy.fft import irfft, rfft, rfft2 from jax.scipy.fft import dct, idct - from jax.scipy.linalg import ( - block_diag, - cho_factor, - cho_solve, - eigh_tridiagonal, - qr, - solve_triangular, - ) + from jax.scipy.linalg import block_diag, cho_factor, cho_solve, qr, solve_triangular from jax.scipy.special import gammaln, logsumexp from jax.tree_util import ( register_pytree_node, @@ -98,6 +91,31 @@ jnp.trapezoid if hasattr(jnp, "trapezoid") else jax.scipy.integrate.trapezoid ) + def execute_on_cpu(func): + """Decorator to set default device to CPU for a function. + + Parameters + ---------- + func : callable + Function to decorate + + Returns + ------- + wrapper : callable + Decorated function that will always run on CPU even if + there are available GPUs. + """ + + @functools.wraps(func) + def wrapper(*args, **kwargs): + with jax.default_device(jax.devices("cpu")[0]): + return func(*args, **kwargs) + + return wrapper + + # JAX implementation is not differentiable on gpu. + eigh_tridiagonal = execute_on_cpu(jax.scipy.linalg.eigh_tridiagonal) + def put(arr, inds, vals): """Functional interface for array "fancy indexing". @@ -123,28 +141,6 @@ def put(arr, inds, vals): return arr return jnp.asarray(arr).at[inds].set(vals) - def execute_on_cpu(func): - """Decorator to set default device to CPU for a function. - - Parameters - ---------- - func : callable - Function to decorate - - Returns - ------- - wrapper : callable - Decorated function that will run always on CPU even if - there are available GPUs. - """ - - @functools.wraps(func) - def wrapper(*args, **kwargs): - with jax.default_device(jax.devices("cpu")[0]): - return func(*args, **kwargs) - - return wrapper - def sign(x): """Sign function, but returns 1 for x==0. diff --git a/desc/compute/_neoclassical.py b/desc/compute/_neoclassical.py index 3eb1c143ee..96a0c054f0 100644 --- a/desc/compute/_neoclassical.py +++ b/desc/compute/_neoclassical.py @@ -19,7 +19,7 @@ from ..integrals.quad_utils import get_quadrature, leggauss_lob from ..utils import map2, safediv from .data_index import register_compute_fun -from .utils import _get_pitch_inv, _poloidal_mean +from .utils import _get_pitch_inv_chebgauss, _poloidal_mean @register_compute_fun( @@ -79,10 +79,10 @@ def _G_ra_fsa(data, transforms, profiles, **kwargs): @register_compute_fun( name="effective ripple", # this is ε¹ᐧ⁵ label=( - # ε¹ᐧ⁵ = π/(8√2) (R₀/〈|∇ψ|〉)² ∫dλ λ⁻²B₀⁻¹ 〈 ∑ⱼ Hⱼ²/Iⱼ 〉 + # ε¹ᐧ⁵ = π/(8√2) (R₀/〈|∇ψ|〉)² B₀⁻¹ ∫dλ λ⁻² 〈 ∑ⱼ Hⱼ²/Iⱼ 〉 "\\epsilon^{3/2} = \\frac{\\pi}{8 \\sqrt{2}} " "(R_0 / \\langle \\vert\\nabla \\psi\\vert \\rangle)^2 " - "\\int d\\lambda \\lambda^{-2} B_0^{-1} " + "B_0^{-1} \\int d\\lambda \\lambda^{-2} " "\\langle \\sum_j H_j^2 / I_j \\rangle" ), units="~", @@ -106,12 +106,7 @@ def _G_ra_fsa(data, transforms, profiles, **kwargs): resolution_requirement="z", source_grid_requirement={"coordinates": "raz", "is_meshgrid": True}, quad="jnp.ndarray : Optional, quadrature points and weights for bounce integrals.", - num_pitch=( - "int : Resolution for quadrature over velocity coordinate, preferably odd. " - "Default is 75. Profile will look smoother at high values." - # If computed on many flux surfaces and small oscillations are seen - # between neighboring surfaces, increasing this will smooth the profile. - ), + num_pitch="int : Resolution for quadrature over velocity coordinate. Default 50.", num_well=( "int : Maximum number of wells to detect for each pitch and field line. " "Default is to detect all wells, but due to limitations in JAX this option " @@ -152,7 +147,7 @@ def _effective_ripple(params, transforms, profiles, data, **kwargs): if "quad" in kwargs else get_quadrature(leggauss_lob(32), Bounce1D._default_automorphism) ) - num_pitch = kwargs.get("num_pitch", 75) + num_pitch = kwargs.get("num_pitch", 50) num_well = kwargs.get("num_well", None) batch = kwargs.get("batch", True) grid = transforms["grid"].source_grid @@ -171,24 +166,29 @@ def dI(B, pitch): return jnp.sqrt(jnp.abs(1 - pitch * B)) / B def compute(data): - """(∂ψ/∂ρ)⁻² B₀⁻² ∫ dλ ∑ⱼ Hⱼ²/Iⱼ.""" + """Return (∂ψ/∂ρ)⁻² B₀⁻² ∫ dλ λ⁻² ∑ⱼ Hⱼ²/Iⱼ. + + Notes + ----- + B₀ has units of λ⁻¹. + Nemov's ∑ⱼ Hⱼ²/Iⱼ = (∂ψ/∂ρ)² (λB₀)³ ``(H**2 / I).sum(axis=-1)``. + (λB₀)³ d(λB₀)⁻¹ = B₀² λ³ d(λ⁻¹) = -B₀² λ dλ. + """ bounce = Bounce1D(grid, data, quad, automorphism=None, is_reshaped=True) - # Interpolate |∇ρ| κ_g since it is smoother than κ_g alone. H = bounce.integrate( dH, data["pitch_inv"], + # Interpolate |∇ρ| κ_g since it is smoother than κ_g alone. data["|grad(rho)|*kappa_g"], num_well=num_well, batch=batch, ) I = bounce.integrate(dI, data["pitch_inv"], num_well=num_well, batch=batch) - # Note B₀ has units of λ⁻¹. - # Nemov's ∑ⱼ Hⱼ²/Iⱼ = (∂ψ/∂ρ)² (λB₀)³ ``(H**2 / I).sum(axis=-1)``. - # (λB₀)³ db = λ³B₀² d(λ⁻¹) = λB₀² (-dλ). - y = data["pitch_inv"] ** (-3) * safediv(H**2, I).sum(axis=-1) - return simpson(y=y, x=data["pitch_inv"]) - # TODO: Try Gauss-Chebyshev quadrature after automorphism arcsin to - # make nodes more evenly spaced. + return ( + safediv(H**2, I).sum(axis=-1) + * data["pitch_inv"] ** (-3) + * data["pitch_inv weight"] + ).sum(axis=-1) _data = { # noqa: unused dependency name: Bounce1D.reshape_data(grid, data[name]) @@ -197,13 +197,13 @@ def compute(data): _data["|grad(rho)|*kappa_g"] = Bounce1D.reshape_data( grid, data["|grad(rho)|"] * data["kappa_g"] ) - _data["pitch_inv"] = _get_pitch_inv(grid, data, num_pitch) - out = _poloidal_mean(grid, map2(compute, _data)) + _data = _get_pitch_inv_chebgauss(grid, data, num_pitch, _data) + B0 = data["max_tz |B|"] data["effective ripple"] = ( jnp.pi / (8 * 2**0.5) - * (data["max_tz |B|"] * data["R0"] / data["<|grad(rho)|>"]) ** 2 - * grid.expand(out) + * (B0 * data["R0"] / data["<|grad(rho)|>"]) ** 2 + * grid.expand(_poloidal_mean(grid, map2(compute, _data))) / data[""] ) return data diff --git a/desc/compute/utils.py b/desc/compute/utils.py index 26ed64c081..e7d496a045 100644 --- a/desc/compute/utils.py +++ b/desc/compute/utils.py @@ -8,7 +8,7 @@ from desc.backend import execute_on_cpu, jnp from desc.grid import Grid -from ..integrals.bounce_utils import get_pitch_inv +from ..integrals.bounce_utils import get_pitch_inv, get_pitch_inv_chebgauss from ..utils import errorif from .data_index import allowed_kwargs, data_index @@ -728,12 +728,28 @@ def _poloidal_mean(grid, f): return f.T.dot(dp) / jnp.sum(dp) -def _get_pitch_inv(grid, data, num_pitch): - return jnp.broadcast_to( +def _get_pitch_inv(grid, data, num_pitch, _data): + _data["pitch_inv"] = jnp.broadcast_to( get_pitch_inv( grid.compress(data["min_tz |B|"]), grid.compress(data["max_tz |B|"]), num_pitch, )[jnp.newaxis], - (grid.num_alpha, grid.num_rho, num_pitch + 2), + (grid.num_alpha, grid.num_rho, num_pitch), ) + return _data + + +def _get_pitch_inv_chebgauss(grid, data, num_pitch, _data): + p, w = get_pitch_inv_chebgauss( + grid.compress(data["min_tz |B|"]), + grid.compress(data["max_tz |B|"]), + num_pitch, + ) + _data["pitch_inv"] = jnp.broadcast_to( + p[jnp.newaxis], (grid.num_alpha, grid.num_rho, num_pitch) + ) + _data["pitch_inv weight"] = jnp.broadcast_to( + w[jnp.newaxis], (grid.num_alpha, grid.num_rho, num_pitch) + ) + return _data diff --git a/desc/integrals/bounce_utils.py b/desc/integrals/bounce_utils.py index 4d4263e849..9b784d334c 100644 --- a/desc/integrals/bounce_utils.py +++ b/desc/integrals/bounce_utils.py @@ -14,6 +14,7 @@ ) from desc.integrals.quad_utils import ( bijection_from_disc, + chebgauss_uniform, composite_linspace, grad_bijection_from_disc, ) @@ -37,7 +38,7 @@ def get_pitch_inv(min_B, max_B, num, relative_shift=1e-6): max_B : jnp.ndarray Maximum |B| value. num : int - Number of values, not including endpoints. + Number of values. relative_shift : float Relative amount to shift maxima down and minima up to avoid floating point errors in downstream routines. @@ -45,20 +46,53 @@ def get_pitch_inv(min_B, max_B, num, relative_shift=1e-6): Returns ------- pitch_inv : jnp.ndarray - Shape (*min_B.shape, num + 2). + Shape (*min_B.shape, num). 1/λ values. """ # Floating point error impedes consistent detection of bounce points riding # extrema. Shift values slightly to resolve this issue. - min_B = (1 + relative_shift) * min_B - max_B = (1 - relative_shift) * max_B + min_B = (1.0 + relative_shift) * min_B + max_B = (1.0 - relative_shift) * max_B # Samples should be uniformly spaced in |B| and not λ (GitHub issue #1228). - pitch_inv = jnp.moveaxis(composite_linspace(jnp.stack([min_B, max_B]), num), 0, -1) - assert pitch_inv.shape == (*min_B.shape, num + 2) + pitch_inv = jnp.moveaxis( + composite_linspace(jnp.stack([min_B, max_B]), num - 2), 0, -1 + ) + assert pitch_inv.shape == (*min_B.shape, num) return pitch_inv +def get_pitch_inv_chebgauss(min_B, max_B, num, relative_shift=1e-6): + """Return Chebyshev quadrature with 1/λ uniform in ``min_B`` and ``max_B``. + + Parameters + ---------- + min_B : jnp.ndarray + Minimum |B| value. + max_B : jnp.ndarray + Maximum |B| value. + num : int + Number of values. + relative_shift : float + Relative amount to shift maxima down and minima up to avoid floating point + errors in downstream routines. + + Returns + ------- + pitch_inv, weight : (jnp.ndarray, jnp.ndarray) + Shape (*min_B.shape, num). + 1/λ values and weights. + + """ + min_B = (1.0 + relative_shift) * min_B + max_B = (1.0 - relative_shift) * max_B + # Samples should be uniformly spaced in |B| (GitHub issue #1228). + x, w = chebgauss_uniform(num) + pitch_inv = bijection_from_disc(x, min_B[..., jnp.newaxis], max_B[..., jnp.newaxis]) + w = w * grad_bijection_from_disc(min_B, max_B)[..., jnp.newaxis] + return pitch_inv, w + + def _check_spline_shape(knots, g, dg_dz, pitch_inv=None): """Ensure inputs have compatible shape. diff --git a/desc/integrals/quad_utils.py b/desc/integrals/quad_utils.py index 5f7c5994bf..2f85f64058 100644 --- a/desc/integrals/quad_utils.py +++ b/desc/integrals/quad_utils.py @@ -1,8 +1,9 @@ """Utilities for quadratures.""" +from orthax.chebyshev import chebgauss from orthax.legendre import legder, legval -from desc.backend import eigh_tridiagonal, execute_on_cpu, jnp, put +from desc.backend import eigh_tridiagonal, jnp, put from desc.utils import errorif @@ -139,7 +140,6 @@ def tanh_sinh(deg, m=10): return x, w -@execute_on_cpu # JAX implementation of eigh_tridiagonal is not differentiable on gpu. def leggauss_lob(deg, interior_only=False): """Lobatto-Gauss-Legendre quadrature. @@ -191,6 +191,31 @@ def leggauss_lob(deg, interior_only=False): return x, w +def chebgauss_uniform(deg): + """Gauss-Chebyshev quadrature with uniformly spaced nodes. + + Returns quadrature points xₖ and weights wₖ for the approximate evaluation of the + integral ∫₋₁¹ f(x) dx ≈ ∑ₖ wₖ f(xₖ). + + Parameters + ---------- + deg : int + Number of quadrature points. + + Returns + ------- + x, w : (jnp.ndarray, jnp.ndarray) + Shape (deg, ). + Quadrature points and weights. + """ + # Define x = 2/π arcsin y and g : y ↦ f(x(y)). + # ∫₋₁¹ f(x) dx = 2/π ∫₋₁¹ (1−y²)⁻⁰ᐧ⁵ g(y) dy + # ∑ₖ wₖ f(x(yₖ)) = 2/π ∑ₖ ωₖ g(yₖ) + # Given roots yₖ of Chebyshev polynomial, x(yₖ) is uniform in (-1, 1). + y, w = chebgauss(deg) + return automorphism_arcsin(y), 2 * w / jnp.pi + + def get_quadrature(quad, automorphism): """Apply automorphism to given quadrature. diff --git a/desc/objectives/_neoclassical.py b/desc/objectives/_neoclassical.py index 2936fe6ebb..302db58cfa 100644 --- a/desc/objectives/_neoclassical.py +++ b/desc/objectives/_neoclassical.py @@ -42,7 +42,7 @@ class EffectiveRipple(_Objective): locations. Defaults to 0. bounds : tuple of {float, ndarray, callable}, optional Lower and upper bounds on the objective. Overrides target. - Both bounds must be broadcastable to Objective.dim_f + Both bounds must be broadcastable to Objective.dim_f. If a callable, each should take a single argument ``rho`` and return the desired bound (lower or upper) of the profile at those locations. weight : {float, ndarray}, optional @@ -68,18 +68,17 @@ class EffectiveRipple(_Objective): Should have poloidal and toroidal resolution. alpha : ndarray Unique coordinate values for field line poloidal angle label alpha. + knots_per_transit : int + Number of points per toroidal transit at which to sample data along field + line. Default is 100. num_transit : int Number of toroidal transits to follow field line. For axisymmetric devices, one poloidal transit is sufficient. Otherwise, more transits will give more accurate result, with diminishing returns. - knots_per_transit : int - Number of points per toroidal transit at which to sample data along field - line. Default is 100. num_quad : int Resolution for quadrature of bounce integrals. Default is 32. num_pitch : int - Resolution for quadrature over velocity coordinate, preferably odd. - Default is 75. Profile will look smoother at high values. + Resolution for quadrature over velocity coordinate. Default 50. batch : bool Whether to vectorize part of the computation. Default is true. num_well : int @@ -99,7 +98,7 @@ class EffectiveRipple(_Objective): def __init__( self, eq, - target=0.0, + target=None, bounds=None, weight=1, normalize=True, @@ -108,16 +107,17 @@ def __init__( deriv_mode="auto", grid=None, alpha=np.array([0]), - num_transit=10, + *, knots_per_transit=100, + num_transit=10, num_quad=32, - num_pitch=75, + num_pitch=50, batch=True, num_well=None, name="Effective ripple", ): - if bounds is not None: - target = None + if target is None and bounds is None: + target = 0.0 self._keys_1dr = [ "iota", diff --git a/tests/baseline/test_effective_ripple.png b/tests/baseline/test_effective_ripple.png index 38bbd3846b2e518fa83d9fec4a4530892612c89d..0a19da2c30efa05e5a51d6cfe429b479cdb4c1b6 100644 GIT binary patch literal 13483 zcmeHucUV(dv~K{hpML{|aU7FHsf@O3>iX`-|sGyV}JrIgI ziUBEt(u+zk)F7dR010pJqt4uS-~HZu_y6nnDV}q7S$+N1+M9dlO!T&I-M1BmLTxwD z|NR0A#U+kHabDiM5xhwYym<}$Q1#cj>~H4d;vaOy*BNDe#s8Y8kH4q;AL4<|zJBgL z-b(VS^2d&dyZQTH^HWn$@cL(fypOM|LYE+IKUieTHGNAz6pHT({L7KAmFJE^9skYX z_fr>xGiHXuQfvbvX6M?oB(!(E6@Pgx`S%T59^s-dnQ3&JZ|F}8ewn2fx%8;QO4a>V zeg$Q7-sy_Xzk5Dg5%#&{A9^wW=GG6*a~B@>o_K96d@3R1rueRyE1QpPOO0S(j^6(u z%%!_Oxo0Y5N$CpT?s{^pOM`3C$ZU_Ba(XME>taqGjRS>>r0v%R|2A*lfkNdda&ADO z@}p7U{o#ZE^ZLIh6QiOql(J=AXn=XKO)Alk>awyh5fI0PY8)-JZ_q3-FKWIgot45` z8X(ZF_v29R`(-VP-@E5!nlpYwp*)+OB-^E#nW8Ib$B>VH)HdS_>>F8HSWMmPbN8HE zlAPZ}?G1=Kr9&DXdNVa-3w}5T2cOK#4D2V2RK1wrpwA=NFmsi?#V-sT0m*Np@gkhoJEh0=E3u=!(tLqh{?Xa0jq_gfVhv{HXIzNNUxJqg*>7$ad_PENzuF+xKi&OSprGf(NVe^t=qH%5K^+Y)kP^c zS824l>2F9_$=(>(EBEXEFZp#%e5u-rLRG%xT(~C*KggOBOw5aRiSIUmRGDw%O#)Q?l#EXN`?N%1Oz;ripBJ)W zeXO3H`yn5!iO5Ozh<0``bmhv=d0NJ_dh1A|NK_CpPs>MTInhr%dwPjEqoSDpeJP*|oEX>m( z9fy2!Q_hUkOoal`cxxLZ$rw@4mb#C8DM3Z}Lj4CvKhoz_%MU*~hm5+jD26wLo8$u8 z!;#GcsF~!n6P!+p$nXU2aXg}*+^7Q~lV&Yu9jU^|FlSa11uEi!ljjjZ;XUGbW3Y@i zPMxMN8mt{@qzt(`M8HV!=a~NZe-Bn;#+_Wh+Ko1FMJ6h0+ptpLo{_~>*Bup)JF)ck zLZWPEmj>{g#P*0^O0i;~G(MC%w5Jp~0vGD&v#0~AU;lE-j!;HdJ@G;DZQuT+{dD}$ zSNA)W@NkXv1TR9p3rA=LLip5|yZChPUE-H^0m1keS^F3(d+zx4^hiH66z3RM z6ox}W|F{(};Q3Tw=D=dn$aIfu9Wxe}pltk~Y!||%`|ykZe6sd7(AC4&*MlwPku5Lo zA%A>EF3Se)Go%9q-;oIT{FdUhvhdFk%7%qdgn_n|MKYsaT6+>2{6%&#Y% z_~{pV#ErMFLT%TWoKL8=pCP&eUd147Bg+$RiOc;R})e0 z+{ly4$;)sNMEoQ`O6yr(EfUQD`mH}{Zc!{ocL7W@R5wx z)(71%mW*Epks=i@?5jV7KWYM3e*soUBddS(RO6`Ii!9Lj)5)!3pu7&oC}>x4`*=J57oc<)#3Pb@AEQq4C(Il0j>l4s6zW6KXCa$|mo!&P{2>lYvzOPJ#xr@8qU zBg|dcLneigr~3j-UnL=;i{jdBoesHk^c-~K2rWSL6FuN#DjI$+vRn~i1(?MmA`jj3 z5_-b%(`BYZ_>vRIdRI6?r;sVP^+XN4t4~Fit05%(CTh}65D+gyWVRva!mR2Z7KQ+`_7c(>!~{wfk>C5@)va zN}KmB&+ktFpp78+K01>)zG>|9&qW4~OLL~btU;~!5T z@|fjGwPOMSHwU4Xyo!Kh#Tq;{C4}t?ybX6OQstOzzWR=cad{VDpb3PvW@f4#X<>>h9 z>B4!5u$e^R2!BA15<-{{0?xz_LUPf6ln=vF1 z?4{iTGMFDSvH4(}gvy|}wi=?)C)mNZ%7+ z0DTBuj{e7q6H*+VWaMeI+Zqpi7^dbcr3i;02K|ZwcqDcxg4Db+?)2^bKSS;SsXWAbLu=D zwRa(p=6sqd+xU1Mc*6{o|LLdmTO%u;IGOS*$)840TE*cgifX$7O!B*)!&sP~$UTl= zk#p&rMD@LEjU)YcLj&+ayq;TM1uFdGQyH9~5-vCPIhb_-&N2qTo*%?HK+Yv(#9Ph< z8jPL!ZxBeo`GNP3r;c;X^O61na7nqS1ASq$ly?WGvTh%S%DsC2XO4HbnDgaAo#O#4 zA!g-UU?E*#zU6^8w50sO&kk@61lo*69Y{SPtGf2RM*87ZD3|6K!n)gWkOH2oF4S}X zSKin7_da1G6V)fQLzQ0tF$<c-c6ZbwK}(S_h#w3i zpe9@pB0((xI6%;lGvGi;T9xtuaDbq~=F5?>P>;J>pK^~wtsVC5u_dz@K%>%d&A*p9 z7wYRnsX2CVsk3`@d^=8!-J(Gc$iKZPCIR1Wy&Tzg63Dpmv312^%QD?aM75FJE6-;G7k$&6d{e4NMSOV$$PR!!NdbfT zmE$`OBkR89^cRKhbo7vDhy?M=MvhnCkoSKCpR{bSi!4z^_>qZv{6U}pNPgHZE0`i3 zE^OnK(i2HsSOUT?LiRoM6PzJv%HdT^tdz>I8*7s(w4Ib|KQn#s0l3<2*{H{P(yB|n zK=z*yvTOiZVwvwVKvFH=hYa!KRI`Srn{rM_Swdyxh;}bTD~-2net7v0d?R)7-z49G zK(XCGOuv%66|DFiPFqH{RU(_b|C!`_VZ3O(O`Le@G^)$-sv8=Ulo80!)8~shZpiHv z2dYYWBooxo@7E3}yo7v@<{tNiWun9BLhEGU=Q}E4$#!(VH|dz5km1S-Fnx1qdnr3~ z&KI+SzbQdI{d^HG9OkL$Wi^EFXojNQ{fg5c0U!=l=Z(WDW|YHJccB5TI(MGa>HpzwB8s2_m@Y^4b*e<(y>q9Nn$P%A;sqK4q0D9qhd2S! zc_Hbwpg5_of83;Q7#OJ~1xV)S>F5)LhGtv&qaOItatZ81f<;L#GRgDW(SM=^0?q>a zfHUdl(hS@Lkhv}#Y6!?FF{ITh5U1ohUbR6jOz#4#??pxcWc$n4XMT?aO%B7iJy_sv zZP4-m9hVu{jwaVJw(bbbg1XAv3mv|s3{s&E;Ge1=?hNmPsi%YWah0UXve)E)d~BbV@75|x9>gFjV)0})-LQb z+PG;5^QIWlW67lgoc-FXy*qv^c^B+ML-r|w0pNHab{%NY7SUh|&>*ZR8hHhUWCe0~ zb_Bv?;&69uT#wve!FhWKp<+LOj?zwKqsRY7fk0?oz5jd~sQ4KLBA~!#YFJB@mz%gZ z*FxVeu?1{^;2Ol+pADC<{Z-cks||$Ftyi9>BUJ`Y^a`cKQ|<_D*u47;oM^1Ie5w>v zhr03ybfHRsS#=LumgLY3#y23NK^_xq*AEQ?)@}l$Qm+(PNFRPG+9zV{O*$WW^g1Nm zI$)}Zd}ll7?RyAA&;P0Ccd7jeDJo#kYo-56^bgJ8X39#SHTT4@y^R#oL>(2w>X;X) zX~G6xe0+S3UM?4GT>=o*^zNNLznq%t>({PMp3@l;G10Ngeq>2mzt$7d)S4-ILjo(f zZ+kId=DUY$X28qbx$75jZ*$v~6_o=o+_{_V%%t;73CHn5>X$2i_V>K!{6Lp)EK0Vx zDENRHv=(uLPETbeH8rL5;&2VSgP`XJo(qKU^Fn^~amV+?7^PRrle;-{Q{B6|pM2=G z%4tOW7~?hfb+ogq%=dow)C5t@Fo{(1_1^e#C?fl?)q(c+35NtEiHFmRZy)dK^1Jk{ z-9Fy0J;EcV>jQ#w8(O+Bu(!QH;!Un;Z@5~V2%U4$xbh25`I#Zf2; z15MTh{rbd{I@QIGez*jlt|F$x0U{^J_cW3cz_WNWQ#rBGzc=lMVauPWCweH?acP=@ zGn~x2V?yOGG#?s{S=!Xvk)|KD~GJ&U4^Mk@AXaR>iBq&1N_trkcEmLPTyJ zRL*_U>JS;rcBRj*XK73r1el&t@#Wrt$}tRI9UzBL6Rm0Wp8;$=Qw6=#_H%NBk%_t- zvvs6=OdjcnhQz){NfZ@KE?Ad}sMvVzs{k+rYcOG{3J7l+!iSEHK=8IzZ9q~g2m?3buo}!YYKvEp z7tv7N58A*R*7v*Q&L`w-wmW^&0?vdLSH~++%teD! zdis2mks{6N<2tK>%%n=shM3X^0SMi0v<=xSCwrX)4B?G}$?#GokUtXftc=Q1Npmse zGklWSYfhk_CM5DV6qnq$Trp$$1`XnLDbkeQD3z485#_cY)c$K@_EJ8>0<8sv=4}z^1MndDjyMeS<)mly`#+D-f&36Fm~MVfOT4h7|1U zbK#rdekX)Rd>K_(b=7%fq@=c6NUKX3Z70g1^LRQrn`KEaum>E{KF6JubN_i@5;7m# z{_t-7f%$JD*)Sf1fIF1S(Rh31ogNO$>Yxx!5uU*UnzW-r`}i;I}<*CT{$eCX{pVE!#F z-^{TIv;-33@-&>g;#}K|_Ql9y4x!2u##->89n2op@+2At?G!i42j%A~P4>*!;jk*< z-$vu?Sy$%?nVM+T&yg;=daf>AiE&LHl!td+C=s(BF0Ex;@zyX!Mv6i;GuAY@V8I zy)AZq-chWs&LH>0J5tD0>D5b0*0B#C-RhCGKnHdwjQwpbf7_R z!)~cmW$Y7yJ(LkNoQkjM4IQcxsl5OTtb=iYEC3os`9PD53C@&Mu$TrG?vPoack!G1 zaV?U`h2j3&-+c=WEpEd}Dj4BFf}3b4JE-g_GfH`JgT5?=oefLRvaX`t(`qWTwB}1_ zS`S{DcADJ4R`3L-1Cs5jViLVSI=1b7x%*UU@Z5Dgo4O(LUY94dYa;t~Qmuc1Zr6Su zIY(IGMMM}xEi4RXGv^ihCq%MELmr3Elj>ymYO&lw|2!Q#BO2z(Z2sF(Y?Elpbi+Dj zxTf4)pwKza@~>ed)AmoPT38+DS_kw(N0Otn2S?eJ*Wc_%+$2_MwJg~qqgum3BSq^v z<;*q*X_u~9a_V*=k;vs%&~J>IKl-z-zj3GxlK^vkO=RfojCwHorv_50drFv&$PB5SLvRC;WK zERz%(j#}7UKMzdEslWVgsT27+h^dpqZtW@Y-~fx+o8?lrt;UPHnZ-!1x*em%SS~mV zEX1SYi8k?D8V{W{n;q{m|H9%=<#j*bxi?fJih@|=7G;f$Bi2c4%$0Zr%b;4LwYkX)Do!rlcWv!z%68RI}mdAKTQf*PSU#U85%t#3r_8t80>b zN-F}6`BAeTsL!rhRWxiy1@_+APUf*Y=iR7y@-iE366VY8#2|8ZcJ52njL>aUC!})m zd3y5y2BWd(r-|}q{bpHXQ;aY3B>kP&dMz-&yTXb}jKQ+S%$KyfAYIMxH(ZjAxI57X z?`82v4h3{ul8dUOdh8W6wIj>RU6uxXY?I_?NGhMFt95hZ(My+;ou<@po!B?!`z^U~n!%;4eqBiA>Cy#Yu`^wqY!~(uPv)c0T0IwCQ@e&2;q>+d5Q=4QOXz7xc)6_EyE~s15q} zuRhJAK7XIxaDbZSaqAlH*STy~R>1|kEe4|NREmM=Sv%IGuXWQZ&9I&|LszJqDmJuZ zHv3wC9$BrUoDn-qtFLXd->w}5Vo^^+J&60ev#;swduw@n`1D(W6Ej4rrs{Z(s_JC^ zYo)~xwL3E{j$b4BGG1%tdGem212pWvH32|mFcUR$JgjxnrZ+$*lRd?7& zqjt}LeTLV*QSRFtB^`MNt@J;_FlKjEz4Oi$?;YZDe;03^&L8R%^lVL|Ocu2KMx1=^ zt5ST6_biA4in?}m92rgJ&JhCLyNHOf#Bn}M3C%r4f0!_%IWm7+xm;cS7Me*mz+&mC z?FrUvc7De6yV#U))|7MP`z%|U>3V!zf@=$^#1Djkb!;GQ(05K1_6KY|`^zaHsSw8> z(I5Ohi`7!4W>5qgWcZrZzBf{s?WS#L6;N$Ui~4S9+|6noNwCn?oENmF!6IJVAVitG z#v~7n*vUII%+^z~*Hyt$DSp0PA*d%!T}73(;562(@#UHk3P-?HtUYZD8#ZIFOMU;H zb3^YV){>zeqsG@`WOY^{*O82)7}}pTQG2$wMwU*F2W4DU`AV2G>ld@EUfBz27=9ax7<4Tf1XN{S;cPpfW?Hqz6O6&>4+ zKGsNC{k*9~_HQwk8xg|``wr7dEXIg`{|X4z$qhDkasZUzN8=P;Q%(>dg7>|yuAN6* z?4}g?7_+}BWH30?Pe-D5gIgK&@rpmiuJuA>qJ(xsa8koxB zcVMt4Ey#!RBuvAsfVY60KLZ-)_(vaPOfPG#PAmRSt5+vdaPTZ2q(N7W+;Gy2*Vx2# zoW=#l-zfFj%vmGM_ghX5jr)sX;ePJkEk5~ner9KKvvmcGkPWOVfb`B{rmC8)i4Z7H zyMp{1#HRXKF?laePZr=A4O-0A2GcCoG=L#X`?VvRTQi6g<%?moLjq+D%vQ2oRYG3r zY?VC`L-~3iFSZ^$rPU7LM3OvO9iodfd=Wa|TGfA5d6Sr9+7>H)?KlW#|K5#%VhBKU zKR86a`l+jo={{9pE-7m+VfHkk%2fmtOf2gLS`4SEgXZm6PfT>Q(}zcRfIjY*Uj5WV z$7hEhD=?3cOSUPMxzC?@62(~oz-~8FFq_Bhsm4+r;^g3tDdkyQYN z`oVE;zX?DR3lmc~Pr0LMM8&5Q;{9JF~yj+~rb zxJ_uNIjD2dXbY-pd&m_~;`vq;^_8_)OTw%ut@SMUT+Tt8hybEeHkifO0%m0b43Vm} z{%B#rCehg9blSuO@$fvA@J^+H**zF&9irV^JCYe4oxt!E^VKKjvoidexpoab;r48g zU`WEyYard{D%3`WIMAv| z=0fyYf08VvCwyL1bZL;9bc+5|Mbe(_TZ)OOYmmxWJf1%GO`z3Vb~riq72Sr zQGu7h7V`DW88+e9Lw1Eas?Tm)8d8Dq`IuT7EC3TNlH&NV)ug=DCkn>=3i3XErFBx@ zWLd6PaAS#`vC6Q5}8SX}7BkF2wR7l6wxk;Of3u)Rqd6&bTU$QZK7Li_e!88c7k z;ZM-Qv{*O)nJ=XL5q<1vUpK)NbjV50^K0|K($#>a2WR{UaWJcayEXMIpD23?GCr@D zheqUrK81cNz{O!nHWr-zv2*jchUa4GlleedUoE*fOSZ#)CQrzOVT;^JpSpmqns|_) z)`GIu*PO^W-i-2nq_4AL&R_ zafbZPpPm!Oq_voTV3x`Vl9mZrHWd##I~N##b%4%LqkNi@G!w67Am2bDU7C7cfP;CJ z$6nv$;86dzlhNs(@8=PIHWmdU?EE@`PBA}^1%=;HpV4@Uhvz^$Fe4ZR)&`h!wK zc?0&yx@=He=;{D|eT6_rm9ZG6cFV7Pk4=3qv8~9tWD_!0At$Nk#TN@Wwx zO`|6i)zn6Nc)D)f1ErT#-;GU3;r2!$i5eX~VMv~9lFnv%(s9WWrXgqPX+09zY=8QF z`;TZj|6;NA>Orfp5060ghf`{Eb81cqFsS7JtrdW6@4R+FPS|F?=mz4yb{SuU3~`m_ zbghK63tZnIyXPfNowZy0N_oYDAMhcBS-JVNx&Y*uyeUfl@&jVo^CN0jO%tR@K$jE%l%1 zOI4X&77iVkRrDF2v#eS0w7%E*0i2eA4ocg}gQ+SCAv8ZcfH9;`lg+X%(rQ{d2d0E~q`1XA^@*?eKiXG0I?*wSBOgSx zP;+mE@XyRzPFj+x!-X-+Pv^Ha%juXJ5lbplv$7HauLpddj?3icjwQZqQ7~mM+PChrXkf2kbOP`R3Mo44aDxXU!M97CV0!lr|EGw}hd zSr#8e%C$M8h{Y&ItT?tS{nqmiRI-Ai-Y%d$=?5x+kBcplLzP#~guJ7UB_~^K#?IZR zl)T2V*%S&MFF+DcRdi;sg}!o0WrG%3|NDn`rF<}c^h82J)m2+e2=ndKTcJ!>Z7S^6 z%x5bc96e-uh_D<1l4bhv#p2AG`8LX<`?%5A0oF3=`IoJmwS|Bf=}CAC(xcQwS=Vn_ zVU?Jymy@5AeBbPCC4J_Od7o-Fky2uYVYbk>j*KNsTi9o=0rc^rNy;b8I^>BXd(C(E zrW`b|Fy0Qg$(7ko!2aG zs3%p7fh%Na#*wNG7|#UZ(EM3NB_hD*q~`0)F>+OdttJC>BmrU|1nY!c7Vikzuvz== z$yt&TY!$l(S2NviN6+32Sj?pg+9VAR`+lH_v9B;s!1Je}@4nGW6zYDeR#XhQfATPw|=k%+NR@pCdEjM|nz8_k|dw1;+oJMa7V z&M|r0k`t{izEefc*K9B$5UnmufJ`b0?%tB7oymok0j66ng0n;a2jH%#P57KDarvvQ zR=dMG3omKdtcUX)ObK6KO}|O89yzPI@R05YYB$hn^x2Aod&kn0G5eUo{)3n>;QKU( zMG!USzFgf0k4ZZoN=<6qnv)d+nSkz!HP>Xc*h{T5v(q0@xWx3^RP)vb&+8=~QxuhI z5;haf(03iEelaK-6mHyvWm$mh2qD&=qV|#v?Jn3|=lX8JiBk8BPyOFU(|5dS8mk-ykSrTt62>Kiecwe{Q%w2p}*63L*3oF~=sH zGmg5>kLvY`(Itvf{yK>3sjm#DJ-my1_|W(~{h6sGcAYUsmOz+ z7ezm4(gn9@D>kUVNm=)2{j`bE-i;7x>dY>`rZKkflDl^q14aEH#%=;J1h0T{qQDdUpGUf=-@3e}Lg`}q3`@}>v1 z#&go{aH}PuIC7pCr(UpQg40sG7ORU+1XmxT{xEwGK7MBG#ytc2LSFJQa0qXWEkkPz zB3+oT_uEqwQqxh#fj?6h{wQoo*~r=dwRX4Muzb45-LczZWz z1~mT3vuh?y{brT5pkX4DlQD~7^h;1Z=cbz8nlv|Q2A=OO&Lka8F#?6%YX8S)V~%#n?x zUi*F*6jI}Ft7F2wb=TJks(4~f6r{ca#hUx^eXuRb1#6>QAw!j4J!%gJtd%r|$!@Id6E$K5}^Q#;m;Q@g^w@)R9)h*b)Vy!G1;J$qMk|zD| zut!gEE4cJH_o} QP~eY&j>+%EryXwoFCS-EGynhq literal 13812 zcmeHuXIN8N*LEl)prSLPA|l`nB|0<(1nHZq37qmOw(MVYO5O;7J^!}$1_P$BMn3`@X}IA zLAttk9+^xY^L_F*`LJzNVf9}{5fO_kb1HSiw8mg|t+$xY|8w3px!0~;t1`_#3C6td zt*&0X`fg6_pp7i=%@p+7zi$110!HyRSy@?!*9#A8rp^aG#U__)<;?LR5ZNj5adGrV zguNb-Yc+8Rtte(`SmMVXJw84}Y}V_Aw^U+WBF21wSYp>IgwC12MYV<=baZr0OXs=^ z|A;{JUGpkU7nR&`2Pxxj3ZPuSLt}sZ^6Cq#&gahK?kH6TK>=)JFHrDnAM(6sg=Sbt%+ z@tL~TET{PcS*M{%ZMTcccjY6$szd>P^&+ov4N}O$qj*U|isBDVkzcQL=Pi2*#SV`8 zo}>vZgx=k&+zbq)e0`|C-e}8~Eteh$BeHdO8+UdQmRZO@{`f=Uz=5Hu@!sbSwp|`s zr((;pA6ukBmNj=2fvB=RBiATHIuP2vI)ZS0b?N-rhjiPn(5yGRln{tWyy(Ozq)4OY zS{XF23Ly~)=@tfL@q_2C4BE!NKIA4FtQp*05{pwdyLfr9M#+*tr>uKjCu!ac)1-2%05FI|a7d{_4&z+T1mI_8;Itv9P46 zq7&T^VnQE>P2-Djv`v}wUoT-e>G3^%04=WwmiIlYh7AXM5Qp|)|LYy*a>4L0fTaq& zu7|)?I7YdQe?wCkgoDtmoE+x)*mIDc+r#eslCG(oel&N^yebja48SA{r?DUGDH`6> z+^d40ojfB)P)NZVJ%Kl&xZ4=BejY}I7Oa_KzS?7u&^p7pi3+SC zStG}5kcgwAqA5lEk7A(RC1`0RTV|wwxPOKv2M={UBNsBR02$UnXgwbY5{17(Tj-dp zR>jbN;_Lhs5>Tq}P;+l($aK~pcxcYtf1%>!FVKrYtGKegkQ9vrU-eM73%=SfZQWB? z%F9iGSF<}Kr>}30!_f%EkY~$w2h(7~)iA;T+~|)TzN!G22dcuSEr+BMaU57C|kAffSeNZ{l&jFv|a5246YvxV4QM96KlC;1x% zTEPli>T?I?^*XkA4|D`Nt-!vfM9ijtItVb93D96|m7rffXWo|ztJ;5<*WvR~Ut|zQ zOR`tjmSsEmAP=!!VmvJ)eR|FoVbT(UTxM2s2A8>;lIHnz%HWU^@=aMFkcmk z+_QrXU-%oUzFk2c``jE~P8hrhh)uy?Zo{6YSRH=ZP2egP(o$J%qqN0bol|U?0Ccv4 zlKX!}|6n!^S#4SEzo}vm{Pdfz*%b~&f%cvUR+8<5w;?4J(jw5zCRzSRF3@HN9KTtZ zP&E{j+k0-_#JeY$!zR`z1+8E1#WTf18!RgX4ttNHpqcL1tIwZ3XgYkI)(}1e3J&3k z#7DsWn)kWDZ=K%_A})na2}8Y|gHEGN83+Iyj&h+>TGXr`|05bKK<899CC}=gE3o}~ zK8Ey<1GgYymG^u%vSK8>4rR;oDikq&bQASchoUPGzNYYU>KfK8n@XFHwF2y|!XCHKe>hgfXaL4W;cnyL z4@W7VZ}sYnU|~mUH^cLR;zexv%Pa0BzyFEw^@gvG>SaE6LK+Rgb`|g~PC;*3!~(5q z0=Vnzrt4J#5uD)h7X_|vz$f?T$~2)Ux@A*l79tq!d+my>H3OJ}k7y7Y27t>EI&kj* zE(&Z#^}ieDl>KhTqNsd%gyLpj3Al(+eg$w|&qrBo1j-$SBj~O>5A`;mRaDHG??b@?11}T-oz};>myP)9 zd285W8Zw}{!WHe3i-@sbNzN=#F8v<;_2?Aw227=%(7sW}Bg2qDuP$i@7|km{b2@%S z3I8J}Xvde%!%n#h29VvqKyq>q{8|Yww;CQ+>T6pP{Fq@+`wg6|O$P8nvyLdEDi|Z0 z#!1V|av@^?mgk`}`j~eJ?@)@T5Wo>VgWq-+xXOZ`{JyzPPtX-KtN41gGR54j#MkKp zuXX4SJyrlenh4bT401M4jL*n@V}CX?s*1}4E|*)?h(PSq z0^1HY3m&=>hAMF%#0o@qKRcg3I52L~a10x^+pL=!)pPOjC!?0$~{0Xeuv!_m_qQlgofu5QHry*#4 zog^JrC0}O?ysg7`0J5&6p;_=O|E5QnnEg7+e-a*cfoxDY3brx(^{Tc??+1(6P?!+O zgHR%{1Eeu^U}@N$?gHK?^7l^2&S0KqBO=*u4!yW15>}$1Fa(ROfE^a`3jzV_Un}0B z-#djLO(?+bf*IMqh>Wkq@gxm;Z}b@;j#l*%i%#856!M##cPKE zEZac|Ds(^}9;8?jomp~akZY@d{gN*$0CLq0@7Gr-QW3&GmP6*gPrPEpm3W=Oe-b7% zM^u2XI(}69^NlKG4$Q?G{spP_)3C&Ew}h5MA^6CcG)}K_T7~az7d-C}oKnC9$L+%& z-I0ib}U;Y9T!Q84Ld7E%QR8*Q&~Y492VU4k(Wvhz{0zoL zx_bVcxev^xj@s^~Ja>Qt{mB9dD{<*LKR3>oguw53Vf4#i2F-vKP=f9B! zd9S$J56u*;Ev=8|fVn*3;b(ROb_(O8ObuQwNI<6=io3;a-<9$LQp#%oSEWEyr{3V- zK;j5FH#bs1U2HK(w01etZd3_H)JM@WMf@}hiU;xSi+(@ru@99GT>hsIclvlJ9{xKX%5#5SJzoPBi z6(BpTTLOChFoUYq)JO(cTV4TYJ$n+;L9f4&HwxZ<^Ea%m#=)};#b8UW`*kw5Qr| z?88%_$Q1ZF$XNXBI(zNGM3Aqmzg)_f6)g{GR{2|e>(hdkqy>M2S~O%Aeae6mX5P;X zI!Pn6!^E%p@gY2y`v?kfJ}`r$o!t%WW&~1i_~c?Ep?ax-vgrS$RMT=)rEp$C4dwI2-US9`(z#3kvxv9x)uS(4!b+60^1w{txELy}R=o%FE z!V*n{Ks5*f2wy@{Qc`AU=k@J7NogG&cL^FL7f-3%yS!4W7%DkPO+%*{B%bOMB7&_f zh>5#K_7fJ5s+x~d+BDU4bX4EHyIT^u$Oua?6@q*v=X<%YVk0CscPV&trNU0{CzyWS zlk?Tg=r8`9?T=x;EHu&dKYj#Mxtm{39D$^lGtccX3!ER2)7={9brsv2*U@(-yB;?4 zm@jea*4C)Dl3A`Zs?93hke7T1`LGMZa11E9ZtL`MPNFyXV|H~~qT1rqwL^SWD(Rg4 z`MF}pZ;C@Rm!@Y%xT5ncD^-dxEe!Ts@fS@uMTPPzB3#DW47JmozX-d+RltMh!7JSq zL#{iszipNAL0#4Q>}AV}@jqUAKykNn`UO28&?08K0l{m?&@(9L2&im1MPF`taO{{>m-c=udx!G4^96REt(n>mmBF6t z8ynNoLLl8Y$Ff>Q zbh>j428?V`dc%*H;Gd5UEV*=^nNlPtBn7q_s*VDKb%vM>3SORKEG}=z3JF7UHPV1`Q3ow@bP9uX7$o1+R%!pchcUd?8)rIRnAlfeu z^SJc947_BZ0K}$Ul=>}F?1^y#;$HF|1$eARiY7wFMAv$v(!bQNGkfE6? ztBOu15F+YI!E)n}7JWjti7O-;5;@akg^gLT4LDS0$49SsQ5p^GVk}9n;GQ3(g|er- z%@KM;3{rdgd9cV{xD3#K)5MgZRd+O9<1l-g96qQ$`wkLJbGeQsPvS)dFSSt0e1R6L z+r)rmnoIG4a=dw-d?40|$xMS}gCc&rNz8ICeWsK!u7_k344-26_l3X#t}^jPCjQW(o&J==*7E~pKpEC7mX8z(B#+-4q`qg0`sNY!u2c| z8xEb@UCfu;%0x4q2ws;OAR7q5ST zh#DB{4UYxfl=-nBtyT#$I0*LJ02L;D_*%R$BaH*EC?g!J&e3WOon%N0%AW@O2xQ&{ zeV(hreoA_V-M`OpU_T<;3zVN~Rth1@ATt%gag18Twlg2Jc+q45g#>hq9Y12W!{9i( zfDFBZQ+R zwniwBmo0!e-3x*@QbRufnWirxvvjt>6tD3g9$pCH4$>u{#{NrU`JMUbYLpirxZ(gT&Q4B=F6BaH#AtNf7H=8hvr8(zMbdA zFh_V<&z;j}$T(V;QJA2x`=w9~3_U5|^!=(UZtYzmZualnl%GsBmU3mJl9f>g<=st* zF`;hUMVSotFl#^f*GwKSSFj$6xj&=im}h07V09@m2*eoQvUk}C|K1W$3YzQ5KI&%& z=6@w9T^$Xjyau01WX&6A5akyZD4 z-crlbEZr3L9cz;99UB}`!Cem3a45~I(O!Mo*xVa{t6?XS7vq%;xj>F>r1DqV&Z}~c zdDevr&L$=%!9a>p zQ5?K$e;*>hJRqgalYd?7X){>I{UqeX8yFrg6!?8^C|~IPPo{FWWL#rx$#8!RI#DDm z*1s@+JnsVu*KEr@4l-`$=mmSZ3KV=h=`>BP*h zsX7FoxLoOZ;^wUkRR3SYAl6!ai?vnJTYPlhIy|;oB>HWDxHz#=Z6^eGB z+&>o2Zc%cL!n~+qeu>^B{SbR)^<|Yg}8Ip02~3^Lm3d!)@*#_)H(H8CA$6!F&Ya7GPNuPlzr1a2<1MK<*C zEnXM*H&n{A@*|tpL`O$gGMn#)FjBH!eqMiya$fx#c&^&{W9rz|6PSq{slX4oPuB!O zCLgza$mnfaB8^B(PV!ZOWQY+k4!UE(j1*voE1$Bb8*h z__vygwV^(y<8J1$oi?xbDB|y`S`V+h08wk2xMix>Ygdn&)f3=Yjz#i!;)(0=j5^ZVZyVFo(^Mv-S0BuKCqQ#7koTxeHsKI2 zH=u#o=oGUfwsxtqO};3~lgVWv$s5&SsCUEr{=KAA$6=`I{@^DA(#J0sw5EdRYn@gT zq89HmlHLHHlTWPjJhg$^JV&$TFG~bfXOQ~Vpr%Q@{>d9zf3y5HT6ZZO?8;BmTB>8W zFwP_yzhq1{WsEqhdin0@%@IHd{c5kUyv< z=^hKa#2n~N<^5a>I8m!El~BGqsLFn4Xi4O8Ky#d+Rr9mN0ZZxVvbo!QQz5-Bl)(`q zXEhL66cYxVtRJ9WeF5<%Y9U_i#nlQVN-v4ih)cG31qd6vud&%MpNykUIwG^^bn$Fm z|Dg*+y-47FJ0Ank7gEK3bSshyGO9eY*~ZgjPuUaZW3(#A9?&b>tRUL@UR8yZU9;ZJ zID8~x#D(=YV$_?0XLi0pMg0O3?CdZMKYv`hx4t4~bDEmi(ZL+rXhE3{ua+}uZ`>1U z%uoP5bK+uLRsNQ~2td_pShbWd?-@yB^$8+Bx;ZZ|CpV&6(+M>jG)j{h^Ra4x_8dz< z=8#qmDYP^7;N-k2PRdA^{dp(#u{|&t%Aevr`l4k9yr)9HcJ!uwAdxuTX!#Ltc;$g< zRh}B?bXG2Rddv;ci{Ig|U84tb^PF2>A)S2%KQcH=RGfD-XLPmZ4Jb_;DTCHkjl^0h z$x0fl#Ac=uDrz2YCya1fm>z{<)@CZN53)JqRMGjt{Sp$Cw||OzHpJP6d_3dr1^$)> zeKwYmka)_~5)6GjU8KsY*?fFgJD zq>{r7ez$QK?X1$RGh0S~l54_zD6ZVJuz-<}sy%KlYyKbUANZW&0nXb9iJiz{zNWF4 zsP=BGNPhmZ%u_pIrdCj7%^k-zQ+K2(9ahzg0-S0z5k=gbp)%tGGeJg)IXo(BzF$Hr zEA8`VW5`GA;I+`fPJ8G6c(voLx7Ak$Rv_1gHU^8@hUOOoUcv)nV% z@$5VH$ugiL3z+y~y-{OdagJUv>q#!KodWEB9~9uK0h9lclVYe&))axA$y%i{Sbt_6 zDj^0V^vj zjAljE))f0xtsJe@PX|7{+o4EKkEBiGtupYS-yd2i8K|Us(qEGxUKd0iP9j({Lj|k) zfv3>p1YwRo!-?0|EeZ-w2xE823=Ra6nJz6LxLmF-jM8D=M#+_Wr;Ts$I#p952sNv( z7)co;n{$0FqGfe$QF$C57f`y&{Ox^O`=#9^p@_NK&uE-{b)5QIKlJ+hsN=|ljJ9WXQWck-xKZ^+0yH0^+7hP(OSFZFh zC!4^jO1H5W6OkZVNho;}9qO=Gd7ELV?K7uK`z`UH1U(;V8a!7?nauCi+a0j zMevhdU^7*1ppIm%;|r)?8+BsheSg$b-x4Ljey)ShtXy{e#P3eor4+mMKvYn0j|0Wj z@BwF@rvNJ7Dzx=CJ3;zFkWj~`s}n%-o)KtGjC?qEJ4Y9i57p3Elf0OOM2|CuvErVK zU#NbWJ9I#g;K_yDq?%7ZFanx+>qQjs9YwE;y$WT0IiI>R)|Dp@Y890tt3>ST+!OCP zP}(?_%<_tq{Ke;J@+i6I-%4e<2|DUf$qY)oz->>mjE{JWm%1T$@kBq9?(P)t-xrj9 z|JG5jox(=jd^PZ5)`tv-*{vY!Vt)~**2T00JcCFa>fS8?lgF<$si`q8k?!t*7|h?@d?Kr2B+rH9yMRB z%3kj!u4vyAIEW# zl`|#G->Mu+o1@RGVjY61$_$zwZ#7WdW9)-sC!poc$e+#sqt>R=cq;+a$c!=`=xL`L zk~#aN0>&KLbPnw_teJi_lS*ixl*ZNQ^?fFgEV;8I_)Qk7^?2@()5;xt|7f1|x9|Aw zild3sMdV*3H4oe$9u5O#KDu-8ZO^+$qy>eIgeE3^>Gc8;b*Y<4rafbw7LuU)p}b*| zZTlRJHW1;^k%YGB7d2~3o~jlv>Z&9b4;;!%~ULkGrP4FlzZBLZawZk z?$-}K0Zj%3fKituhO$gmGUrZ{f0f?V&<>b*t#%PE8W3no5%&_aUVuX4VJerW#xel; zxj4S~QhAXN`fe=rFL{gby42|IXG=>9?!;4-SG#BlL9Cc7=|C=U0W5gx`6dSr&b!ky_S&9Yv!GW-IAP1(zNid$;?KHR_wVcanyS>u~%`X^R z2vP$MQx)0kY}?RerIj+TiJLe@$<$`ZU!IThR@0C(*(W*On}!`Y!oG z!?vvX>JeJ-_3q{rhG*=yy>c>ZMofI#UuP%0Ky~ei7wA{oS|27|AHv}zWXarViXktA zT}JP3vmR!q21oN6;=@qX3maS(HGUM>c;(+yF7u}o465Qfgk&sV)ahuioLl6hfp*Xl zm8(~qYky8d6TvsHm{Ti~!#1ulYs|DXlnKZ_n2B2cPq2=wJz}V2Q+nOlF@_ft^wKz- zTcE_^wvyhKj51KLdKO`w9R`%#;F5cW_+Z5wmpT-JD*4@*9Y2*;m;52aK4x@1gPa~U z7eYzIIkU=Un)9rjte@&|7J@OfO%@WeDswwD=Xa(PIR!)+8iBvC>9jNA*KVI=$;9yw z=G?VJLQ%{yHV6=jelHFuxdK0vi;esegC{lR-3q&F!^i&S!E%~F=<8fINEZcj!#pAQCTr6lqSXqqt|hu;-X z6Cc=YE`q*fmXIf{pHovSSt*kIFG)!MDr#oe_wJBql$^ONX#X%<0#%xwMy_H^OsRu& zGEqm^Gv3rDF`4=cOhLiea=<-PSr(waPoYJmnbyawcQH&RZi@jbRumn;{34#{Fh_Vx0OV(!wemQKmPilg&C6!YE>(Z>nqi)as8Dx9v~{wOx)es^t9ZBs$VLG&NAlr zO_=ZtUeY*z{66%}d}H=6b;IuSuSc8ZPXZT%>1`U@2g(hP^2bZ!5{PLY!?rQT?904` zAay7PDHheN`QlmtnU&e}86kh>P1y^DEO!WCcMF4rHZEa~&D$Y#z{cS5@lUZ?(Hfvm zWmB%l?z00M7IJ)6p2p!xf=2GO{p_7~;FEr7H2Cm&oP>@6mI!(^6=REys>5?6UC=s{ zS*zbux!90@tA5z|&TAeQ#kHe6{V9e-J~6794!j_gGn z63JBkhz1Gd__Kmk{YoJZKYe{IPcFyy(*3QV6e$>+(B-4iIam!!OCz(mG4ti)_ikB_euk;vUQmjZDcgU3+a!;b9fmc9W_^uXBsU|&eQ zS84xqtnlexS8M?EGYsep8#i%-#oH**>S%j`jVqTka%{>?G