From d681976ff530eb5968817a4c5075149ad20e2aae Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 5 Apr 2026 19:41:01 -0700 Subject: [PATCH 1/7] Fix OOM in visibility module for dask-backed rasters _extract_transect was calling .compute() on the full dask array just to read a handful of transect cells. Now uses vindex fancy indexing so only the relevant chunks are materialized. cumulative_viewshed was allocating a full-size np.zeros count array and calling .values on each viewshed result, forcing materialization every iteration. Now accumulates lazily with da.zeros and dask array addition when the input is dask-backed. --- xrspatial/visibility.py | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/xrspatial/visibility.py b/xrspatial/visibility.py index ab4de453..65139ba0 100644 --- a/xrspatial/visibility.py +++ b/xrspatial/visibility.py @@ -61,7 +61,9 @@ def _extract_transect(raster, cells): if has_dask_array(): import dask.array as da if isinstance(data, da.Array): - data = data.compute() + # Only compute the needed cells, not the entire array + elevations = data.vindex[rows, cols].compute().astype(np.float64) + return elevations, x_coords, y_coords if has_cuda_and_cupy() and is_cupy_array(data): data = data.get() @@ -217,7 +219,16 @@ def cumulative_viewshed( if not observers: raise ValueError("observers list must not be empty") - count = np.zeros(raster.shape, dtype=np.int32) + # Detect dask backend to keep accumulation lazy + _is_dask = False + if has_dask_array(): + import dask.array as da + _is_dask = isinstance(raster.data, da.Array) + + if _is_dask: + count = da.zeros(raster.shape, dtype=np.int32, chunks=raster.data.chunks) + else: + count = np.zeros(raster.shape, dtype=np.int32) for obs in observers: ox = obs['x'] @@ -229,11 +240,17 @@ def cumulative_viewshed( vs = viewshed(raster, x=ox, y=oy, observer_elev=oe, target_elev=te, max_distance=md) - vs_np = vs.values - count += (vs_np != INVISIBLE).astype(np.int32) - - return xarray.DataArray(count, coords=raster.coords, - dims=raster.dims, attrs=raster.attrs) + vs_data = vs.data + if _is_dask and not isinstance(vs_data, da.Array): + vs_data = da.from_array(vs_data, chunks=raster.data.chunks) + count = count + (vs_data != INVISIBLE).astype(np.int32) + + result = xarray.DataArray(count, coords=raster.coords, + dims=raster.dims, attrs=raster.attrs) + if _is_dask: + chunk_dict = dict(zip(raster.dims, raster.data.chunks)) + result = result.chunk(chunk_dict) + return result def visibility_frequency( From 62a2a3ace4d644e8943f184fa25b3e5d04904128 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 5 Apr 2026 19:43:11 -0700 Subject: [PATCH 2/7] Tighten viewshed Tier B memory estimate and avoid needless copy The dask Tier B memory guard underestimated peak usage at 280 bytes/pixel. Actual peak during lexsort reaches ~360 bytes/pixel (sorted + unsorted event_list coexist) plus 8 bytes/pixel for the computed raster. Updated estimate to 368 bytes/pixel to prevent borderline OOM. Also use astype(copy=False) to skip the float64 copy when data is already float64. --- xrspatial/viewshed.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/xrspatial/viewshed.py b/xrspatial/viewshed.py index 96b095f6..0e282052 100644 --- a/xrspatial/viewshed.py +++ b/xrspatial/viewshed.py @@ -1560,7 +1560,7 @@ def _viewshed_cpu( num_events = 3 * (n_rows * n_cols - 1) event_list = np.zeros((num_events, 7), dtype=np.float64) - raster.data = raster.data.astype(np.float64) + raster.data = raster.data.astype(np.float64, copy=False) _init_event_list(event_list=event_list, raster=raster.data, vp_row=viewpoint_row, vp_col=viewpoint_col, @@ -2167,7 +2167,9 @@ def _viewshed_dask(raster, x, y, observer_elev, target_elev): cupy_backed = is_dask_cupy(raster) # --- Tier B: full grid fits in memory → compute and run exact algo --- - r2_bytes = 280 * height * width + # Peak memory: event_list sort needs 2x 168*H*W + raster 8*H*W + + # visibility_grid 8*H*W ≈ 360 bytes/pixel, plus the computed raster. + r2_bytes = 360 * height * width + 8 * height * width # working + raster avail = _available_memory_bytes() if r2_bytes < 0.5 * avail: raster_mem = raster.copy() From 3e773dd0ff1e00bccef2731a250a4036cfbfe56e Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 5 Apr 2026 21:22:58 -0700 Subject: [PATCH 3/7] Add kernel density estimation module (#1143) Implements kde() and line_density() for point-to-raster and line-to-raster density surfaces. Supports Gaussian, Epanechnikov, and quartic kernels with automatic bandwidth selection via Silverman's rule. All four backends: numpy, cupy, dask+numpy, dask+cupy. --- xrspatial/__init__.py | 2 + xrspatial/kde.py | 700 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 702 insertions(+) create mode 100644 xrspatial/kde.py diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py index 1eff69a3..643f3ac5 100644 --- a/xrspatial/__init__.py +++ b/xrspatial/__init__.py @@ -35,6 +35,8 @@ from xrspatial.interpolate import idw # noqa from xrspatial.interpolate import kriging # noqa from xrspatial.interpolate import spline # noqa +from xrspatial.kde import kde # noqa +from xrspatial.kde import line_density # noqa from xrspatial.fire import burn_severity_class # noqa from xrspatial.fire import dnbr # noqa from xrspatial.fire import fireline_intensity # noqa diff --git a/xrspatial/kde.py b/xrspatial/kde.py new file mode 100644 index 00000000..753331f2 --- /dev/null +++ b/xrspatial/kde.py @@ -0,0 +1,700 @@ +"""Kernel density estimation (KDE) for point-to-raster conversion. + +Produces a continuous density surface from point (or line) data. Each +output pixel accumulates weighted kernel contributions from all nearby +features, yielding a smooth density field. + +Supports numpy, cupy, dask+numpy, and dask+cupy backends. +""" +from __future__ import annotations + +from math import pi, sqrt +from functools import partial +from typing import Optional, Tuple, Union + +import numpy as np +import xarray as xr + +from xrspatial.utils import ngjit, has_cuda_and_cupy, has_dask_array + +try: + import cupy +except ImportError: + cupy = None + +try: + import dask + import dask.array as da +except ImportError: + da = None + +from numba import cuda + + +# --------------------------------------------------------------------------- +# Kernel constants +# --------------------------------------------------------------------------- +_KERNEL_NAMES = ('gaussian', 'epanechnikov', 'quartic') + +# Normalisation constants for 2-D product kernels. +# Gaussian : (1/(2pi)) +# Epanechnikov: (2/pi) (product of (3/4)*(1-u^2) marginals, integrated) +# Quartic : (3/pi) (product of (15/16)*(1-u^2)^2 marginals, integrated) +_NORM_GAUSSIAN = 1.0 / (2.0 * pi) +_NORM_EPANECHNIKOV = 2.0 / pi +_NORM_QUARTIC = 3.0 / pi + + +# --------------------------------------------------------------------------- +# Bandwidth selection +# --------------------------------------------------------------------------- + +def _silverman_bandwidth(x, y): + """Silverman's rule of thumb for 2-D data. + + h = n^(-1/6) * mean(sigma_x, sigma_y) + where sigma uses the robust scale estimator min(std, IQR/1.34). + """ + n = len(x) + if n < 2: + return 1.0 + + def _robust_scale(v): + s = float(np.std(v)) + q75, q25 = float(np.percentile(v, 75)), float(np.percentile(v, 25)) + iqr = (q75 - q25) / 1.34 + return min(s, iqr) if iqr > 0 else s + + sx = _robust_scale(x) + sy = _robust_scale(y) + sigma = (sx + sy) / 2.0 + if sigma == 0: + sigma = 1.0 + return sigma * (n ** (-1.0 / 6.0)) + + +# --------------------------------------------------------------------------- +# CPU kernels (numba-jitted) +# --------------------------------------------------------------------------- + +@ngjit +def _kde_cpu(xs, ys, ws, out, x0, y0, dx, dy, bw, kernel_id): + """Populate *out* with kernel density values. + + Parameters + ---------- + xs, ys : 1-D float64 arrays of point coordinates. + ws : 1-D float64 array of weights (same length as xs). + out : 2-D float64 output array (rows, cols), pre-zeroed. + x0, y0 : origin (lower-left corner of first pixel centre). + dx, dy : pixel spacing in x and y. + bw : bandwidth (same units as coordinates). + kernel_id : 0 = gaussian, 1 = epanechnikov, 2 = quartic. + """ + rows, cols = out.shape + n_pts = xs.shape[0] + inv_bw = 1.0 / bw + inv_bw2 = inv_bw * inv_bw + + # Pre-compute normalisation + if kernel_id == 0: + norm = inv_bw2 / (2.0 * pi) + elif kernel_id == 1: + norm = inv_bw2 * 2.0 / pi + else: + norm = inv_bw2 * 3.0 / pi + + # Cutoff radius in pixels for compact kernels. + # Gaussian uses 4*bw; compact kernels use exactly bw. + if kernel_id == 0: + cutoff = 4.0 * bw + else: + cutoff = bw + + for p in range(n_pts): + px = xs[p] + py = ys[p] + w = ws[p] + + # Pixel range that falls within the cutoff. + col_lo = int(max(0, (px - cutoff - x0) / dx)) + col_hi = int(min(cols - 1, (px + cutoff - x0) / dx)) + 1 + row_lo = int(max(0, (py - cutoff - y0) / dy)) + row_hi = int(min(rows - 1, (py + cutoff - y0) / dy)) + 1 + + for r in range(row_lo, row_hi): + cy = y0 + r * dy + uy = (cy - py) * inv_bw + uy2 = uy * uy + for c in range(col_lo, col_hi): + cx = x0 + c * dx + ux = (cx - px) * inv_bw + u2 = ux * ux + uy2 + + if kernel_id == 0: + # Gaussian + val = norm * w * np.exp(-0.5 * u2) + elif kernel_id == 1: + # Epanechnikov + if u2 <= 1.0: + val = norm * w * (1.0 - u2) + else: + val = 0.0 + else: + # Quartic + if u2 <= 1.0: + t = 1.0 - u2 + val = norm * w * t * t + else: + val = 0.0 + + out[r, c] += val + + +@ngjit +def _line_density_cpu(x1s, y1s, x2s, y2s, ws, out, + x0, y0, dx, dy, bw, kernel_id): + """Compute line density on *out*. + + Each line segment is sampled at sub-segment intervals (step = bw/4) + and each sample acts as a weighted point where the weight is + proportional to the step length. + """ + rows, cols = out.shape + n_segs = x1s.shape[0] + inv_bw = 1.0 / bw + inv_bw2 = inv_bw * inv_bw + + if kernel_id == 0: + norm = inv_bw2 / (2.0 * pi) + elif kernel_id == 1: + norm = inv_bw2 * 2.0 / pi + else: + norm = inv_bw2 * 3.0 / pi + + if kernel_id == 0: + cutoff = 4.0 * bw + else: + cutoff = bw + + step = bw / 4.0 + if step < 1e-12: + return + + for s in range(n_segs): + ax = x1s[s] + ay = y1s[s] + bx = x2s[s] + by = y2s[s] + seg_len = sqrt((bx - ax) ** 2 + (by - ay) ** 2) + if seg_len < 1e-12: + continue + + w = ws[s] + n_steps = max(1, int(seg_len / step)) + sub_w = w * (seg_len / n_steps) + + for i in range(n_steps): + t = (i + 0.5) / n_steps + px = ax + t * (bx - ax) + py = ay + t * (by - ay) + + col_lo = int(max(0, (px - cutoff - x0) / dx)) + col_hi = int(min(cols - 1, (px + cutoff - x0) / dx)) + 1 + row_lo = int(max(0, (py - cutoff - y0) / dy)) + row_hi = int(min(rows - 1, (py + cutoff - y0) / dy)) + 1 + + for r in range(row_lo, row_hi): + cy = y0 + r * dy + uy = (cy - py) * inv_bw + uy2 = uy * uy + for c in range(col_lo, col_hi): + cx = x0 + c * dx + ux = (cx - px) * inv_bw + u2 = ux * ux + uy2 + + if kernel_id == 0: + val = norm * sub_w * np.exp(-0.5 * u2) + elif kernel_id == 1: + if u2 <= 1.0: + val = norm * sub_w * (1.0 - u2) + else: + val = 0.0 + else: + if u2 <= 1.0: + tt = 1.0 - u2 + val = norm * sub_w * tt * tt + else: + val = 0.0 + + out[r, c] += val + + +# --------------------------------------------------------------------------- +# GPU kernels (CUDA) +# --------------------------------------------------------------------------- + +@cuda.jit +def _kde_cuda(xs, ys, ws, out, x0, y0, dx, dy, bw, kernel_id, n_pts): + """Each thread computes one output pixel.""" + r, c = cuda.grid(2) + rows = out.shape[0] + cols = out.shape[1] + if r >= rows or c >= cols: + return + + cx = x0[0] + c * dx[0] + cy = y0[0] + r * dy[0] + inv_bw = 1.0 / bw[0] + inv_bw2 = inv_bw * inv_bw + kid = kernel_id[0] + + if kid == 0: + norm = inv_bw2 / (2.0 * 3.141592653589793) + elif kid == 1: + norm = inv_bw2 * 2.0 / 3.141592653589793 + else: + norm = inv_bw2 * 3.0 / 3.141592653589793 + + total = 0.0 + for p in range(n_pts[0]): + ux = (cx - xs[p]) * inv_bw + uy = (cy - ys[p]) * inv_bw + u2 = ux * ux + uy * uy + + if kid == 0: + from math import exp + total += norm * ws[p] * exp(-0.5 * u2) + elif kid == 1: + if u2 <= 1.0: + total += norm * ws[p] * (1.0 - u2) + else: + if u2 <= 1.0: + t = 1.0 - u2 + total += norm * ws[p] * t * t + + out[r, c] = total + + +# --------------------------------------------------------------------------- +# Backend wrappers +# --------------------------------------------------------------------------- + +def _kernel_id(kernel: str) -> int: + if kernel == 'gaussian': + return 0 + elif kernel == 'epanechnikov': + return 1 + elif kernel == 'quartic': + return 2 + raise ValueError( + f"kernel must be one of {_KERNEL_NAMES}, got {kernel!r}" + ) + + +def _run_kde_numpy(xs, ys, ws, shape, x0, y0, dx, dy, bw, kernel_id): + out = np.zeros(shape, dtype=np.float64) + _kde_cpu(xs, ys, ws, out, x0, y0, dx, dy, bw, kernel_id) + return out + + +def _run_kde_cupy(xs, ys, ws, shape, x0, y0, dx, dy, bw, kernel_id): + out = cupy.zeros(shape, dtype=cupy.float64) + n_pts = cupy.array([len(xs)], dtype=cupy.int64) + x0_d = cupy.array([x0], dtype=cupy.float64) + y0_d = cupy.array([y0], dtype=cupy.float64) + dx_d = cupy.array([dx], dtype=cupy.float64) + dy_d = cupy.array([dy], dtype=cupy.float64) + bw_d = cupy.array([bw], dtype=cupy.float64) + kid_d = cupy.array([kernel_id], dtype=cupy.int64) + + xs_d = cupy.asarray(xs, dtype=cupy.float64) + ys_d = cupy.asarray(ys, dtype=cupy.float64) + ws_d = cupy.asarray(ws, dtype=cupy.float64) + + tpb = (16, 16) + bpg = ( + (shape[0] + tpb[0] - 1) // tpb[0], + (shape[1] + tpb[1] - 1) // tpb[1], + ) + _kde_cuda[bpg, tpb](xs_d, ys_d, ws_d, out, + x0_d, y0_d, dx_d, dy_d, bw_d, kid_d, n_pts) + return out + + +def _run_kde_dask_numpy(xs, ys, ws, shape, x0, y0, dx, dy, bw, kernel_id, + chunks): + """Dask-backed KDE: each chunk computes its own tile independently.""" + # Determine chunk layout + if chunks is None: + chunks = (min(256, shape[0]), min(256, shape[1])) + row_splits = _split_sizes(shape[0], chunks[0]) + col_splits = _split_sizes(shape[1], chunks[1]) + + blocks = [] + row_off = 0 + for rs in row_splits: + row_blocks = [] + col_off = 0 + for cs in col_splits: + tile_y0 = y0 + row_off * dy + tile_x0 = x0 + col_off * dx + tile_shape = (rs, cs) + block = dask.delayed(_run_kde_numpy)( + xs, ys, ws, tile_shape, + tile_x0, tile_y0, dx, dy, bw, kernel_id, + ) + row_blocks.append( + da.from_delayed(block, shape=tile_shape, dtype=np.float64) + ) + col_off += cs + blocks.append(row_blocks) + row_off += rs + + return da.block(blocks) + + +def _run_kde_dask_cupy(xs, ys, ws, shape, x0, y0, dx, dy, bw, kernel_id, + chunks): + """Dask+CuPy KDE: each chunk uses the GPU kernel.""" + if chunks is None: + chunks = (min(256, shape[0]), min(256, shape[1])) + row_splits = _split_sizes(shape[0], chunks[0]) + col_splits = _split_sizes(shape[1], chunks[1]) + + blocks = [] + row_off = 0 + for rs in row_splits: + row_blocks = [] + col_off = 0 + for cs in col_splits: + tile_y0 = y0 + row_off * dy + tile_x0 = x0 + col_off * dx + tile_shape = (rs, cs) + block = dask.delayed(_run_kde_cupy)( + xs, ys, ws, tile_shape, + tile_x0, tile_y0, dx, dy, bw, kernel_id, + ) + row_blocks.append( + da.from_delayed(block, shape=tile_shape, + dtype=np.float64, + meta=cupy.array((), dtype=cupy.float64)) + ) + col_off += cs + blocks.append(row_blocks) + row_off += rs + + return da.block(blocks) + + +def _split_sizes(total, chunk): + """Return a list of chunk sizes that sum to *total*.""" + full, rem = divmod(total, chunk) + sizes = [chunk] * full + if rem: + sizes.append(rem) + return sizes + + +# --------------------------------------------------------------------------- +# Public API -- kde +# --------------------------------------------------------------------------- + +def kde( + x: Union[np.ndarray, list], + y: Union[np.ndarray, list], + *, + weights: Optional[Union[np.ndarray, list]] = None, + bandwidth: Union[float, str] = 'silverman', + kernel: str = 'gaussian', + template: Optional[xr.DataArray] = None, + x_range: Optional[Tuple[float, float]] = None, + y_range: Optional[Tuple[float, float]] = None, + width: int = 256, + height: int = 256, + name: str = 'kde', +) -> xr.DataArray: + """Compute 2-D kernel density estimation from point data. + + Each output pixel accumulates weighted kernel contributions from all + input points, producing a smooth continuous density surface. + + Parameters + ---------- + x, y : array-like + 1-D arrays of point coordinates. + weights : array-like, optional + Per-point weights. Defaults to uniform weights of 1. + bandwidth : float or ``'silverman'`` + Kernel bandwidth in the same units as *x*/*y*. + ``'silverman'`` (default) uses Silverman's rule of thumb. + kernel : ``{'gaussian', 'epanechnikov', 'quartic'}`` + Kernel shape. + template : xr.DataArray, optional + If provided, the output matches this array's shape, extent, and + coordinates. *x_range*, *y_range*, *width*, and *height* are + ignored when *template* is given. + x_range, y_range : (min, max), optional + Spatial extent of the output grid. Defaults to the data extent + with 10 %% padding on each side. + width, height : int + Number of columns and rows in the output grid. Ignored when + *template* is provided. + name : str + Name of the output DataArray. + + Returns + ------- + xr.DataArray + 2-D density surface. + """ + # -- Validate and coerce inputs ---------------------------------------- + x_arr = np.asarray(x, dtype=np.float64).ravel() + y_arr = np.asarray(y, dtype=np.float64).ravel() + if x_arr.shape[0] != y_arr.shape[0]: + raise ValueError("x and y must have the same length") + n = x_arr.shape[0] + + if weights is not None: + w_arr = np.asarray(weights, dtype=np.float64).ravel() + if w_arr.shape[0] != n: + raise ValueError("weights must have the same length as x and y") + else: + w_arr = np.ones(n, dtype=np.float64) + + kid = _kernel_id(kernel) + + # -- Bandwidth --------------------------------------------------------- + if isinstance(bandwidth, str): + if bandwidth != 'silverman': + raise ValueError( + "bandwidth must be a positive number or 'silverman', " + f"got {bandwidth!r}" + ) + bw = _silverman_bandwidth(x_arr, y_arr) + else: + bw = float(bandwidth) + if bw <= 0: + raise ValueError(f"bandwidth must be positive, got {bw}") + + # -- Output grid ------------------------------------------------------- + if template is not None: + _validate_template(template) + y_coords = template.coords[template.dims[0]].values + x_coords = template.coords[template.dims[1]].values + rows, cols = template.shape + # Pixel spacing + dy = float(y_coords[1] - y_coords[0]) if rows > 1 else 1.0 + dx = float(x_coords[1] - x_coords[0]) if cols > 1 else 1.0 + x0 = float(x_coords[0]) + y0 = float(y_coords[0]) + use_dask = has_dask_array() and isinstance(template.data, da.Array) + use_cupy = (has_cuda_and_cupy() and cupy is not None + and _is_cupy_backed(template)) + out_chunks = template.data.chunksize if use_dask else None + else: + if x_range is None: + pad = max(bw, (float(x_arr.max()) - float(x_arr.min())) * 0.1) + x_range = (float(x_arr.min()) - pad, float(x_arr.max()) + pad) + if y_range is None: + pad = max(bw, (float(y_arr.max()) - float(y_arr.min())) * 0.1) + y_range = (float(y_arr.min()) - pad, float(y_arr.max()) + pad) + rows, cols = height, width + dx = (x_range[1] - x_range[0]) / max(cols - 1, 1) + dy = (y_range[1] - y_range[0]) / max(rows - 1, 1) + x0 = x_range[0] + y0 = y_range[0] + x_coords = np.linspace(x_range[0], x_range[1], cols) + y_coords = np.linspace(y_range[0], y_range[1], rows) + use_dask = False + use_cupy = False + out_chunks = None + + shape = (rows, cols) + + # -- Dispatch ----------------------------------------------------------- + if use_dask and use_cupy: + data = _run_kde_dask_cupy( + x_arr, y_arr, w_arr, shape, x0, y0, dx, dy, bw, kid, out_chunks, + ) + elif use_dask: + data = _run_kde_dask_numpy( + x_arr, y_arr, w_arr, shape, x0, y0, dx, dy, bw, kid, out_chunks, + ) + elif use_cupy: + data = _run_kde_cupy( + x_arr, y_arr, w_arr, shape, x0, y0, dx, dy, bw, kid, + ) + else: + data = _run_kde_numpy( + x_arr, y_arr, w_arr, shape, x0, y0, dx, dy, bw, kid, + ) + + # -- Build output DataArray -------------------------------------------- + if template is not None: + return xr.DataArray( + data, name=name, + coords=template.coords, dims=template.dims, + attrs=template.attrs, + ) + return xr.DataArray( + data, name=name, + dims=['y', 'x'], + coords={'y': y_coords, 'x': x_coords}, + ) + + +# --------------------------------------------------------------------------- +# Public API -- line_density +# --------------------------------------------------------------------------- + +def line_density( + x1: Union[np.ndarray, list], + y1: Union[np.ndarray, list], + x2: Union[np.ndarray, list], + y2: Union[np.ndarray, list], + *, + weights: Optional[Union[np.ndarray, list]] = None, + bandwidth: Union[float, str] = 'silverman', + kernel: str = 'gaussian', + template: Optional[xr.DataArray] = None, + x_range: Optional[Tuple[float, float]] = None, + y_range: Optional[Tuple[float, float]] = None, + width: int = 256, + height: int = 256, + name: str = 'line_density', +) -> xr.DataArray: + """Compute line density from line-segment data. + + Each segment is uniformly sampled and the samples are convolved with + the chosen kernel, producing a smooth density surface that represents + the concentration of linear features. + + Parameters + ---------- + x1, y1, x2, y2 : array-like + Start and end coordinates of each line segment. + weights : array-like, optional + Per-segment weights. Defaults to uniform weights of 1. + bandwidth : float or ``'silverman'`` + Kernel bandwidth. ``'silverman'`` uses an automatic estimate + based on all segment endpoints. + kernel : ``{'gaussian', 'epanechnikov', 'quartic'}`` + Kernel shape. + template : xr.DataArray, optional + Output grid specification (same as :func:`kde`). + x_range, y_range : (min, max), optional + Spatial extent. + width, height : int + Grid dimensions. + name : str + Name of the output DataArray. + + Returns + ------- + xr.DataArray + 2-D line-density surface. + """ + x1a = np.asarray(x1, dtype=np.float64).ravel() + y1a = np.asarray(y1, dtype=np.float64).ravel() + x2a = np.asarray(x2, dtype=np.float64).ravel() + y2a = np.asarray(y2, dtype=np.float64).ravel() + n = x1a.shape[0] + if not (y1a.shape[0] == n and x2a.shape[0] == n and y2a.shape[0] == n): + raise ValueError("x1, y1, x2, y2 must all have the same length") + + if weights is not None: + w_arr = np.asarray(weights, dtype=np.float64).ravel() + if w_arr.shape[0] != n: + raise ValueError( + "weights must have the same length as the segment arrays" + ) + else: + w_arr = np.ones(n, dtype=np.float64) + + kid = _kernel_id(kernel) + + # Bandwidth from all endpoints + all_x = np.concatenate([x1a, x2a]) + all_y = np.concatenate([y1a, y2a]) + + if isinstance(bandwidth, str): + if bandwidth != 'silverman': + raise ValueError( + "bandwidth must be a positive number or 'silverman', " + f"got {bandwidth!r}" + ) + bw = _silverman_bandwidth(all_x, all_y) + else: + bw = float(bandwidth) + if bw <= 0: + raise ValueError(f"bandwidth must be positive, got {bw}") + + # Grid + if template is not None: + _validate_template(template) + y_coords = template.coords[template.dims[0]].values + x_coords = template.coords[template.dims[1]].values + rows, cols = template.shape + dy = float(y_coords[1] - y_coords[0]) if rows > 1 else 1.0 + dx = float(x_coords[1] - x_coords[0]) if cols > 1 else 1.0 + x0 = float(x_coords[0]) + y0 = float(y_coords[0]) + else: + if x_range is None: + pad = max(bw, (float(all_x.max()) - float(all_x.min())) * 0.1) + x_range = (float(all_x.min()) - pad, float(all_x.max()) + pad) + if y_range is None: + pad = max(bw, (float(all_y.max()) - float(all_y.min())) * 0.1) + y_range = (float(all_y.min()) - pad, float(all_y.max()) + pad) + rows, cols = height, width + dx = (x_range[1] - x_range[0]) / max(cols - 1, 1) + dy = (y_range[1] - y_range[0]) / max(rows - 1, 1) + x0 = x_range[0] + y0 = y_range[0] + x_coords = np.linspace(x_range[0], x_range[1], cols) + y_coords = np.linspace(y_range[0], y_range[1], rows) + + shape = (rows, cols) + out = np.zeros(shape, dtype=np.float64) + _line_density_cpu(x1a, y1a, x2a, y2a, w_arr, out, + x0, y0, dx, dy, bw, kid) + + if template is not None: + return xr.DataArray( + out, name=name, + coords=template.coords, dims=template.dims, + attrs=template.attrs, + ) + return xr.DataArray( + out, name=name, + dims=['y', 'x'], + coords={'y': y_coords, 'x': x_coords}, + ) + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _validate_template(template): + if not isinstance(template, xr.DataArray): + raise TypeError( + "template must be an xr.DataArray, " + f"got {type(template).__qualname__}" + ) + if template.ndim != 2: + raise ValueError( + f"template must be 2-D, got {template.ndim}-D" + ) + + +def _is_cupy_backed(agg): + """Check if a DataArray is backed by cupy (plain or via dask).""" + try: + meta = agg.data._meta + return type(meta).__module__.split('.')[0] == 'cupy' + except AttributeError: + if cupy is not None: + return isinstance(agg.data, cupy.ndarray) + return False From 5f53eb0f069c426c335dc1cfd486c969ffb5cdc2 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 5 Apr 2026 21:37:00 -0700 Subject: [PATCH 4/7] Add KDE test suite and fix CUDA kernel cutoff (#1143) 35 tests covering correctness, edge cases, kernel types, bandwidth selection, weights, and cross-backend parity (dask+numpy, cupy). Removes hard cutoff from GPU Gaussian kernel to avoid box-vs-circle mismatch with CPU. --- xrspatial/kde.py | 6 +- xrspatial/tests/test_kde.py | 304 ++++++++++++++++++++++++++++++++++++ 2 files changed, 308 insertions(+), 2 deletions(-) create mode 100644 xrspatial/tests/test_kde.py diff --git a/xrspatial/kde.py b/xrspatial/kde.py index 753331f2..1f8a032b 100644 --- a/xrspatial/kde.py +++ b/xrspatial/kde.py @@ -8,6 +8,7 @@ """ from __future__ import annotations +import math from math import pi, sqrt from functools import partial from typing import Optional, Tuple, Union @@ -263,8 +264,9 @@ def _kde_cuda(xs, ys, ws, out, x0, y0, dx, dy, bw, kernel_id, n_pts): u2 = ux * ux + uy * uy if kid == 0: - from math import exp - total += norm * ws[p] * exp(-0.5 * u2) + # No hard cutoff; exp decays fast enough and each thread + # loops independently so the extra iterations are cheap. + total += norm * ws[p] * math.exp(-0.5 * u2) elif kid == 1: if u2 <= 1.0: total += norm * ws[p] * (1.0 - u2) diff --git a/xrspatial/tests/test_kde.py b/xrspatial/tests/test_kde.py new file mode 100644 index 00000000..7281a93d --- /dev/null +++ b/xrspatial/tests/test_kde.py @@ -0,0 +1,304 @@ +"""Tests for xrspatial.kde (kernel density estimation).""" +from __future__ import annotations + +import numpy as np +import pytest +import xarray as xr + +from xrspatial.kde import kde, line_density, _silverman_bandwidth + +from xrspatial.tests.general_checks import ( + dask_array_available, + cuda_and_cupy_available, +) + +try: + import dask.array as da +except ImportError: + da = None + + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +@pytest.fixture +def point_cluster(): + """Tight cluster of points at the origin -- easy to validate.""" + rng = np.random.default_rng(1143) + x = rng.normal(0, 1, 200) + y = rng.normal(0, 1, 200) + return x, y + + +@pytest.fixture +def simple_grid(): + """Small 16x16 template grid centred on the origin.""" + return xr.DataArray( + np.zeros((16, 16), dtype=np.float64), + dims=['y', 'x'], + coords={ + 'y': np.linspace(-4, 4, 16), + 'x': np.linspace(-4, 4, 16), + }, + ) + + +# --------------------------------------------------------------------------- +# Basic correctness +# --------------------------------------------------------------------------- + +class TestKDEBasic: + """Core KDE functionality on numpy arrays.""" + + def test_output_shape_from_width_height(self, point_cluster): + x, y = point_cluster + result = kde(x, y, bandwidth=1.0, width=20, height=30) + assert result.shape == (30, 20) + assert result.dims == ('y', 'x') + + def test_output_shape_from_template(self, point_cluster, simple_grid): + x, y = point_cluster + result = kde(x, y, bandwidth=1.0, template=simple_grid) + assert result.shape == simple_grid.shape + np.testing.assert_array_equal(result.coords['y'], simple_grid.coords['y']) + np.testing.assert_array_equal(result.coords['x'], simple_grid.coords['x']) + + def test_density_positive(self, point_cluster): + x, y = point_cluster + result = kde(x, y, bandwidth=1.0, width=16, height=16) + assert float(result.min()) >= 0.0 + assert float(result.sum()) > 0.0 + + def test_peak_near_data_centre(self, point_cluster, simple_grid): + """Density peak should be near (0, 0) for a standard-normal cluster.""" + x, y = point_cluster + result = kde(x, y, bandwidth=1.0, template=simple_grid) + flat_idx = int(result.values.argmax()) + peak_row, peak_col = divmod(flat_idx, result.shape[1]) + y_peak = float(result.coords['y'].values[peak_row]) + x_peak = float(result.coords['x'].values[peak_col]) + assert abs(x_peak) < 2.0 + assert abs(y_peak) < 2.0 + + def test_name_propagated(self, point_cluster): + x, y = point_cluster + result = kde(x, y, bandwidth=1.0, name='density', width=8, height=8) + assert result.name == 'density' + + +class TestKernelTypes: + """Each kernel produces valid, non-zero output.""" + + @pytest.mark.parametrize('kernel', ['gaussian', 'epanechnikov', 'quartic']) + def test_kernel_produces_output(self, point_cluster, kernel): + x, y = point_cluster + result = kde(x, y, bandwidth=1.0, kernel=kernel, width=16, height=16) + assert float(result.sum()) > 0.0 + + def test_compact_kernels_are_zero_outside_bandwidth(self): + """A single point at origin, quartic kernel: pixels beyond bw should be 0.""" + result = kde( + [0.0], [0.0], bandwidth=1.0, kernel='quartic', + x_range=(-3, 3), y_range=(-3, 3), width=32, height=32, + ) + xs = result.coords['x'].values + ys = result.coords['y'].values + for r in range(result.shape[0]): + for c in range(result.shape[1]): + dist2 = xs[c] ** 2 + ys[r] ** 2 + if dist2 > 1.0: + assert result.values[r, c] == 0.0 + + def test_gaussian_nonzero_everywhere(self): + """Gaussian kernel should have non-zero density everywhere (within cutoff).""" + result = kde( + [0.0], [0.0], bandwidth=1.0, kernel='gaussian', + x_range=(-2, 2), y_range=(-2, 2), width=8, height=8, + ) + assert float(result.min()) > 0.0 + + +class TestBandwidth: + """Bandwidth selection and validation.""" + + def test_silverman_auto(self, point_cluster): + x, y = point_cluster + bw = _silverman_bandwidth(x, y) + assert bw > 0 + + def test_silverman_keyword(self, point_cluster): + x, y = point_cluster + result = kde(x, y, bandwidth='silverman', width=8, height=8) + assert float(result.sum()) > 0.0 + + def test_wider_bandwidth_smoother(self, point_cluster, simple_grid): + """Wider bandwidth should produce a smoother (lower peak) surface.""" + x, y = point_cluster + narrow = kde(x, y, bandwidth=0.5, template=simple_grid) + wide = kde(x, y, bandwidth=2.0, template=simple_grid) + assert float(narrow.max()) > float(wide.max()) + + def test_negative_bandwidth_raises(self, point_cluster): + x, y = point_cluster + with pytest.raises(ValueError, match='positive'): + kde(x, y, bandwidth=-1.0, width=8, height=8) + + def test_invalid_bandwidth_string_raises(self, point_cluster): + x, y = point_cluster + with pytest.raises(ValueError, match='silverman'): + kde(x, y, bandwidth='invalid', width=8, height=8) + + +class TestWeights: + """Weighted KDE.""" + + def test_double_weights_double_density(self, simple_grid): + x = np.array([0.0, 1.0, -1.0]) + y = np.array([0.0, 1.0, -1.0]) + r1 = kde(x, y, weights=[1, 1, 1], bandwidth=1.0, template=simple_grid) + r2 = kde(x, y, weights=[2, 2, 2], bandwidth=1.0, template=simple_grid) + np.testing.assert_allclose(r2.values, r1.values * 2.0, rtol=1e-12) + + def test_weight_length_mismatch_raises(self): + with pytest.raises(ValueError, match='same length'): + kde([0, 1], [0, 1], weights=[1], bandwidth=1.0, width=4, height=4) + + +# --------------------------------------------------------------------------- +# Edge cases +# --------------------------------------------------------------------------- + +class TestEdgeCases: + + def test_single_point(self): + result = kde([5.0], [5.0], bandwidth=1.0, width=8, height=8) + assert result.shape == (8, 8) + assert float(result.sum()) > 0.0 + + def test_two_identical_points(self): + result = kde([0.0, 0.0], [0.0, 0.0], bandwidth=1.0, width=8, height=8) + assert float(result.sum()) > 0.0 + + def test_invalid_kernel_raises(self): + with pytest.raises(ValueError, match='kernel must be'): + kde([0], [0], kernel='unknown', bandwidth=1.0, width=4, height=4) + + def test_mismatched_xy_raises(self): + with pytest.raises(ValueError, match='same length'): + kde([0, 1], [0], bandwidth=1.0, width=4, height=4) + + def test_template_not_dataarray_raises(self): + with pytest.raises(TypeError, match='xr.DataArray'): + kde([0], [0], bandwidth=1.0, template=np.zeros((4, 4))) + + def test_template_not_2d_raises(self): + t = xr.DataArray(np.zeros((4, 4, 4)), dims=['z', 'y', 'x']) + with pytest.raises(ValueError, match='2-D'): + kde([0], [0], bandwidth=1.0, template=t) + + +# --------------------------------------------------------------------------- +# Line density +# --------------------------------------------------------------------------- + +class TestLineDensity: + + def test_basic_line(self): + result = line_density( + [0.0], [0.0], [1.0], [1.0], + bandwidth=0.5, width=16, height=16, + ) + assert result.shape == (16, 16) + assert float(result.sum()) > 0.0 + + def test_name(self): + result = line_density( + [0], [0], [1], [1], + bandwidth=0.5, name='roads', width=8, height=8, + ) + assert result.name == 'roads' + + def test_parallel_lines_higher_density_between(self): + """Two parallel horizontal lines: density should be higher between them.""" + x1 = np.array([0, 0]) + y1 = np.array([-1, 1]) + x2 = np.array([2, 2]) + y2 = np.array([-1, 1]) + result = line_density( + x1, y1, x2, y2, + bandwidth=1.0, + x_range=(-1, 3), y_range=(-3, 3), + width=16, height=32, + ) + # Centre row should have higher density than edge rows + mid_row = result.shape[0] // 2 + centre_density = float(result.values[mid_row, :].mean()) + edge_density = float(result.values[0, :].mean()) + assert centre_density > edge_density + + def test_mismatched_arrays_raises(self): + with pytest.raises(ValueError, match='same length'): + line_density([0, 1], [0], [1], [1], bandwidth=1.0, width=4, height=4) + + def test_weighted_lines(self): + r1 = line_density([0], [0], [1], [1], weights=[1], + bandwidth=0.5, width=8, height=8) + r2 = line_density([0], [0], [1], [1], weights=[3], + bandwidth=0.5, width=8, height=8) + np.testing.assert_allclose(r2.values, r1.values * 3.0, rtol=1e-10) + + +# --------------------------------------------------------------------------- +# Cross-backend parity +# --------------------------------------------------------------------------- + +@dask_array_available +class TestDaskParity: + """Dask+numpy results should match pure numpy.""" + + def _make_dask_template(self, simple_grid, chunks=(8, 8)): + return simple_grid.copy(data=da.from_array(simple_grid.values, chunks=chunks)) + + @pytest.mark.parametrize('kernel', ['gaussian', 'epanechnikov', 'quartic']) + def test_dask_matches_numpy(self, point_cluster, simple_grid, kernel): + x, y = point_cluster + np_result = kde(x, y, bandwidth=1.0, kernel=kernel, template=simple_grid) + dask_template = self._make_dask_template(simple_grid) + dask_result = kde(x, y, bandwidth=1.0, kernel=kernel, template=dask_template) + assert isinstance(dask_result.data, da.Array) + np.testing.assert_allclose( + dask_result.values, np_result.values, rtol=1e-4, + ) + + def test_compact_kernel_exact_match(self, point_cluster, simple_grid): + """Compact kernels should match exactly since there are no cutoff boundary effects.""" + x, y = point_cluster + np_result = kde(x, y, bandwidth=1.0, kernel='quartic', template=simple_grid) + dask_template = self._make_dask_template(simple_grid) + dask_result = kde(x, y, bandwidth=1.0, kernel='quartic', template=dask_template) + np.testing.assert_allclose( + dask_result.values, np_result.values, rtol=1e-12, + ) + + +@cuda_and_cupy_available +class TestCuPyParity: + """CuPy results should match numpy.""" + + def _make_cupy_template(self, simple_grid): + import cupy + return simple_grid.copy(data=cupy.asarray(simple_grid.values)) + + @pytest.mark.parametrize('kernel', ['gaussian', 'epanechnikov', 'quartic']) + def test_cupy_matches_numpy(self, point_cluster, simple_grid, kernel): + x, y = point_cluster + np_result = kde(x, y, bandwidth=1.0, kernel=kernel, template=simple_grid) + cupy_template = self._make_cupy_template(simple_grid) + cupy_result = kde(x, y, bandwidth=1.0, kernel=kernel, template=cupy_template) + import cupy + result_np = cupy_result.data.get() + # Gaussian: CPU uses a 4*bw box cutoff, GPU sums all points. + # Compact kernels match exactly. + tol = 1e-2 if kernel == 'gaussian' else 1e-6 + np.testing.assert_allclose(result_np, np_result.values, rtol=tol) From ed68557c3f7efd5b51194a130f97edbba180ab07 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Sun, 5 Apr 2026 21:37:28 -0700 Subject: [PATCH 5/7] Add KDE docs and exports (#1143) Creates docs/source/reference/kde.rst with autosummary entries for kde() and line_density(). Adds both functions to __init__.py and the docs toctree. --- docs/source/reference/index.rst | 1 + docs/source/reference/kde.rst | 24 ++++++++++++++++++++++++ 2 files changed, 25 insertions(+) create mode 100644 docs/source/reference/kde.rst diff --git a/docs/source/reference/index.rst b/docs/source/reference/index.rst index 3b3eef32..a268b5cd 100644 --- a/docs/source/reference/index.rst +++ b/docs/source/reference/index.rst @@ -16,6 +16,7 @@ Reference focal hydrology interpolation + kde mcda morphology multispectral diff --git a/docs/source/reference/kde.rst b/docs/source/reference/kde.rst new file mode 100644 index 00000000..0f832d6f --- /dev/null +++ b/docs/source/reference/kde.rst @@ -0,0 +1,24 @@ +.. _reference.kde: + +*** +KDE +*** + +.. note:: + + Kernel density estimation converts point or line data into + continuous density surfaces on a raster grid. + +KDE +=== +.. autosummary:: + :toctree: _autosummary + + xrspatial.kde.kde + +Line Density +============ +.. autosummary:: + :toctree: _autosummary + + xrspatial.kde.line_density From b6aa0d75029fcdf633824a57172d27ee24801707 Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 6 Apr 2026 05:34:10 -0700 Subject: [PATCH 6/7] Add KDE user guide notebook (#1143) Covers Gaussian/Epanechnikov/quartic kernels, bandwidth effects, weighted KDE, and line density with synthetic earthquake and road network data. --- examples/user_guide/49_KDE.ipynb | 177 +++++++++++++++++++++ examples/user_guide/images/kde_preview.png | Bin 0 -> 86095 bytes 2 files changed, 177 insertions(+) create mode 100644 examples/user_guide/49_KDE.ipynb create mode 100644 examples/user_guide/images/kde_preview.png diff --git a/examples/user_guide/49_KDE.ipynb b/examples/user_guide/49_KDE.ipynb new file mode 100644 index 00000000..3266e327 --- /dev/null +++ b/examples/user_guide/49_KDE.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bfd6d891", + "source": "# Xarray-Spatial KDE: Point density, kernel shapes, and line density\n\nTurning scattered point data into a continuous surface is one of the most common spatial analysis tasks. It comes up whenever you have event locations (earthquakes, crimes, disease cases) and want a smooth density map instead of a dot plot. xrspatial's KDE tools produce density rasters from raw coordinates, with configurable kernels and bandwidth.", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "ee7bb1d4", + "source": "### What you'll build\n\n1. Generate synthetic earthquake cluster data\n2. Compute a Gaussian density surface with `kde()`\n3. Compare three kernel shapes: Gaussian, Epanechnikov, and quartic\n4. See how bandwidth controls smoothness\n5. Use per-point weights for magnitude-weighted density\n6. Map road density with `line_density()`\n\n![KDE preview](images/kde_preview.png)\n\n[Gaussian KDE](#Gaussian-KDE) · [Kernel comparison](#Kernel-comparison) · [Bandwidth effects](#Bandwidth-effects) · [Weighted KDE](#Weighted-KDE) · [Line density](#Line-density)", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "92a81141", + "source": "Standard imports plus the two KDE functions.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "a1764753", + "source": "import numpy as np\nimport xarray as xr\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Patch\n\nfrom xrspatial.kde import kde, line_density", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "965e4d4d", + "source": "## Synthetic earthquake data\n\nThree clusters of varying tightness and size, plus some background scatter. Each point gets a random magnitude between 2 and 6 that we'll use as weights later.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "78a82a28", + "source": "rng = np.random.default_rng(42)\n\n# Three clusters with different spreads\nc1_x, c1_y = rng.normal(2, 0.8, 120), rng.normal(3, 0.6, 120)\nc2_x, c2_y = rng.normal(7, 1.2, 80), rng.normal(7, 1.0, 80)\nc3_x, c3_y = rng.normal(5, 0.5, 50), rng.normal(1, 0.5, 50)\n\n# Background scatter\nbg_x, bg_y = rng.uniform(0, 10, 30), rng.uniform(0, 10, 30)\n\nx = np.concatenate([c1_x, c2_x, c3_x, bg_x])\ny = np.concatenate([c1_y, c2_y, c3_y, bg_y])\n\n# Random magnitudes for weighted KDE later\nmagnitudes = rng.uniform(2, 6, len(x))\n\nprint(f\"{len(x)} earthquake epicenters generated\")", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "3ff48b8a", + "source": "Raw scatter plot. The three clusters are visible but the exact density pattern is hard to read from points alone.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "ca17dc8f", + "source": "fig, ax = plt.subplots(figsize=(10, 7.5))\nax.plot(x, y, 'o', color='steelblue', markersize=3, alpha=0.5)\nax.set_xlim(0, 10)\nax.set_ylim(0, 10)\nax.set_aspect('equal')\nax.set_xlabel('x')\nax.set_ylabel('y')\nplt.show()", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "498b0083", + "source": "## Gaussian KDE\n\n[Kernel density estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation) places a smooth kernel function at each data point and sums the contributions at every output pixel. The Gaussian kernel is the most common choice because it produces infinitely smooth surfaces with no hard edges.\n\nThe plot below shows the density surface with the original points overlaid. Brighter areas have more points per unit area.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "ca293129", + "source": "density = kde(x, y, bandwidth=0.6, kernel='gaussian',\n x_range=(0, 10), y_range=(0, 10), width=400, height=400)\n\nfig, ax = plt.subplots(figsize=(10, 7.5))\ndensity.plot.imshow(ax=ax, cmap='inferno', add_colorbar=False)\nax.plot(x, y, 'o', color='white', markersize=1.5, alpha=0.3)\nax.set_title('')\nax.set_xlabel('')\nax.set_ylabel('')\nax.set_xticks([])\nax.set_yticks([])\nax.legend(handles=[Patch(facecolor='yellow', alpha=0.78, label='High density'),\n Patch(facecolor='#1a0a2e', alpha=0.78, label='Low density')],\n loc='lower right', fontsize=11, framealpha=0.9)\nplt.show()", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "1612b72d", + "source": "The three clusters show up clearly. The tight cluster in the lower middle produces the sharpest peak, while the spread-out upper-right cluster has a broader, lower dome.", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "d9a7148a", + "source": "## Kernel comparison\n\nThree kernel shapes are available:\n\n- **Gaussian**: bell curve, infinite range, always smooth\n- **Epanechnikov**: parabolic, drops to zero at the bandwidth radius, statistically optimal for minimizing mean integrated squared error\n- **Quartic** (biweight): similar compact support but smoother falloff than Epanechnikov\n\nThe next plot shows all three side by side on the same data and bandwidth.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "cff98669", + "source": "kernels = ['gaussian', 'epanechnikov', 'quartic']\nresults = {}\nfor k in kernels:\n results[k] = kde(x, y, bandwidth=0.8, kernel=k,\n x_range=(0, 10), y_range=(0, 10), width=300, height=300)\n\nfig, axes = plt.subplots(1, 3, figsize=(15, 5))\nfor ax, k in zip(axes, kernels):\n results[k].plot.imshow(ax=ax, cmap='inferno', add_colorbar=False)\n ax.set_title(k.capitalize(), fontsize=13)\n ax.set_xticks([])\n ax.set_yticks([])\n ax.set_xlabel('')\n ax.set_ylabel('')\nplt.tight_layout()\nplt.show()", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "69773ca7", + "source": "The Gaussian surface bleeds out to the edges because it has infinite range. Epanechnikov and quartic cut off at exactly the bandwidth radius, leaving the corners at zero. For most practical work the differences are small; pick Gaussian for smoothness, Epanechnikov if you need a hard cutoff boundary.", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "d3e9f621", + "source": "## Bandwidth effects\n\nBandwidth controls how far each point's influence reaches. Too narrow and you get spiky noise; too wide and real structure gets washed out. The default `'silverman'` rule works well for unimodal data but can oversmooth multimodal clusters.\n\nBelow: four bandwidths from very narrow (0.2) to very wide (2.0).", + "metadata": {} + }, + { + "cell_type": "code", + "id": "4fee443f", + "source": "bandwidths = [0.2, 0.5, 1.0, 2.0]\n\nfig, axes = plt.subplots(1, 4, figsize=(18, 4.5))\nfor ax, bw in zip(axes, bandwidths):\n d = kde(x, y, bandwidth=bw, x_range=(0, 10), y_range=(0, 10),\n width=200, height=200)\n d.plot.imshow(ax=ax, cmap='inferno', add_colorbar=False)\n ax.set_title(f'bw = {bw}', fontsize=12)\n ax.set_xticks([])\n ax.set_yticks([])\n ax.set_xlabel('')\n ax.set_ylabel('')\nplt.tight_layout()\nplt.show()", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "9dac2c7a", + "source": "At bw=0.2 you can pick out individual points. At bw=2.0 the three clusters have merged into a single blob. Something around 0.5 to 0.8 captures the cluster structure without too much noise.\n\n
\nUnits matter. Bandwidth is in the same units as your coordinates. If x and y are in meters, a bandwidth of 500 means 500 m. If they're in degrees, 0.01 means roughly 1 km at mid-latitudes. Always think about what physical distance makes sense for your data before picking a number.\n
", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "2465dd33", + "source": "## Weighted KDE\n\nNot all points are equal. An earthquake of magnitude 5 releases 30x more energy than a magnitude 4. Passing per-point weights lets the density surface reflect importance, not just count.\n\nBelow: unweighted (left) vs magnitude-weighted (right). The weighted surface shifts toward wherever the larger events happened to fall.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "8cf1850d", + "source": "unweighted = kde(x, y, bandwidth=0.6,\n x_range=(0, 10), y_range=(0, 10), width=300, height=300)\nweighted = kde(x, y, weights=magnitudes, bandwidth=0.6,\n x_range=(0, 10), y_range=(0, 10), width=300, height=300)\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 6))\nfor ax, data, title in zip(axes,\n [unweighted, weighted],\n ['Unweighted', 'Magnitude-weighted']):\n data.plot.imshow(ax=ax, cmap='inferno', add_colorbar=False)\n ax.set_title(title, fontsize=13)\n ax.set_xticks([])\n ax.set_yticks([])\n ax.set_xlabel('')\n ax.set_ylabel('')\nplt.tight_layout()\nplt.show()", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "9fdd26fb", + "source": "## Line density\n\n`line_density()` works the same way but for linear features. Each line segment is uniformly sampled and the samples are convolved with the kernel. The result shows where lines are concentrated, useful for things like road network analysis or fault-line mapping.\n\nInput is four arrays: start x, start y, end x, end y for each segment. The plot below shows a synthetic road grid.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "eab23452", + "source": "# Synthetic road segments: a loose grid with some randomness\nroad_x1, road_y1, road_x2, road_y2 = [], [], [], []\n\n# Horizontal roads\nfor yy in np.linspace(1, 9, 6):\n offset = rng.normal(0, 0.2)\n road_x1.append(0.5)\n road_y1.append(yy + offset)\n road_x2.append(9.5)\n road_y2.append(yy + rng.normal(0, 0.3))\n\n# Vertical roads (fewer in the east = lower density there)\nfor xx in np.linspace(1, 5, 6):\n road_x1.append(xx + rng.normal(0, 0.1))\n road_y1.append(0.5)\n road_x2.append(xx + rng.normal(0, 0.1))\n road_y2.append(9.5)\n\nfor xx in np.linspace(6, 9, 2):\n road_x1.append(xx)\n road_y1.append(0.5)\n road_x2.append(xx)\n road_y2.append(9.5)\n\nroad_x1 = np.array(road_x1)\nroad_y1 = np.array(road_y1)\nroad_x2 = np.array(road_x2)\nroad_y2 = np.array(road_y2)\n\nld = line_density(road_x1, road_y1, road_x2, road_y2,\n bandwidth=0.8,\n x_range=(0, 10), y_range=(0, 10),\n width=300, height=300)\n\nfig, ax = plt.subplots(figsize=(10, 7.5))\nld.plot.imshow(ax=ax, cmap='inferno', add_colorbar=False)\n# Draw the actual road segments on top\nfor i in range(len(road_x1)):\n ax.plot([road_x1[i], road_x2[i]], [road_y1[i], road_y2[i]],\n color='white', alpha=0.4, linewidth=0.8)\nax.set_title('')\nax.set_xlabel('')\nax.set_ylabel('')\nax.set_xticks([])\nax.set_yticks([])\nax.legend(handles=[Patch(facecolor='yellow', alpha=0.78, label='High road density'),\n Patch(facecolor='white', alpha=0.4, label='Road segments')],\n loc='lower right', fontsize=11, framealpha=0.9)\nplt.show()", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "43046e13", + "source": "The denser grid on the west side shows up clearly as a brighter band. Intersections where horizontal and vertical roads cross produce the brightest hotspots.", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "d1036afc", + "source": "
\nTemplate grids. Both kde() and line_density() accept a template DataArray instead of x_range/y_range/width/height. The output will match the template's shape, extent, and coordinates. If the template is dask-backed, the output is computed lazily in chunks.\n
", + "metadata": {} + }, + { + "cell_type": "markdown", + "id": "5dd32bc7", + "source": "### References\n\n- [Kernel density estimation](https://en.wikipedia.org/wiki/Kernel_density_estimation), Wikipedia\n- [Silverman, B.W. (1986). *Density Estimation for Statistics and Data Analysis*.](https://ned.ipac.caltech.edu/level5/March02/Silverman/paper.pdf) Chapman & Hall.\n- [Epanechnikov, V.A. (1969). Non-parametric estimation of a multivariate probability density.](https://doi.org/10.1137/1114019) *Theory of Probability & Its Applications*, 14(1), 153-158.", + "metadata": {} + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file diff --git a/examples/user_guide/images/kde_preview.png b/examples/user_guide/images/kde_preview.png new file mode 100644 index 0000000000000000000000000000000000000000..53bda77c3a64204a769f4bc4536df2c9031ec05b GIT binary patch literal 86095 zcmcF~hd^^d6?=N!>h536Qn*RFM{uW+N{y~p?9l4Ai`FneK`FprNz7Xi>>*wm_ zd0k%dn!KXy1s8vRZ$C{11^B=JA@Akutl;OhJH2S>enoi5_-HG2?+w8rMmI#j5??nero<1vL zyDj3??@gWG&eMbYJDzFo(tdB*e~s8q{Jl~E>$}AMb1w6DZTD|(u4ghYIOn;(zvwLz z=rPn7=sW%Jp*+{{gNv_T;|<_3LS(~#l!VA!&q|(!h;X%W$DZimXBS?4D*#z;5E%73 z(#gIgRQFK%N)W}P?~*VVSI_Zq^@Hyvtm~tm@yw_7G3vsE8n>(b247g!c~&g*Fu5jf z?o%btjf8T25MJ&E06U%245E6Uul!;#r|gTm7>B~jiyT<@hd2`k(@G!qHQ>OPnjf&^ z|K0jQ?VzMEm-oq@()Q4?wbs$^kq5`4=Iy*vt(?kn?X2g_@c`41O_kT7Fs=>m*u95z zXmlppZmZLcXEo}#@_uQRVGZS>71m;>7sFg@mi48s54ZlqZuXh58h8A^IwG|vuN0+k zr;QHw9&gKGa!he=Cezy_zK$e)Z{uoH4#X%PvWqv+b@fW#6PHpb!?gB~_c5 zPo$=pQ~;!(=&}R+IXdqAw>3bg5TIJ zI+6y#RB~XOwV%QP+1V-dC1azpdiWnz=f2*b+WAXI-2fK)Q#$2Roj6-8PY@0{s( zgtV+4%qwE?G*FZtoDC3u$*-5riBd|e9;T{~gz+jr9198WsgvL5BZq5G^4HE@zJV8G zZcfaCR3q_iT#_}J19oafX0E-8y(dw_c%-<6)-)*>aT>1NP@loie=j^#c`dFw5)~>BZrf?1i6s#)VlO-rvak zY)=5cd^f^|f`T1cuO>5t9zR6|)^GCT+Xa4y13+o(~mbluxaj zOZ10ArTnzO!ns`!AN2Sf{?{Fe{yA--eSOvIIM~a;GB~)4e{Okr3DLZ=CJKBfHWlWM zt=)Jyzl}pp#sj><7%ZR7+SNL~&06QlAH7RthRL57DO|WC+{}s{y&W!k{G%N+FHmtu zWreO5M$rCd*tc@$*V2}st&GmHpYJiL&a~QTMu#nF>-oY=ZJ5LdvW-m>9brl=5+&QhI||j+&B~h%ivb<4B{jOuPf}7NKBj+_1rqs z%;Fo?7AoT?ST;?snFnVkTyvamm^;OQn$AQ=$k#ehoYM4OFIh;{<%qFOaWh%p6aK{7|19-q$y%aJP)nR7j`9}8?o4yN*`Of2 z@2_s4he`EXZBcCqzD{ia^zuVMPqk&xkKlIp&+{m{AB>if-D6oIE*D>!Gafkon1X)d zrLW3jf7+~aHg#rOK4Fh$moM}-Tl!(~fN~Q$7zAmv)?N*0D1>FK1-x#hu-fR^3u+X# zOB4B76c(Y0z8_pd1xxzY^oh z;gTYQp|Fqi?;~w1CNW#yn0`GbuLWbm`|*6}d7Qq1&ZgKuebi!n+{y6K2!VW=YB?tD z!L+u`NYUdIZ1eXNUT4?gN_=~N=FO72P*E*p+EXmYKD(uKTUDtDxIw#))L;!iwLfes zSR(Sx4)#2Vp;fO&z7HN*s!wVBH=NuOK#pH3%b84Ag|?)(z0MxTX3E*Kd1;6PA*tVn zaaY8C`Tv>&%!8nj-f%Ly9}-Q;n4YRdQ3F~4iJBhM z_2{=24eT&bqkP2*8@i}sUZN!MOG5;GT58#PDtszV=U4CU#fG_qzdT~;4eC-45YC6E z=On({dno%UUzR^+UIyU~8vHyjmHXKl@$kzj_gt`OOug^!X5N?Uh2@MzU%RHO$TXKT z{ZPs>n-zQgAM)9%=6Dcks9q;J)LiA!21;^pC$pK?<9{ zylnn`ZzoNCtwhX5ny(*yCxfjEk*PI+G37uh+P*%$s0ghWZXx_AYRgI+Vvx+xQxwJC zF^TP^_5t9>*9npkW|8#AAq^Qbgf^fU1#$OuPcDbJzX|Yfl>ZvG)VopF^Y~Nm_Z}wA z2M%}qRCR3snVG?tj^b89KETf4xZT4&;R{TsI zPlv-*ECVXq`7g$;zx+KXk(Y0Md{HDD&o?m-o#3GOR!L3G&4JsI3pHzXFAyp2OLlnjaBN0h3kjcFMi9@mBs)yv%2j$#p-6lF=!sNS(G#;o)0t2SZQYEI z;r_}I^s|&3EQS5#mn-Luf|ajbf{lEw|Gl5)lwLRa6ZS`w&!UofDd@2Xp+?f_+D2tnQM$ss@sbw$w>*@Kc z0RfJI+Uasxwlw%|2-jg)s1hL%4Iq>?O8T!tXw*S@P561nr^6XD@c8 z=KL^ljqaT`jvWi+d%K2;_F43bmbjTudoIVAxbLp81PdZE4{j9e&;PjC-T8PRC{w}V zE1+$Vk!{I0aZ)!8UvZ=Lc3?!KmrrQcg!zzjxep!mJ|W;ZypD^p6twZbecLkybKkb3 zK{t|}<@h?*cv7=rz2RF8O;#p6qX5edJ zi1>F+3Cn-HhEx?3Y;~O70ZT~a@X1FX?!K*dvQ{ydsw(`nshjqC{f$%NeXvqf@&uSM ze&vkq<_{@1g40~_ykXeZ|bA#p{c!Cz*=GHdxPaKX07)gK6ciHulC*(qud z+XMGR2SuE&=S-pO6x3?lI5s6+zQCRU?8U)*0C)*|tI^``&Q4QN*jXLH8-S(GEo^S6 zm<^&X2(x3hS!<_aGsX^8-j+(1V}#O&8A+F6UAK}BM^{v(!$q>F1FA8ntoD=GXPwJD zK$e-?j8E_J9^9mDTrTU9Ns!ASe-CQgjp!F4s)WEhC}Olyh<^0kiP7S!iyl4z9Xq_7 zMJgYAVmWNV2WCL}MnNipDv*1^)$yv*73Dj#4RkKF6HP8sES-^&&ncYf4tz=9E=5g#Wy` z?hKK85+eUFMAp<$xX|sM?{xi**zcn}-W|;}rP9hayP3huNDUcbdxO` zt|z4?&TTx!>A$ux68z?y1B7y+8MPm(p@*c%-6I|9Q?Y{3j@=mfW;C1#c+yb5m>)M8 z%p$4YsmeV!a!+Ht^gYW1IOs0mQ<@Ot;l*Rwe>JqO2?%+7Y+7v0FCGaf_)J;5?bv6_ zqOHk!yc4vTemA8n%%XWbIG^1(%EMCI6l)x++XlZyeWkg;o5V zSrifh2ei`&t+L%TSYboir%j{6o98X>+v4g&e10;qCCObAiSam3avmhyPNBZNLHvl> z2Db&k`wB(Km?Pg~D@#eR9b`4Up$C@*%i6McKN13x+!D0E98@KWd0!l>GG03Jmf#^$ zcFnO5c)W5wS~>a+tqr5EFfjqvOQUUdDl&$bzFY+^4dlzmK&mos)9$w577Iby8no|P z`>cQ901P`~L^!`yLNRyH00H_;~+ZQ)quUn>T2Z!Y5q z-GhvGz$=UK;C_P^^H5sCpjYmcRS>4KPLq_7+a(zgm;z6^+%;k&L3Os)_pm^yY?LV) z#}!`)!WuwDDB}Z2w*{?p`;RyMz=Ims2KP4h9?6;XJWM8sCuUhd&8OYYx2+~8)|f{y{m6IGv}Ht)CEN*^z1(7DMc znH>zE{y=Zmn)ex*HrnAs=>&(bEw{?~^4p`2r}w2<0j}^Syz=W|LdNftNs55AUFZsQXwjDE1qigk67FjIW zd5$=nx;HTM?!1HBFmQ(<>sp4>RC>(LTvAu&g|LNIlHdpmb@wMrb!lcXR?jaKEn7YKk z17i*qwL!CQ#oJ6*cT-SLnb=6izk-O4yOp$J52QLbN_6w3bXX- zbLQ4=(NPnMKa`UYMf+LywNLC!`8(Z|PVlsYUXnHr?nMYJY1VlBkM~NrYIq|MJX$mH zl4|RI%&XYvI4bboX@Egt5@rL9t2d`poByF+b2Q0Fd8O2Lk-UObzf?ArwqtqoyWMx3 zV??%KLqiLq6T1^iNGnqNG_GO^m-2g_$6*|Y6@?$FkKOknTmiQ_b^;TWmawwt=uc>M zNWpJRLEL(KoG0;dU)3=)6CRQQ11>||&_-V+)n?m{>88>3p3_T?S(>1T{t|^>2k+)f z8EO;Cyb^&BT$ccLD)?ICkjgp|bbImXb^Rf1DJAKK(Q5ULx+( zEn45?udFkqAZK|nE<{s`=n#Si}wyK~8M*hkC_UkK};8`i54Hy&R`hVYl+3 zU=YpRF>tBhLR_`?o*B&0_o+6b?&8#(=Q*0ARLN9pXkRGIC1CXmL3DEmb^F_aV*Mo!*SxW@@Ch?jx7DfA}MpkYF^_c~Wa~%#SVuB#CjHsZ3agz;> z4-)=27z3UsbKr4j3S~z$JP;Dx0It?)$eD5c#f}i^XKoxyNT9oZVEI3h_%-p{c-U4d zmZ(j&wYul`C@q)n^Vg1blQ`6%C`y^x`k?DB}jW+hGZ4udS*sYzL9X`~Q|F4FH z=S&GbcRt7^N?{}&say22*DE{~JQWt)2Gz=9AV~OCcmpnv8M=V+_H_I!UpkK(qHLZ3 z$J>tOcYR?!eNlFdK!Z@)s_e1F!Pmv^JJMzB_b3bYTpS9euitPge?;y7UFAa^eVW}W zFP0zJd~~!14{Y@I?Ukx_OEq7RxEB>mF8f~K z)&W!Bk5Bze(|EZ0E?2XG)yipe*vSoLROOkL2GCgEXdQp~tTw2@T;DQ(VS5lu=r|jX z9BZJC;@ya;oS!3Na{(nPA%CY7txxdGoK&A{-#FsR5~GL;f^k(Pkw>)Vz7L_rzYP}G zrEA)_LopTL6*;g|NYpKoIl|EJfE+SXscz<*3Z|B>cO+%6`Eh8^A^*?%17GZc>&6KN zk-iz)Dij*sO049%g(gL>ppjIQ zcMG|z@*Ma4jqr5Cj+=rK;K{amvV%DtqJZW*V=nG22Rr@xTKnOSC;g zdK;^fuw+fqP4CoA3o8U=otUwF)H%VH*^jv&9G@pc{?ay6U#U6}?5p(f$lqKOpVhNJ zHt6PsAdgnY6|ROMr3p_5=9Hm_+{Cn*i{5f~Lt>Ie8U^A6L>anm45KOY zR!=y^qDGf$A zG$m*>5CG+4PjNZ+uk$+FJTX7Q#Y;r@x-9=)d!3|lFUL>g&>8^uHPHrG0R(jn@Q!JD zc-ovNsKj?(MjI3?Fe+6^0F-K4epUQ6{!*=0)~l;x=SDanau)dBK52cXcf;GxGqS1v z&F`R>L#rpLB?yUxdKkR4K;N;^^H5NXkk?Bz^E?>Ra#W*Z8Uuw&+FSL^e;(AvPTFkN z&WhfwlPRBGblVwfeY#|i=oJMN-rQ&ec)>gJlIvlKc^N&kYZ^ab=YE8h-PoPyvz`(- zNPn2PrK^E_DMZo_*i?P~S?h;h;`_e*S(_`xkz`QEH=;opgTXgq$v1?5R5sl<5`NKr zeav>INO`;F?}ZxG>XEvphKk19Rr;pOK=YZTt;6zg=#M98e4263LHP=(9 z%Cp zzktE;p*_Jy>&T5<_ZP|S$dM72Vh0>_>n;T9AxY??0Deb?n>@97+f8KRwE5QN_D|)& zIc>`h4$T}-*k~ejnP6=|GK3+^SETP|F^2iY=t|xB8r>wFjV}=7IZHUdqaBdr@~9uH z0X`C&2k6EQl9zXmN>a~LGEw5bQA*A`e7flfoq5g&rLyGj-g`v^2J`fB)4>4J;CqT_ zU32D~_TlW?;2!7eILk@l3;$JyXPw?BC29Aba)papZgsN{e7OJG1hrD?9b}v7xB0U+ znhe4Q*MSLgV0U*Eby4!4twj>Xc7B#A$G$9RoM`JFZ&UVHMB(7+sqm7W1{Ikq3v8Y{ z!9F|I?qECR3I2!IozJcaQu-ys=t+PLav_9OlzEWc^a%=_IULkOgINUaIBez4 z_5&ZfeEl;tKVLM4-k4tEp?yfEDA0?DtIC3o3$@L6(BRlfPq1eQaw_R3$SGxmKaq5Y zOdYg!j}Q-f0)wwJN6&WiMiHIy_e`TyP)&?P%ypF zBV%(usjv*iRk_u%TR#DBut#%uu^`%O0S;pch&s^{Su+nmQvolmR+RyLU8czd4ryByS%MevT_;nl^5$UYe$ zSK}ci4VDu9^$E2#szUE2W2ub9afvl?Xkt!$$B&DoY~jYhdzBLffLUjxxRPG7!!Jd9 z>%yIW{i0)uPGd8@DCMI5)(10>UnynqjVv(IP7AlSzUZ`l4~P81IcK#w zH%>h6=t&r+woR_V>xgHWivPBo-=CBsbL%})i;kA#=K8=jh$~8eJjtuYaFNGUq`sKH zUg1`r<+P*u_7Q>!b302I>OCF4?F#9BJGubD5wh%s)@RJe-ao*Kd(7C#$PLmktg>k3 zdCWgIqm?o+(aP8JhZ9EuFP9z~uQQBr1%c7miwNXi6~)Wg5RsC4XKV;#n+Gyy1rvNwnC?TB>+s#Gz_jq@6Qx z+~I9WG!&{V^Z`4%_|9If`N^E8S){4{eIS@3^^y+{mi*}VACl1lc|dT0hhm;S2(_^k~NMa7Z6?P zNrN+C>bfZti~V7#QCDby8Y1uLO(Z!w4TTvBZXI5t{G*ZW$d7hTox?`&`K~Z=o z5U!SK>aJpRb*PQH6^n4=|NEmH?NkBz_y|c~3b_S&-8n1iTu_Bi5VROaAKC2w&e1e_R6c z#Z%#jO`ALtIzLB#9VC6s4j*EsZY2JtaXqV2TiZ9Jfn0-Mma5SPkC`|;oiSB1?r(f) zTlaI!Npu?g0wQAIdcBj%$h0L4<#%oDVpiTVABOv#cUv;lDz9wlZ|xdGk;_hdu{G?i zn(XXLe&zdAuzxFH!0B3Me*R>G8)o(Rlu)*3dA z4+2;%Ai>UKC2ng0spKv|Vcrinhm9^89ruHkTRZp5m#*EgUyVGo;Fct5Z~#iC{=$<5 z$9!)_h!+)n`l!)4@nB_^y|O0ZDlbIFnW|2p9beoKpkR1C;To5EcY&=dQq{W}C@f)! zxG4EbbSHB<{Kh6R4Y=S^r8+bkuA+RHp5-Fn_6oakPJ_Gx$+_mG7E=CnMBJc}m6tWM zca~(d`PC_96E^0~iTnCi8FjP6rMbl|S;=~@SRb}WU}Wb;`8H4--<}z4epZeA176VC z*1GZXrv3!-AX<2}Yj&XC?ZYMZrO+yAK_SvZ;WO9Yk|6})-^?uLcFjlG zY&=#IR;WDT)SW3soGw3VfrT)_N22=`SgG*S>j9g2jz<$WqLm;spxKr9(-cu=bzd=# z`16EX5$@O6U(H~$GOseLx_jo`f*hVZRd{$uF6%TsZ z4L6D-Bw+NVR=0-7dOIxrq)m?eM6ayvtY3d^yTt~>7ZpW!;!US4Tg~jK#5;Hu!39b- zFT>eD?Y<<5WswV~Um{BUM_UBpbJ!<(HyLVqB(aTa7F86EuA zoSJU7E-mAK9Mcrvvl+V&F>{&Lyicy%{cCZo#=MbLGq)JvckWsCeZ-~4rQKeyJQYoS zfLCs_Tf(%dcStb_?WSd(fG0GaqFSDJC~+%UY)kR#Y^!J_V-A+c^xlnd@ZdN2d{Gc# zpaeMhTCB_~*l4wZa#)25Q#@L>5WE^db&oyM|DbO=#v`j-X3gLJwnc*D*an*emD=VQ@S5l0b{Z7~>)i$d1*m{pA?-8t z{68}~|Fb5Qpw=B?ufsZdZ8BBtgP=UZ)p=eTvi7WH@o6txvaFX#zTo9<47z`BQXzT1 zb-rN1$5bBd?T#G7YY>+dw`52RUmf=yz4x$WcgoJIhs-gGqJMyO9L*HJE>#5s9u<*9 zsQ<*8qb9nht+f>9=R7xdE8aa_4Zv{w9a#uI8N&+6@F-Wm$*0fWvB&zB0zT}S`55J0 zss)wmYd19b2>Zirk$Ot^?nu^y1A*$-I#(QZuyaIN?ds^flK4j&3XskB+CNe<<`2h? z!$h1b46S3^@`XNFe|eGXbJh-i9Kv${;V4(lmIn`NeyW?U$M@gHJZh6%40&rf>p21!880h$QB<_4jGqNFuihoDg4o1x;WK;freD?Q1l`gQYMjfX zou-s#PRb%}Y`RQ*`qE7pN|<2jrof~Lr6qo6w55GYINUTxKs~HvQG%e)|o%EGLu7m3z0(w#$5&5 z_ue-=N;tBvA~+(DOM`%|g407uAaC(3Cn@aZd~BIeH^{@;#W(&1#PdgWA}~`#(`JXC zsw&pvrZDOWf=?_y$;jBMCH`iMB87G0c3#VW?YMGgw~nC;4$;4Nt;yle#i3kQ|L}d6 zip$8fGxRB|cSjtJ%M=yu7uAy5kN^##q`CYw^FDSy4GLt@>Wms&mm=K@L8HT#QdJ;- z$2#`#-Fv<-O=nI;UKg1#^DS0=g?O^yM=Aq&O+%Y__uVbmO9fM30|iNCyDU#>nz{1< zZ(L5DbUU-&pjy@$ly@P~GYXv9IWKj3>Ut@_!V)q~-vGgN;vGfq$?YJ+Qh-H^ZAmGH zJKDI#x_4#c!Qin|1frw#MR#M?>MdUZMPX%w&OJa&dd9pO|KWX1;$^U-oV^vGcJ`*P zA0Sa0C3_v&xM9<0G$Xs_PnJ%4B(Ja6;AzW{Q&IHwHA1PXi1u$O!2KdY7J4Oyhp8JS zT?6qBL8kuhoqW`Dw-5<_GLW^fJZR#qcnI#(HG zUHK2OgW6(mCpexjl(8u12%X$M<)*P6}Q zA2Sphct+a7PYo6X&O^rEWoo(_&=VO}lx-oC4Y4|`4?4B4)J za^hWGhsiuStF>5?maE87c33$H^O6gyi`(0vj*nhT13}F5|hqI=&9k7IZb%a?+%bnqHig2 zFZ1>QBS_!Cuh-0P0?)&@rRnII!G5yPj5#pE$OgtJksfnOGuH=ONh?fsfbGh-oty+o z^v^QS30N5gBmw>&Zq&!0+l8c10;oT7e+-$I{9rvtUxsN9eLMuRp2~`;h!{Wk<($&C zlxohLX$uj&IG+4|`UEjrp5zbhX`8M4CG%0Uy!`A|;P_+mkUFj)d!_)|G@<&nhw)}= zxbchRrtHNILaYN^7T+7vl^NP}GdN26v>nM&Ll~b;SSIMEY0rRarKwxH+JzwB3J!}s zg0QKTjL;1alx0p=Pt?Bj^~is9{`BajvQHZ1Gi064FIayc53}<`sA1x{T7AxhH7StO z*y;JeT6*8}{E@9mr~33AE_F4rOMmn3MMh$A&KQeZJ|LZ>V1^J>u^p30vm56_en;b0 zjsU#m93S;GE>E1}v65h}iS}6Y7OVbEez!-ID!(C0UZm)7`WE&pua$ghNu5wkW+TnD z!L}1gDS6-e?D1zLFzQ-rs&F2B{%cESS6mFe3L_%e#B_54b*eVmmJslA3Qt%?7q;2w zJNXvDT2_??Qr9i-?W@5bbZLNFU-+4}G+0M%sTY$~HgVW3J_pDlnA;c%dImR#S#@D; z(;vQzABBz+cEcsrKSNHnQ|btu?nSXV`6V552K;^_{;Dv^5>K~jHiZ2XKf(Au8flej5;Y1AN{5nt+OEMWV9jWZrU7XNr{)}vvkIse= zB_i%;qscw^UyWYj`$VQr^C8DCXB)AOW~Q>EvzdxnwlLHV6upylQs1*vV#f4aO+6KE zbN;0cW?gD9M5aI+cofk?#+*|4{;TJmY_`we2nL+9bU*RMd`A z>Bbotb8#4y`7m^YM#ArK28Of38!CK#wbyLNqu0V`RJTZ9b7TeSy*sZCj_1_B@rN3> z9-IC92Cv-J&TjsT;Xxha<3SwGczBj!{y3Cj=jSPJCP>{;mhHBHCA7EeEt=te{1X=U zDN>j`QCUO16$bPE+v-T}lJ~JXq_ErX_nlDDl08%5*Q<~k5{gne1oG9%CL(4Zc@qTP z1TI|DXxQH(g%M0oP$-q$uE z%*9>Fw(ruXpY-)Z=cIT@jS;(}N~HC0NPN*&y#|P;dje?@q0Us(GjzPkEFUrX% zOS0EqEdPkYAc9}XkzSdH-Z%upU|stvuMJK*3n0Ak_m@ES!+{ZMC_nrpnz=T zU_LLq{j|JdCP1ZVb>>?c!RwJcfl$zE$x}l-W9rI)mt>BIK1QS?*Nf{+wP(okNMMkPN=sAOY8o-MdmRUy6<`UUB%#QgdJG@-0V^ z&|*T^v%E3k#rp6i3v#ZMyv!g%{GzfSlu=&Lt4<~p#(<;C^>Q*yHDc5B<`d_9S<;ZH z%R9#x=#7A}QbofxXKmj1>50&s$&L0pLTlV_)iauU)^G}P`pwpr?>@#uXFO%D(g4@1 z)f5FD%*;fc=rurRoY$Q>_%4eE!WiZ2&xDRqLg}`t(^buvL!V?V+VOgKD&3p-@G%xD z`!MlvA&9f56KD{vC!Sjrsz>sc*Y4z7FUyr-mZgfXW?%QryQr7iz-PXd)DKYJq&H5U zJ%`lPjOkCleW%g%`slf_r#kK{&I*mLIweM3pM5+Wif1eGPeRg7 zZ6=!mSJ{>-S~ixV?px1`p{(pLU@-nd$RBdo{hTeA_$ zN?`jd;5k!Q53Aada52DA|B$o1*T}wB}u*-tRD+)M9k&}tk)6aj9(u>1#S0} zky16|rPD5JfRY6h>{rnE%)MiDH>Gt(ZyWuZ;JV+g=aH~3%+Q%}G=k)0x3YD`DuUlg znh^B{;FTPd!JEkng=z4}v>{4XK*4m3V`MdyVLq*9_pphKOz0zHaEpMxLq# znBMI@ZmGCW`kopCj~rOWTeE6zHD{|YF>6G1v5yzc{Y$8hj>Q7g;Xv>#bZ^|r2tN+> zwMEN!E`|Gf;w$54e#UI3XT#7F4pCIasvTiK<=&XEnt`Urr_mC2%+INBXie$VspV0D zu2E>5q8R)Jue{ZD{PyG5o8ij!cWyYi%6C{s?6qX~D(X<9Y1mj|+w-*ENYH4%7gKU6 z1KWRIe;@>e7;#_Gyhup(GL*SF2>6(vo;t{l6h0R#6QvxX z;Aw$NK?&nvpm|L6&(GC_du8-Ov&XmAJxyYbN2Av+-n#u`EzR(_jkV;<$LHJBU+1Pc zHzg?D9{u&%&qo$#|IK^jiAr-V!0=%U{aQ;6M1;r#R-O=|5oTw#5ery z%M+1ovF7`56xY)#DyRfKAF8)sk^KtZYcqPgo#-#^>x#pve3;Bbe# zm~$prH@yAlg_A_lv1|AzW{=}5wSBs#RkAhf?QB~H_J+XxWNzFwX9Q-Aprfd^PG&=* zH=;2j_YlpI+LW@%cI6Y<}*H-anAm3k0Smsy1rj;V8Z0fLi@#&rAO;7j)vw5x9tZSK0R zX*{*S?X@lKxAb3@9OT9!40gjG!0TjqZ5y+LLQ2LHi@KALGlN4vWx($E@Ga<1MZ<1P zXcV~h27qeI2Y+Blr>+J#X30w^;!1ryZe9rh9ZOTiRrxm6KX}7DA2ia04)gdaGkloK z>T$Z=BS74f4|^^4&X zINziiO;8<0tjM~NqByP-%0taokWU>vSc2@-65TARVe{aS+Q~VlBJNvu;#J+Zun)OU zWjC(n>;6CU`NC%>McM;K#@^%PLrH>0Za-Q-mcfBwWi^Wfc3#=9G*aK%ES?}PWOtU1 zUb*;&hv1j;9{8K{Q7^4ix2W?bw*iedJ=pWGq!?;0?_$`b!ZVDUhq%i-VoBPp9)|j* z^i+^PB4jToERB?DMkBOJgjlv?gq_ll9K@bE%rE{jk@|Hv0l)`3_`j_7&lR;>zY6Q= zPQ29i7*;eQ4d;CMIv~hyg);tYyHo?NAcJJh4`iaRzYl+)q{gOaZCE&U38KdK-s_3C zM+M(yA22+=sE?IKIhTsh$xA1o55nB#3E+=-@MnU**?i9>K^S8Xs^A?}br8Tyb>#Cu zf7ix`>Dy=zpaSuqrQX&|I1f1~*d1oo@#N23C54-mugUG-@|c|8LyXQWFL1x|4hSHT z^UgP8GE+m)!sc`oHGFHk+V8zo?Q?j_PAH;7In{cjEqAn$d5sh_HRD_0mZZK@y9!sC zuKBw+m3&z^R2-vKNcYaSnFnF=IG# zTxMqY#gw{2K|k@1g)^Q8VdY1eQG3ipRze>o@o`B#bVoHh@E0pT1Uf40+V0WAhCRLk zw$4@)xBB^RQdM**0WM6~5Y7 zwdnYon_zQk)Zbb&ShgYa3AR&#Xq0Ydw2f_O^9H_>i=I@@GzTU;ZB9@T{QEjxTKP@P zl8TOk{@0a&?4#!4WFVA-&C_Dp05Lz`I(NLa{b$vZhwb@F+r9n_l{)J6+7yctHXj26 zDjHf3cx_lpmFN5uJRF=^ej51~z%kQ4b>EA8xT`8fcD9kQU+TC6=~=xw)Ysxt;NP`m z7Vf?l(AGVm-GhwzuA}4d)?C^_gTja1B<@|TNf?J38cmaJENw0V#ByPaS=-RGjPnPt z*8tkP!`Oy?35NV7Xbo3(ajm&7z6_$&9slaU7hl6allCZpg|zuj>F-d&Lsy!e^#-j;NZ*Ojuexks}02%M@Vy)=sR{nZ1w9V%wm4q2lKNuXHT&3BrrxO3K>+!n;?dXbCR zu9l~Od}7aZdt8-^OGk0+=&311s9zrh_-@nOYbwMKy;Tb_!tN<8fZq?zNxBP~%Zz&i zud}ylfQS&c!{IEe6U45ze-2rN`asVFmd4pACa7LnvNsRDw}+qCSAi{4^%v_i)F<9G zD!Ij1A_sh>_d?U(9wOD4&tAuBuRe`SDE;i*>#Jv_l5wtLs;`~7^~*kC|GnS3oHOOw z%!zxrR3u6Z3H0Hi4!sP_G!4=_tAQkU5Xj`IFEpn>Rh576tg9&BPvTkQ8iZrMQAk z_XYzpItd?3+R;+sR%)4=F%2LZ{!ur+Ij?Q=^Zs2%E|0%}Ecre+b^ z$Y0Hm>-R1Y&46BJn8%yixxnQuv#u!wdP(E@)hQ_%=bL>uv;pvBZrQCS&-*P1h*k?A z*klPsPBPPm9g6!M;C1Hfn=UJ)WdE^-4P6{TZf2GD36T(YgrGCDhT4O^=% zVnBzZi$ZLe)2d`@uyJ{={P{L)>QUgB56=Z$0MA~|cuk+0D=_acO+ys+uZ9jbZ$mFx zxu%N9a_GNy`~I=;?VoVyEvE-23){!s308xqk$T)T0*Oti+wS1Uc8{YetT?8dx90`` zj9jGkIlY=VJn;<&h4~gglFq`nZO@qsxtw|5&J+DoAFd&#{qyzsO=pv8cK()Td0Dwu z!X9ifj^4CsHM1P{#fpG4KenE$^_93N`>0j&QRl_?%Rk<1n6#Ect4E{2Yb66j9uGIS zlqg2ysu9)6a<1{a9yjKO(TyPGULLoJHBPid@Lov28du}kP9!v8vbk zY-J2439b8VHdj9}1J`n;eC9v$d{h+74uX}$o7xll2NpOvF4XPb7 za%E$(#&kxY5aje5+6RZ(wJW76=@s*T9-T|=3bf?8fqC}84}0hJTDH5@@M_58C3trW zI@Rj_Ivk;s$b5}8S4ch^d4L)+!MW5DPmTfTSvT;`pHiYTp2g^w#^x*v1!t}e z#VG!-lq(mT9D>w5ygqm1*2hJ@F6@oCr1G(;DpQbeaaT?!$X~`#422dR)g&sT`d7n* z^;v$NERGHhK zD~BP@*HujPF0CCa7&ndHswG_P!gunK4aGI~hJK>TsJovTcV+xFtfEul=8h#Qe+`^@ zNPiEIlNkhfY-*#WWRRH>=1<0=v|} z^@o^Q#yInBH=NTjAw}kmQR09qSxT)$qFH0I(_KLh-}*W|9Oj6DDi~?6p93|p>u&HQ z0xQtmfDMb1_GGvGY5JjQK{t|$nwnXulmU1TMtFk~1WqJIUGbV_f%IwlLvE$(1AguU zJ8OYXS=V;7L}cg(tOB)z>P9nd`4L90bNCH8NXAmFipkKNC$5^{kk0W#`o?5&Ei#n5 zuergSaT%Pu)U)N{JZ#QUz7uUUsmxKx=O0CYH`V6a6xtn*b2#X54e{>{e&#XiN|jky9%1%ml1vYJ8PPn zResbjLY)IqJDWQ+4)^;}KE`7{oR*{;@mm9|7@zcMkq3X@#eoBr2MxThmhz8;3(-^0 z-}z$yoa*)C+Nn=55#PsMCF1N}c^_}X2!JmA+Za(49!(BT+;Vm56d!4e{1?~;xz~!* zFm>mAb3CXhu7Wp;HjMc%ciN92o!E48cf4qXnK$k^x5usc;rmj(_-fj3o&u8hs+vh;`Yc^_3j52Z` ze64KWAN7-aAXxv=^40p(JDhBT)5?~CI*(EdV7Om67dW|cYwtP{(s5SBkW-LS(gpN` zsZ?W%%8l=kwDiuTd@HQ|_r5m?QjMA_!Vsx5k`Qli_qpsskW})e1c`uv$x(HexWL7g zFLt0#zYY2qA6fYs_Y|p`i`7O&GdGG|I@IXtM$I_rBCncrl-}t2(CBWe6gUh&=}i*) z;Jv;FTG@qo+TMI57y#qEXezh-0$DBF`Uq)7v;=8?zFkQ)3jR29#P97UO{m!76;l5oaJb*?DKDvYA16gC}8ZG|tQt-oJ&FJaX z_}TTFA{YSWrFw3}+ z_{)dbTgX02@(Y$OEjvNYf9^uA|JyJv>Y-`F zd5nX{Zx6rJ1)eAWjKnEg$6XVMIk#g}^eWf~B>)mHjygI1aeQ~$usA2ZimoTeY{YVR zww$}NL~XO<#|!J^Z0Y?rOH5iKPw+RVKm3N2IYX-RRMg3yR`X32AW$}DTtH>OIb>Xf z+jfbVtkMzAXNenfv^{5=(-qCv=ee==1Ai-BLq$OVwG}q@>fymfTPtfJ_yDa9fAW50 z%lr0Dz8pc5h+7k<=_57O^$C@s{jR;EtQ+Odf8UF*U_)*6H%wt2<_GMY9ZF!^e_N_X z1XBZh#0&4PM+>6G^fH#MEfIKj6&EZrsj-v?gF4hlVYh#xZ`x$>O~W&GBXf%@*GNSw z=R(I{gYE_mNj9Eqb9z;BEwARP$9U+yRobfdd+oGzLpr;1Ra?>SF`}A~zv^)m=6?vA z)0L|7%N%J~9^1#Ns1t|W3!R9PH5%q}C#}3U?g&<2w25^6qq0_agP{)~ekHG;Qug+` zfh7*045Ot6F9D$j?wc*_yYm%OBJ%fv!Vgd4XP~#~qR)NIeGlIIZ4w73|oE$^;4w6gDs#yQR{8XrPua$!rET05nEO<<-njnogeF{nU zZ~4`ug!Bjc`h6*D-2Dil;k>|aipSz|N83H;@7Wow^ejB1!JkjKWbW4|KX?m!&eCB% z6F7F@!R3({#Z6`G{2BibayNqx7lk)w3e^z2ZZc-Q)K3`Ez6~Tg*|z@yAbn)3FE$-X z#$?H#OFt1}0*~mS*F1hfxrfGHzmkUxmRWC~$&wA$Y0tp;P!5@LhooNY!jaE)*XBu~ z-Q`-|EXf>$OQX~TElVVcTcImr&tvSyJaRf9kIMO8)iDSq^$g|YdG{N1Kk&q;3@$io z;u>9<6f22YE=#TN%4hk5C)8cPgt{*@2`dWz-*x|&M^VMM%BhJO4?Uf`>l!ygF2scC zji8r?^?Og=>8YrUEmh;wKRHg=3(xdwN}Zg(rH#=Y@;~124Q_F`5}I#ln@c$TOBb7v z@-QIRIB_lg{P;|mZRC zhKzbtObn1N{#oa7v0KRETrh@Hu{jcdnnjM9BH?9#dH%Q(R&hsdwy!Yo4TMQj8ms>5 zOHp=6`eavdQ_cCEzf_;EhjgBJ*^y$-zp(^V?%wims$e|l><)0br|eVEvCiO-OZkFe z8mZ@nfhEF_!sUL&JlOVPl=pIN{r-^hhnn^9dzMlrX$I-k`+C`r$n`a=DAM4Xk9(rs z#e}n&RruNQEBRKrb>HmW9X$+)uRruuk+Q7cZs2MNc&j)c>7@zbK_)lUHW#Xnlhto{ zCEF*}xk_c+2mhQvcI+&- z1wUZ}=eOuwg0^~;EF-9|(YBU%AY4iS<(biDr)7#_ww0ZMVH4g%RY<_#A7u4P>#jD)Ti3E-({UFP=`O)^s}1F&RoBiQ3JvYwCKm> zRUP>c_v8R3=WTKriMaY1REB+KcZVMH?AAF+(Lne}%(W~|jm7FkyLCzF%GVlIeuboP zxqpwPG6pYB-cpyrA#C^6yO8c-0tN+OhxQZ3CJx5n!G}i8C=j=kicwMD{csg5 zlttYZ&S@t!x|&EB@@(+HQ@^tickcfpxXNLP#B*arHexUpz%q?9-I}GSRjV|@9i~EG z-yJ(nIVsWaHJ;zLZP+IYpsCz9M|^QX!`YjARSC8ujo{-OA?)qH#fYyv&?>XG7w?=F zx~HL0;l~GT0v_k|T$r6hp`@y676v}^F<~56^2YrF*xjZjt8g;!hTBXH*T#lrSke1@ zhD9dt?SE7^*+11&6--C+lKZ81BtWKltJ*b(GBAf$-&wTyov_3St7g7>UY2gr67guO z#-IT&0=r(243B}%ENI^U+Eb96Q$&mxT0WnB5Hhhe(umZes$ZS13pHPBwq9A*NpiLk z4jA_9?7LIS;^10wc#U!iMNNHwhkd^N9Q{`;;Yp`i(M)9jh|Mk9(@clg8O_`q{7;?^ zE!T||)R}3mj*SokKIwZekNF14s9~&g^|m{?B<39!nPvzVm?kiKmuYh|1&Rm&0+a7M zpeA38OnJ4o>VIPGz^gSvzAu1oDmRWfe4QyUx}qSeM*hy|APwMZHGa=@`^#`xd-wYz z;Bww)XX!)x6~zQ}V>+qeMPg)fC^x#-@@C;zFcKkF8ge4ZnH)Bpm9QfWG}5%Fcp zUVag1_&UM;K*nQK@vHG8Znuo6VsXw}``fpa(eM}o2}q$@OI4&2vkrA4Fy)#4sWoIk z;reYxdgQANU>{^hA8OUl^v@)zUvwVrMj6!;o)4j!}9#*Uni5nCP{} zFljMs7l?jpHt4^J=3m%JdQMF9$rD!A2V6td&kVMT+#d&~MszM1&lh(+u;SAQogCxo zqRVE%<989}^-3In?jJ6kF6o9#Nw1T$@_EuJAi?%Wpe%Pe-Lu`Jp%iW)!+$8?r;JUY zKs_&61<{w=)XKe}#{+1lfF_$z0Pb}4Qxq?k=ww`F)J;{mR;<%N0V!(**MWF!1;c;f z_x4Po9?A&<4nIdBC(HL18 z!LT6nE=tK?2IdZ?pc>|2pSpdS)PQ2JH>k3ZXxX48b&#+NQ#MsKxzas2u=b!o9M*js zcf(P&=lT-$cVTU+v+es7L;{^GjbRO4}B> znMr1_d@tixWYpu6#U$K2i7QYCJ?lgXed#a%brc6Q5^67GNn|YCKT$2FfAjjYx97%H z{vGMg%#!kh8sx=mZr@>P^Spz4FB2S7Lxclv$ikJuUOFrE+0cx1jaQ!J=l%_oh1WQB zRmV+BU(=B_nzMAQI2-2Ya3z^wJXgmHulh7&CQ{+jungUQxX%+*6=LGOWCNM-^0X}c z!Fo6>UsOt*UZL9UkT(O(Qav6P{dAWq50m;Z`r$iNO+c&Dq`-=L+h#4emx#B~`qDY=VLxoUY@K?-QiVOf4DBxcKi;`TkHxDz^PV4`+N;cNaGa5Ix;F&pxQkdRDm&>YC%2YjvdGFKflq4YD{Qw2? z0cYv#AWgvoA6L&eEauea&|K8LqtKLqXhs?Kw54?(B~EPiHvVVi|HSn~&I}3SWw@yj zONsB(aqN&-(<*&zR4>@^8pmTHUg~+9Puo`T5=Ho=Gc9e`oQs z?~a@Fo$I`pKGgD76rI-VRG#pW(~B^R1_E9N)6sN&y-O6J*BHjH9Nd7IxCeVp`{PYD8` z{sUZ2j0y}}NtMCD1$S_~2C(Q@?=0PS#^-UKF-pZQAhvo=y9R$t@Y2}Ev9qg7_XIvJ zu*1UBDGP~?RsIJ@o6O1j$5+i8Kd7c(jj%Ny>nk^;QYGU(qNV@R(yC9yfyiON(M3Dp zc`o;Rh4vKL1D+Iv`E}ha!=x1_qGJvekt(1_LVxesh?Xr|y1GowS%o@v_PsW^o=j+d zy5j7Gs8gqmNgsxPSKi)ow7<=p<9&j-uBpe9vM#)()Uu%5GuM-TAy`>Chelo zRQ4Z<%xuo~4)$l8t(uO2@#l;KuheY5=*pGeEvh(JsxuotnXIU{%ee=2&op)EP1#E~ zzuAsmeIMYoRi@AwputqJyq7Qr1$tI0`M?~Is+D;0oPC4WTvqyFcYvUiss{FO)hz=k z*YK?>RN2n!^r73i;?z+z=!&(|^?`x6Rl-ei0d9qL(e55tvsF<>N9{=``!CAfMahE} z#9*P@Z^tR-T(T?u{l2&xOx6*~e|N5(LL@?F)a~S5W)=W3oO^FVxEnH>4LO%ytDYh~ zGSUg6zsW0Lyf9jOg{+{2Fq6}ZJ>y-b1iByIY0 z$wrmvyWw-7)u6L*>uS&>vO{RDx{GvQl@Yil-NH@{Xgqe@dg8#ovzBJG zp!h2tL0`Z4Yjo4s=SaJv+)ZZyh^j?%abJw*sl&KGd6K_RLQmadMa#irEnV-e?S|L} z%1uV~y{2#lkob#JrTE$e4r$0UUYog^6q~trW`GnbUbR|oHzce0{ZdD%k%D3sKKSkE5T!gtK4(Ce&(*Gawa=bfF118X!-ahKJ#Kbz$CHV=YfYn{Zm_9 zJGS#e=CXTl|9|j;sxxhI*P2ltswR66x$jkrC_KaMz`xBFJ-=^7k}l|*6Yli6-Lh#f z#9GJ!wqmrpdKzY^oh0az?!hQoU0Kz(0QXUjAOZBVf_+ZKqEL!NTp#(PN3kl(R>*ER zrz=(sJQCJ%O{B~twuY8j&f(IRk?wg36%JS0wNtr& zx(pBA3rkhBYhy>M#VgL|#Bq6j$l`4-az0$7*XX|FmxyHM0F3(oz5Xp-gahhwg%7O; z3j10RSwre$0htvjjn%-nQKQLt&?Cc=@X>h%=~Xu^FN6lk;mR*|6I!wCLAm>zA+o`$go$7~5>ol~Q;LJ)h(-Nv^-sOEn;9?xU`F=k{KoO zA(hPQTRXQ^!Hk?!X=qF)Xe(fn?Dv#{-80L(iilSKXy&#dhy4RsMa?CPiI*G)H}ML= z4c?wivk+Np*OR8mY^7sEQDy6TGPY|e=}^gPWV#fhA0%KcEg#5Xuo@gT^CuCdYL^_h z^D=iWG-vWF5YbN+wBI#ma7R2+Hj{gFs;Bw^$7s92poVtx^pM`KprlW9_?_BMGLjm~ zfm-S2?x#pzhl@5&|H9&*O_Nczr3;^Py2@qu3(in>Gx9565^REG!$mK47q11FL9bUm z)JjdWV%=n3)u~U|gUH&&r*3%9-7$9G7mJ^fmaTi)i};xtW#x57Jsdy-sx zsZFZJ(Z=D|N7Kp!H*QHd#a(iWIhD@lvh-DtHCv9Lay>l(6s*>ed$)RsTkdpBYD9E+ zy3l+&Xu#ock)0@nu?jE==Nhy>0eU47lCF?)wp-93!iH~07*wWLv?^TiLuiysihr6f z_6h?2Y*j_Qa>&#ov``U>__zZ5^Z~LESAlE5{`Y8s2kfs8Vy|e}eS1A`n0e3rH99`* z3t|CrV+A>TC*Mp%s(X|;lCt$eDI1iTk{j6?aV_XJCyg`ropCS|8)xVL-*=N90 zp2!(SUzlI-9F)*Sqi(xfYgLs( z3J^{hsT;?u5v;q)HFyo3v?#t5vDxZmBUtrB*WDzbu24xXZRpZpSR=CBy zWC=X#`&+(YXspf5Kb~91PeVgpTe&Qc&SJ-Z{`QQuls%bETAaLfOxkl;F`K~h!tkwr zZvyi_bv9p8ehE^*5~?jEuRDf&%zD=E_fcguP^#9O?B&_(C9l}r0{E8rpL``To&>Ee zAN_j-(x2};5h@qfqD2w(7#a1fQyZcCyn2pbZJj8y${9CV8CmM9r`Bvpi zw{yB^?{GyUC}{#-J8ffc0Z>`J8&|WS0^C(6YS@Ava*|oqGctyM^RT!Pe5Jj(J+uj> z;juK~)OLv-8F5lcAi-77EPp;sg@~#|z-4U~N?P}+s zEI?E~<~y?hoots~4Z@U*%$N3m@|fGVMbE&Whwg5#bbUz)K_s}*Cqs?$*S&ajtDX*{ zG~yH@PbK2ImYr@8G1}h6O)i;H^1}S&-H+?s_Ig8#2i!PLnYI4d<)%!#Z$=czC$g{D zjAcX{6_VcTr0uN9Jmj0D1M8Ck)bE5w)x{yj_SGf%g3Bn|p5Gb;UI1rPb|P<$~yj$8=0GE2t0-;6QnF`67XcYazPLVZEUjd_bNO^An0jqlw&Ki)2V zWcUBIH!cdrNcVe)y;>$ebvpfvqZH4c+dL7_=zk|GY3kMFN$9(PxsV(gX)7oI@$B&? z1_QlQmA2&VnhY^)+)$#{a+wOsmNtrZsjouXTE{JwgH9H&FN4I0RUD-&jyhnLWc7%q zhjD&aq~1G%?ZQU9;W67qr@rt0;%;W9`iK5YfW}v7y2wXH{>bMIUwH9TXBN_dZ2P7D z$<+6IW||mauJa{*a);cNl&Z66fsi&XdN>+rWPH(hjR8f8y;`cU4&SQr^eCg$OXC$C^vKhKM5ahoK#Y^9!)CfN>Z5=>*@3t1G;9S@6e`C)lEG&kbQRXc}L{4am?%_M2NV{KX?p7}q5Pw=%yG>CY={-=}6-KiM? zfQiceuj7Uq`%s|7h#WK5;uFbNQNxG^FfXNR7`N=N zmnd0Zwdhx$F|Fp@zdP!|u_-cZkSnO@C1z>9qO@TZIq7%(tSv+m>}+`*1Xi#!57g>@ zDXENowWf0J$Bke>uB9f+YFxH^ka)wep;Sh6##}j}-fbGQ$?~&aw%#KzJZ`vUnDEh{ zle=Y5gv*;-IpDHx;6pv7jQQUa_I`N39!I{vOZKEG(O|j3?tZW6^QWC$XC&oOlKzwU z`OweT8!Du#OLw+z!=KmLy!tg!f)*gGK9E%&U8s{8MR56!={q%)j>vaLZO|85*^9pllR6Zf2;A&}NtOy`W6 z_4y}Zf-AV=Pip8*i@g#bx#=#;$6db-KiWNiOHbGkwT%$MH9 zv)%?pD~@tMf*)vW!c1z-{Il%sY4|(p6vAl$CMBOW&FRF#2yQDKGnKWAKA&{Qu(=5G z*%x%I(d_Rs{SPjmZ=@t?Qz4IQY78FRs+28MSNGAtCer(c2$iVVyZmP zC5A;3>W4YJSOJ==5mJ-@cg~#Z$qzbdE{%7==Ug3uL*l@m9EY0z4k=|`4t=4o&8`u6 zc#B3#I*BRj0OO4v$9BERFBD#)Z_LXI%qW!fbzX?zc`7!iTExF){=G3z5$2%+2NAAi zfPRKg{`q^8pSd>F$h-HO83o%rzyvu%)pRpeyzjKF$2MQr`=rH6Ee+goetCXlw5(#neKi!;OR{Z6-r_fJ=307{%)gsjc5 z2PmJxJ|-T5^SXJ=#YkUn`4|0JD5H^A_FVg0t49dAwpvhCr<%8Mx~`mY^;}YPgwDH{ z714^FTM2^OKZuPtIBgCHnAaR=&b4vF0HsJXWOJ;>?-jE?aVzoB4fsGRUNk<*B?5z_ zaLc;P!egPbLPYUf$Kf=AdI_ht`;wJIS#5vjH6|%B@>73I7R|2|9?R6tL&{~ipk=%V zjtL!gXqNh9;y#?68f%Fc4Cwxr`MY83Z=Qc*FcK@R)2s2?&$#HafMp0LRrwTg@L{pD z8i9A^<@mse=OI~U*)7TH=W+Z#fR2vr09dqmrPk9JI}UYiuSnXAB*iEOu*^cjVige# ztL>qz=EQnEaOa-|o5_b_Xp8|L(TbPItB*$;ltCAUqulSUD#j5aB`#J?5bolKrDYyR zuS}OmWjQ7x_5pL7|Gieg>k643P=^;ogQ)qgj7yFf%bl)24^ff)*}MGZ7$JjKzo+IG z_-iC{rlqPcF@CbXLY{O|xCh=_$nFfKEv@IzmP*4i;Y^S^2^$?Ol>4Q6p}3oslZ&^+ z5*0#EiMxiOsfnENnwXXk=cI;(e>n9}yW5SS4Y+KmE&U9B^5&uQ$FL&my88~fqA3=C z+CBgVA1b>W;J{kCq>5c#Ul62psf`qRd^BTVy%Om?i?+Tnt=4)#NV8lGCh$>gMCOpF z9$4;raZQZZToYnJfeC*=)UWo7k3gL_S*p+W+!>A9olKwHcqjaSSLri$xXi>`Kv&?q zvy{}VIr61EBuqk9%h4=`+I{MO{y^soKo=$kLx{IkiPO*Om=JI0)T;0Hq4+ysu|$mwXuA1TTPG9j)I> z7x_?@zC6lo?OO9$;#JGqAyNhhL}HaU>H+RugVzT!ki6Shkuf}qK_r3nt%@gDJ!bz937%t4dK=2$_kOS&TFP6=-XQxwB;8fOq~Nel^7MhA07*_u zP1p-UEYZezM488lMsa?9o>G!91+X~ct=(_H9o6auE zQ_{WVJ4w2ZOU@0mKWP)pZJ$s0Q&$q9;Vc_Z4C*8<$IX6BI$+Ny{O=k(V9c~y3hlf% zf7;*)j_WBt5)f+F>ByQXi+kp#p&{y^pnl`Bn&8K9xg9j8N27Ss(;OInTxep^EO1&K zPz~!-nAc8&M;V?RXk9rE{vF^OYsLuFDyk`Nu&w~jo#bN`*52ht4{9&_d|XwLZ7wlH z8Hd>ZlXPfoCEOeArhL6f;~S%%82IhDtG&Lki6~w%Y0vOK3a7r7_!)oUSy;7KI){J2 zovYB)u$owq(QVO*uDNCFB}Q26?&!zKxrt&$k52B6Jp^EvvgK*Q_pEe2l;KQY2;KHk z#u0kH0BYGG5^lB{!A@4-_hrJH*6udgM6fMy`^|=mKJ6T7wPW7s!D@ae{ik<#lxT+% zuEj19qOzbiTUW;J{wv4!qxeyBiXqPXI)EuSi@SWxkF*vzQefi+b~-eFZMJT(^<8D7 z>g_M%S$JfXpNNj1%y?GXj@0t+u-S$DfZw$WYkNSR!SZ(P*uqgzKnR~5{g2}>ip=KQ z*WhfHzn;Z1D<4uWX91BO7=8O~H8{!QURVuiGtf%2EZ44;Bcb30>!_6mzZ#2&x_?;K zh3O)+`+RKs{mo7cipPU@u%8(agaRn}lltR_cbJ1*`{22wTEV)qFQJeJ7TbR6oC-y32Em;)zqm^Y_B2R1RzNzoFEC&NuL`_(b=ad4*d# zzG#1w@u1=Q$C`Bwu}$_`@1{fVqR`EI_4V2>2)aqY)X~hf(O{Q+*P^per!!;zHIBst zrzN1*&0PDx8*suM(zTdODb97VW^-*y4lmRdf_-c?Qjumw%mr_+dGBbN+v1CMS)(F~md`)9<&WVa}$-WxWVlZ%;BJa1mG5kKy;1D+82K0+p z8_i!!gKxZ#CAMLexQR5<&n(axC!z-;R)b{qJ(fhQ(5FGZ&J4E3e5Fckro0{GX5mbx zefn&zYHE$ALQ^@aPTXOG1eWubT3x%A46@yZ!7pQP7I=HhP9{NS-NtbY%Qq#sWNeXRTRhOT>OrR+S3=DheqGX!J_} z9F`bG7`tQ*3e$eEB9lr~C~1icGGf_@v|7@d-|}oOib=U%l7UFbQa{W_l(5;(axZSx z?PvDg31jmkYJTb1kAYfug8$Qe1GIW;n~pn$`xVw3KRU{_Kl;2r(|6;sc!+gem4mP- zmAH9*$Lp2+xv_u>>*E2&jY26nVmTnfF7NWTDYO(2TuA}p8&zV1j5=%*33r8%?LoC+ z4F%EPs>UY`r6RP{nPHu$;pE183b*@aY<-mDz#hL>NVnrvqqIK_j{iwJa-A{gv6_vY zUM}xbM5gdS5-;nf7JxD7FKW&V7`L_E3~#04js{90Ri2Qo#Pu@zZ^HQOO^h47U)~2F zA04na<;*?pA2diA)qwr`ADZSU+`Qv}{I2Fx^v|ixu&?(dqSbM$6P_AcSmuYIg!Qf# z4N#+1{AIy$-+5So-%DFbC%H}ZTz9#Z;_h?uq4_Abk`@u0{ z;^ysBZ9S}P@8V8^S}^wCN#&PA!FXyxB$3-O$sFDG5fr@^r#M(64RZm8z$FUPLorS` zfG60Ma|{+0mQeYl6wAt}au_2$W>#FAJq!rJMPfEeUM|cIU883*JIRr1e>e7lqh|ii zN0mRWNU?2jvY9?Z>*`udnNY9#D^SoEiNTPRt0Bi?G@16a!eVIg@~3X%7`Q1sUFJ|LxRfqreDy$sFC#R7d%LweN_Ly3PC zGA1XAgkeb)vth2wn-`j;>W}m`IQ_^LwfV0OmVlLih)osq;}4qFRvdCb?0a#xMh_f2 z+;+#{)8_T^EzZKA@LT5<7fYLHie?=Z#S@57#$0Q}4L*}W8$M($mR*O$W-4g3Eq4#1 z13s+=875>H01iL~WON8D*=GSvRcaZmXyi3uUJ8>=#8MO#2srCd9(#JsYXolHuASQa zc9u|}+D5apn`9Ce(A<61jI#e=!W=uQ4VHOK*W$B$1-gpXr5u5Fji zRSK{Eq4#t5x4HlGF~ptPjkx|_@!#=_5#gD-L7n2YaMj*EHvjv#R(t1F+Qz~ThGXTi z2X&s63FR+-sbT5!_$-rGG+ZPr>{h%_!o86^p+Y9esQpFwmAfMKgx;fe&>>!tZsw0~g}b)#nDb0=QT;X{$iZdViJd`6GY+o@ zmeonaC+m#t!Ozdofa%po4=)}j5+$_=1d*GGY*9TxgiGAKAs;sGsEQGs_sukHT%)_W z$rF9cfVkGPW<(0J8`LP57k3-~6W>qL%28tv-T|u&Kf#vedK=L8Z5mL=BYmLtq^zFanQwt zuc=|HCiOmAId$-pM;a_FMJo5Qd~kvWID2BNMSqsF>T|#%th>mgO=;&AFxH}m|BmEh z3~Ng*rNxH@uPpp$?cl$1XC|8PZJ)r^=gMMvC@jl07pu|vy75G=&O&my>-`u&0XqS} zo*k(&i&(mv(1xqS&HK?-@0#g8%`O-TL)8>Z6Gnm1ZX@V7l+vPx_uYc3DwWLql#IS? z5UTPBtrPQey>txJ>=w0~9>Z0tYf8@Mf3z8!jAXSpFWEl(RboOBa*37>Id17qtR{Lx zWVMquT-I&cX~K{oyFbh8hsPauH)IBT%QP9Zk#t%9uf+JNdV$YzvfTt5bCZsJ*ZmA8|w$9R+9T09WruG^a z9iv_eFazvOL*28QuDQiY4bZO9Q|Zq}-Wx3m-xO^Vx~wiVx%>AW$~B)iDCPFK_Z#rM zA-#E6=E1Od60eiRUfNWNFH(kD=zI+kOdT}iXxr$KNetJE!<1@ULsjmlTpl9tE_ zQ=W=zy0iVAj99SaiLU$}9aD7G-C_IDNj?D0uhvB=uL%4$Vw2hjE?y0tJF`??;6K(J z-D&Z}kh*2wtAw{2eb}YVF(_hKIrxu+!P~5zdJ2b^_}aV4E0^}JG|B%}od%G9MaeMF z8w0=G_aGazorfimFXOne_iSObI!N@oXq2N)G;iX2;iyAHmzX4myLa*=QNgha5&*{=oPap_4-hlPoXp2g-)k=sy0)-8qaj!L}f$tynRMAW-h&y|Rm}LgJq*mogN;qfoTqy#!l(ZN z5e|r}-7#5i4f>J4o@-(osD9na4<# zaDd}gx_#uVXNFsfuDAI+|L25ToskKrko0FOT=63i?%_UN#N?_rTy=d{zqi40Xfqa@ zixB9Nh@b8=`_<=v56}M@aDsruf^||a)Et=Ly|(jYQ=%n6lceE$Y3$Uu*WZBiHt{$0 zJVj*FMedX4;16Oe>(gL!X&(RJT|<}Bk@sISsVy&nyUUy|u^;nbmrOn;)UyO-PV&i- z`8{+{Z0=q*e^qVl71bZV90BF-|GjQs$DzkGuUPiH8Rg;#k>3cuJV1LNuaA{gE%7+3 z*+2cVyy!ykW5pIK8^TAv%NHr}%Ja+4)L7ir_2Il#27^~Q(cIVEK z0<4nm?SdPN)p_Ov$8);wdfi)7ygC_T(hu1xW6bd=6o8=m#6XFI`?yTP z;ygZ>Jp%eova{0=7<%tWwHu3=SiAnt(AG-MA0})(Ap4VLGnOiMZQwNH8Iol8ZDjVd zdDbBEIo>uSad~LlH@1{7QAoSc<~Bygj@uZg9`9WB8?8X2aeNZAcq9|9q5QZiR&nsX zrG$&!H}r7))sO4xg7jma>6`juhT}<(nYmResNpJ`aehih;>NH|=h;5b?r!IW&EFHP zi5idGX%bX7FE3|iWYV4nI5bjyHZRZ*!SsO*{k9?cj&{R|)cZOq$}9ix6A1rxwsb7p z$PB?@9rbF2YGuG|ls!jQ+8qj2oI|d4d_I#!piTH zHOn%2($ZbF){cV%$5SXZhlOuEk*GGl|By)#Op78 zOp}z|4~Yrun(zq3ufl}4#|V|Z=Id6z8h|Vh$*m-2bn?i80hiow(Z(x(qcfbt3q2MK zfA%UEmJzk!(`Y?`Cy8$((%+W-D*#{mt2AjN&8Ks41Lx}Y1Mtz|E^%Si zXTj%}`BsDQ<*5o<#eH}0^Ee-;Q}%7Q0ll8Unl*oNe~T!)F<9`h7(a8g zCNg70*$XZ6iYM>Ywwh~0sUsk`d#86JAXa_c#|nXvGPk$fE$woMsRTZ*^;&+5Xk=~Y zlXBQpM~ zdk54){b0_L#BJ!Ed;anMT7-I1wybgufj`~#SHKqm59%kEC6UUyIRv+LZu@NtvN_^0 zzp`rEyH;%y(Egj->w?zj*7H;97FA9KES;g$H4L51DFC3QZ5GpnHOs~u9<`IbH@4Z6 z-YiRgLamuV-`Fh_*)!^2gAOnHp`qMBzi?L%e%6Rj)kw_B#8gk1C13<%eYFY@%!OWg zw9h<9ZGQ4pJLTA3rutsW-=#jFl$Bsf zjG6SFoO+dS^-^22J)pOjkA@e)e+!+<0eM-6Pz8Wz+zhx^p-gfU=O=Ql@U0XMIc!+( z(TY!MEuUFdKSTtYnuZ*M8PADXasflLVRhr#E{<*&{|ilV5riHyMNF&AFV#jwaa8F} z;6?h=LlA=ZM6ooMj7Ff(A3Ll{RN zj3NHnRaY`meczaY*=a@1dwb1YG&H9V4!xHefaq3ngaN98F(Op`bj6xgLJX8oS};%Y z>dqg)tBq#2iZ9=s`oTWM_zfWSwj$r-#1PavtO+RZvsf= z*7a(Iuo1bnq0%?$)2Fh>9{*CK;(2Hz#9jls@n(hn@A0!#Kjdl_tf)$SGk~CY$G8!} z({RdjgY!{kyH{1Uz1RG4W7-J%z2Sx=$K7}s!fB3E7eJ~Jfe(30@kV4Qe-U-{!9sY&;fz(GznO@ZfpXoFO{Zo%zPqOU4KAV_bjZY7$f_ z1uky!KE9J~(ZPlS#I_pRApC@i! zrbbl5Wv7!ld~32;M}1*M5&KEdu|%80o1;1@fcJW%-rqN1m)JT*;dxzP?%;yOGb|wY z2%i0lmuXabgu_^7uWm4PuLB3j|aq1hU)uKJq?=B#w0GS8*2R zPD!OOSHCY*jC9T-Vs0VyojIoE)~p=H3foyH;3ft$R=|{2%Wqw}W+YO7c7A@@iEF^N z^e`fG0P8j5<^vcgsTsml5IJsEv!=j zt%|35OZt z6LwVn5%^sTe3oXaW~{aQIMvq!zK^DcS{%kvxk|A-FoSCwgFXGZw22Bed~tn=Mt677 zDuY~y?Dz*$gP2)r0JwdzHvz4}&GR2=3=c#;4K(iyi1U1edoS_*@=ET1-waMFc2G^Q zn&^A<^il%OdB>m`L>?$v8U?LS&NEAuJRSj6-NJeH$?{wZ9iat*I zRWxSsp4K57)W7jg9f%~`d<*|ZV!j1TCpk>!QT-M3O0!j&2LM>!K%gu6Ow^)>QERSc zvtl`v9}8H56l=t(l039Y8x3t!I@9UH<1zTx4g`cW5?M7iIdLuztGy~f;(+*HE- z13eGHC*b|8yhg8E_-idq1=nvLdRJ&44mgKDQxPrpo|+eI)O;A4{cYy?#HNNorBL;M zJgTct7{4ws0KnVVFdnu{2&md!E9xEr=;E=gNw)Mz{{ghG3N*W^862`$jSfmSI{9Q!*VAY zM#Ik8Wl>;xKSI}h*3_r;NqG2#mn$!J2>{5$=|dzdgocZ)L5aL-h^NG+rn|p0r4xC= zf{%3Wf$>JfeImp>j5y|CNVi9Q4>w9{?%cdGU|(7T1dMcl!t8) z0@@hxl53ltA7b3)p|<~6gbZYiz456qICWou0hmC{%;jgPc{noaMB*Gey6u5BY!sJ& zfa&4PRO@~#bV(pw&d-Hr@N(GGI)uG5ACU0Qgep*T#jsShjO*%T)mvaksiFEneUFMg z7QCte$4~6;Xhj68n4dfxSG;}rr)^2;z62rA8*4~Y3-NMY`Zhj;2Y)0g0H?>O1136WP*(vZD5PwPq;wU{0uiA`R?B~n zycH*iM%NlLe;Jxv-It;bOsN&9kmw$W5r@u$iMXPi)z&>vZ-|^vYJ6vRW^_%uB*A&V z$O$}R+AoC#_;2>1l8hBDT}%}fuxL2gFnYtN=mr3igsPKEWlu#Dw8hW=s#{z`g-1_s zbUBT#Z2*m|=LV-^KIx95EBmkbPD-vTKYDP4WgF^gk#UJh$8LX|OoMIwR_4l#bqE?M;l z1W@0O;{cmOsBs#4XKkJ!h7K9l7oC`OQ4FyLd-;4Ef$msZS;J+J>QH2kMAS&$svKK` zTvzr-5Zd#`DkLRiY*gC4OtUfvnTM4yiN!NeWvrsaH{$gC1|3SHRjMd7yP z>AAUjB0}M+c_u5&ATn|7txI>Uioi)2u>n+8(0I&C zT;tL~%7>R#sfYR=(HZKop~6jqe0E44O3(eEYQm(QN5sLg5j#-yMvGTxXca($0lW(V zoH6gTWW#3}+jgASYmSmJXIZm>7J{@VS%XHsurR&;VC@5{<>$!x@Ny$=IKs~;-b&MS zzOC#sxh$FEn&V@wpE`mgtqf`YhK_ESj2Q8oGtlW9Pt+_{VgQAo{0IMUkSSGAO-D7il91Lg)}$2ryJ>p(F@|s!~*nQWAQR5+I>TNgxowMhi`P2}OYf zsgXoL`uEPd*ZD3#S&Oi8);as^^X#WcPW@OYm&cHKI=wIRHRxT*2{UaC`n9?IF>OHB z_{?FXJah@VJK@&#`4xLfU{womv+`04-l(dt!A`^8v%iP1BNu(h^FHq|wajHV$OXFz8|z5f7pK~P z*>wHPZ98zB>N9^eSfuVoltz9WYt?l{Rrbo>H~pFC!)cBKf@0E2Nh#m@y0`9(%vdod zvqz0hLEmSj!VM2U1985?aE*}cMFl{zkhzn4e8;2ug+K{FrIvlLvlWr)KMHsKl}8}3 z+4^{7;mI*05TiNag>lu#l;kDM+HTkr5yXwB>uqI=1~$5gv9<@nYc|*F9u95a-EXn5 zABcEV^2W30Uzv<w>J6=`Rol`kR7f2A2lL(TBJ6_OF?_Wpe1> zarg9<_h9wVSf5UnEyn^|I*+*{XBEUVv<)c8=WLta>?{sP4S~liVAd#Wi%KqJlb%u4 zRa-d>WKUd-J$wUs@duy{>VuM}N~i-MIWaGT_^Aj(WvtWt$z#1=F>~coVhh`(8P$ea z>>6aLVhBLHwOltgS=13~-yI3{uY(+7A(RKyz7EkoV|v%vzBl_QAs}pV#;4#oobCTl zYI{uRRcjJgI4QiwZ+zlzdu>Sr;ku$9Ga;0_^f7WJW}7{5%Jtk|wm?|iV%*WIm8PmO zw?i)cn3)HLjBA#AES+f#QRRvbDQo@{h&OEyx-8EnXh8x#0hpGI&h99h+0MBynb3YR zc8qY%SXQx-srRE(Kf7&lU9FFHUlH*@B(6TTaL&~vxwIUI=gm#Wt2&6M8*La5?5t)o zak(reJyZ)ykwt4FxilZPvOy!wT*BD_$Z)+C@_`!OR9^83vZWGal54nCZnZ!uC+O`* zaZT5?yb;|G7gmtXtNr*|BmG=0X!kP<<`Z;f#>ioM!-B4urEd`{pphbTc-pL6(eXz6 zxOESK=s?vRgS_`_uvIQ5+?}~zGCeQ+x8+9Xgm0}i*SZklZn5~ADW(i{{+_}M2w!SH zSEAarZ~c(6YUnTcDRJTtDllG*wE&1aT8vBI)sE;t?SE}6V#wUR-JzIBtiaAaM;%Ni z1U>d;-R4%=veUu$c2rc4l)JIthDDR~(=)Oaai_qs0=}fQZHPdSLfvZD7<1U|@>Il0 zs7CssJ9CD{SJWC{>ZlPdAorRmoe+fo4>CrNTr;jzJ;%6ef&^l?o|k9pBfmRJZOP>; ziVN7t>f9PVfA1WZ!k>ZmT<9Ia-ik)ZeuCJ#;wK!ySLp|$ciix4 z`S-6sQ%6ky)WdxYH^5F1(03nuJmh*_*-Ef*!xzzO@)7XqVPBNnlz|pkz_-6Y;qJY+ z|3-0Y5Vo-1Dg1d*FUB%zOO?#*`Q@FhqM5SelCXOI|7_BC{KDfVE{U=aJet8XqrDgj z2|>$^h8Q9W>^D3ONlx`zoPlqx3Y)(kXZrBHLXI^(L-jULA3vY+a5Ig+f1@_Stg^DN zdlccni{%T%Pc?>MRl9S{uxXTergD zjR_8K9Jw}Ay{{m??pAI#<%Jsb2FLY3%dkD$^SEjRB_Pu&qPG9-dB`R@ z%Z5JeP&D+K@?o>3?R3_~DBKpXvmJ|tT?b<9o6%tp8KSBH9+ag68c0cY>{+FY4rs@(DyN$OX+y73Z1`u*H34 zkSO^{KKqUo5W|PBAuPKjOT9L1=KWsxAtRzTaZGxtb!foMZOcrPnDtlBIpioMq~#x{ zW8b1I0}o~5{9hvx7x8K75)in5dmH*tRaj8iYVikS72x{xy^DKEcZOeLdNDVTYD@D2 ziUv_|8~96geSQgC6_6G}0n|%J$YLUCUZIVgddF9xv@kTw@_CoJWtn~DD3DN? zvEV1lvh#lDA89JVt<7#^_QR;%M?|=T-TM0OHG0U#^|V+|ETcn4oUSCv^;NAk=hh6Y;9HxkxVldKCR6-KFP2rnpf&sC-@8!Y`FW z)3iNrX#lK&?J5CPUCmw^b0dq*<ovkX-4tHn?{mKk zl$<%}O!OB1ls@OYeQ>$c^Ye8!`)ef!W?3C9b+z}#Jikemn2n^pAePpdj~QlWa^pt$YMDBBeMpsF+Xw;%y{Zh+zI!j=$s|2MJ^JR57*#s{G*!i&!8 zPHimMuR?L&-dZc-R!FTT4ME_4U-_xqUQnD16b%)_?n7c_IbKk#(k&L>n>g|Smd$`X z^DI~&(3DGJ)s(R`MAe(PcxB(!osYwTqa3`=jIoUVGa1nWY?#e%!py^Lm+bRaZoc2; z($0WBQGE)zH&aA5{Ch{o7Lms~i8tKiewMske3A0^So{}$r3q6^(;**sc zn!Vq5?1tFmUl2*|$5MSs2#DivFsvKD`Vyy?38V9>YvMv ztwK7Q`06WzXA)IZA?Cm42b>iIm%R7)DWmPm8^xcVytSj8nmaBSi(p4Tfzv%Dr|11Y z&s)k>+*WmcK1(qx67`V9zySLRXV@2c`OJp8UW|)7Bh(h~krOce$rEG7J%Hj?hrhcb-}+YxrR?}% zr^LG%b4;sN`}`d2r8O?k?wZEta1>=6D}F4tL;3Yrlwoj9F!>YOj{c}Xbd+()wN=TW zhleLdG&Omeo;Vo{BLPc*G$8-xa5<+ zlYi*-rt6kil^e*D=LrjkYrdh?K{jX8WxocU=B};Wgf=YPys}a8yK#`x8uoLwkrIhk zHP|hQ7U((R|L#cc$><<(u~^CF#ZfYbwqH8-H&i%O5=V@~VCqJ<{xFQ2o~eU9>lo5Z zL7;=+u+$FIZpF*bb8uuV2X9~HNDTsQEmo`KW$V}o_sFRI%xXKBx&qZw7ki4KcJsR~@kThW^!Gd#Q%@MuuW-wx7~8@YR^I zic4*Nr^ALKxdF}4BcR58^5O4W>1We0o_lZah;MejX_dW;#$_~%x1kfN`$7HNEe%_H z%~U5v5~k5 zuKc%ij%;KNVb^OD7*O)jKW^Q5?46yj+jF!mb2Y94%P^Xhsi2N{dafHDP*w&u7Z-t? zyUI+Kk|yzl;`#_ks}DGIQs#3_*INAgjws4Dq&lynPXJ{d5vV(ZSv^emE!uJRVXj)p ziQt}&^t!axE+!#*%lv3=I^ekFR;Yzdn}eCoq+(`isMTKbHH$rURu(H#Rq`|?`u)0@ zR#ask2%k2`hPh^e`%b&o;pP0C9eSj?bx(aDqf;ZQ0sM0YY(LK!pSqb;=-V66?^rfK@i1F%??4?mX2eIx2K z)WnK+4_CMl%*hQ}=q2n@eWG}UjE2=BjS&IC3iP0D#Bg`+@$*n=v2VpTUH5wsud91& zkc~wmoQDXW20!NdSJc7yPa zI)A3f#gu z)yo+_mDu&tbj%L$FHB%s@6Z#j6{=9A>C@E_2zF=|jK1gMzFwzzfy>VDgCMsy000jk z0ctBlZ?;K7>5#O!-_p2fz4D)dxw9UmnD$_?L*C~%$)srUQW}{5E~fC|;Dugm{rCli zQTA(ptcJP?Ti74%I_QvRnz|C=z*B;rF8n^`)2)Ci_S$St+H={2%B;+Xn?>J${Tb#w zA?x0Uo_)cLQMzJ+Y)kqyYTu9EIQ19xSLmD~;ZRX|fABvbA^Av|$U%+f-4YKR7R~`{C$T2~1yO0qU6U1_qmn zL}My00nl2~QezdB*4lspZtCS39KmztU2oJV!$)Dbgo9zOT^nDA9xuWVr{(_6m=y zOq#PP-84UHaXrNQxyoh%m#gaqXbY#wb;6C098Cr*33hBs0K{9+O_vcanya5R7O91r zq83k=ogfTM7Yq zXVFe-(&^UN8m+7k8Nr-i`85l?=$j@_KfL-d`rBm#kldTXfx)4LW|z4XfLHgmqmhjV z3N6tlxFNV;5E8cWiRrU&Z`q;Ma69tN=21DX@p9zwxO%p}|D95e!j>vVM1Qm4ev4&8 zH{dP@j@#=z=-v-G&Mxc%X6nX3(iIhOLyRX{2NIHGqd;e_7C`wfO*iTa?$k7ydDb=**tg| zzKJ+On1&vyeHl(*eHxKqmv%Pily58jQe5?-$iSZ>T`XUrmCxjU*Nm^SBSi5Rp#k9S z`U?*hll!bU$Tf)~=RQ*g74i99&7+{R9%8-{HwAJ~T0hX|Ak*kw#Y1}NT*Miw#;e=q z4d1G@@#9_1K$CBdPFs30i}@obK`+i!m#{fKRd!YEP)^XQyL%KoDDO>I3W%+@t4vG>J$>bv{O8qY9qpc0cF(4yMat1t{%#HbP%At{j(+{nn~Zdi5f zhx6U&2_{duWg1{HQ7D!QZ{zp{9(I*Hp2V>~?n${%8Xhapz4Y+c_u5?~4JkI+XI|#a zEpKB^Kt$w5IhtyQv};wHd?7%IhlWVLS->E4^4VZ5hHl{DdWokX z&+l;aW=w^NEdeCD^{|U85&1G>ekk#Spnj9^%@Gd*lD`L~|; z9QKSJWML%J?dgzBb2jZOE_oQUqb)5ht0+y&07VKUl!Tz)8s&V@z$59lBO=XM{oh05^0x6Fcmj?YP{MvQb^;Xr5h8ZTz z>{a;7FN-^orPVV;viQhUbK!{|GfhU@?!6n5P-1aKQssA9v~kxPQoh8^&sr%CG6G0T zmi(R!nduCaLZDWtM}-3cgIjkCH78`R@|KbQ;RD@y@b=*Tw>j7gY*V8Z>S^8Y9t1c}*Nnr1pPSlh zelPJJX86=84mW|y)A#j8r_Pwe>4Sb^0+CwQoYuhcy>;iX9O#}NW&>ZHRcdbgQCU?M zaip+I%ya!`IGq<+M4m z(_l-Z4gC!LnM@8|hfwtt0>1>Y!R(t6!$-?=YHPe}{Lli#E}81$>6VrX46{~Cz}01A zVyDf(+mH196Ap=v7M?oo`ez3pXEh`XVNGiz!jiTp5Jt`n9jaWDoup;ej`oMi6+@^W3}my8Ohd3UF?_fm$cow%Z{Ed6fDpF*#FEDw@~` zwL#H{noT5)?5~Sy*?f6X; zhMe6ft}!&uA-I#Lj%iV`n1u-4p{Ko+>NYq;0U5Tq#s2YCIr^g&Q7gFKOn8Pek@?u& zr)Qv>j<7_Q9+0N(JBYWH+!chC!X>e)41;FX(86O4w|EqzFMWH}rlxqyflIwkGbewA#*>iDsENy>qABGSh zK;JeLVuzY*${D1X#MeXr2!rPG?dEUVvkoz$6nwHoY@lPnE`#%ya0Rs((FUNy3-%=H_7M`D1 z83UbNrVO>iu2&s^?9Epfhn%czPK0Wv0Nj7*bk)uq0Gi-T>HT;wF-9WsTEjK|?pNI- zc0OKKDOw6eagOAjip2AdMhT@f7I%-vTIeq`UdWf&_Xe#Z&`jU#wwxoeIl86id)ZB( zO^%DVWz}%&OaFkb8M^CQNiu={S?eJRp$dVu7U%s6Tqt{Y+ayWteWNg*z~ zcBoi?)PU#(b6TgEA*RC5wGXu-^zjXF4Cg`*l;}!=utwRvN2Cs|c+B(tA}r|*i4PkU zd266J@~pltgB=#CS{G3hkYp0LTQkS!k$oxc-2xygFT*Ppd}jn^fkMqc7BUs{^KX+c zQ7?_12e^XiJMrr=Q@kHJw)J6oLEIlfVwIPx4mb`_T*{ko@?;f_0gd&m8 z@KQU3qmv@=0%($taY17cl z--)y}-gKut-ip_Vif~mX8rX0l>kh!=9EiQWZkDNFDKx_xyW7NU+jaeMg$sKBG|krYUsIXFwTU@S-I{h1zB;g_^k9-y!k}w0pf^%F%{bZ%RiGeWE12ZznM!U z3RXk9gb=siPNmH;#&tg!qNXM2X-`8M-6K@NA)Qu8 zqwU!EdLw`ZBdY6{k$!w+@8>`w&Qqorhm{0Tux3pGx z%+vKOvhwYPUP!8Uyo&#M5OkZ%!ry+3%%m7g7 z*I`*ecCfvrE-?ET*F9&HpVgm#GdX+TCmq*pi$J7?^DNRlt33pf-mZki!u{Dfkjl|W z*dl2>mbq5`y{g&e)QNtb<;TD7hw#ksC0A!NV60N64w_{QEyr|TpiW7+7($P-Nbl+5 z$2T?iK4nJG$$Vc9@Qsq^^h(&dVLkdM+mF2GG#ct#nz|GS?J4SjQM>&QK#K)@ED3P|Ad70;UD4Z{)$VmcckCwEX*NAZ zpe~BI44YlwdvSVvY!(K3)@>esVa6Bo=Pka%3nwP73_L!BvN;~z{gKT>nrmtl?0+)) zn-F>8^b^5oO9L+wCdw(a^UCS+=%h(wnVb}gysJfP$Nt>sIROEiXi=bByogs6-MRBA zC$m=Ud}&I=PQ4m{RCHnu#A|L zs_m!?+4R1MLOQ6kky$H#+rdcgHDp_*;ai8f(+z@XXrBZ>>hnQca;yozgcwZz6A1ov z7&;5KFd9{^H29K|x=$u(j{=QV%;reN7sj{?@$DtJxoE3l+KCk9COuwPiecKNd72qR z0g>$#aD^FopXM6~1*;U+1oz!&iVPr@;*U8M3($)=&zC-dg^l&oA&u7R&@5aH^kHUCmhq1P43t&s&C^%DhN`Vk zs4~PCaX0J@YrL_j5Ce?oK9M%5D|*wPU$tk2I|hfEM4Als&%hS^ zSgcfd&hTml`dblxL|;dME>tc4aqnnils0eOHH+HzPp;Tv^8msbDq?tfDzk{-TU7?t z$>j`oITuJmXB><|_bMAvNYkEICK2*2464y~l^Eck|5Nm``rYVb}yW-X1LF#fl5#{bHVt7}q9YR#g@cX*=u6)lZ3+ zr_HRC4tf8uz1z9hw081?zz~;pBmBIh<@9K29N&Hr^}`HPv2DE23~-y*n$!yE-#w*e zqY{)AR?sl^Q5In}rkTk*_625CP_*f&5MaDOGH_Qil?Xd|f_qBMDJO>CÔVQ>B2 zi4OFCF}AW`aU$ebO@q3Rw+59zyed3?=s}R+&%awgYH=ta%3Rm#BO!FlpGzQs$3xH- zoG|kO$byct?G(LKY(O?g)s4OaU=?*gP(3TbehTKID?C~$NL}WevwAd8)%131IqBJm zW7QZdTG7-6^W006I7ZzuLtS3~X8yA?3vDCy=qRP)9^WbPEAlO~o%de>uN638@0ng+ zk5DlSSXu}_>eCp?V<%d5qpJ@cpVXI;jOEjv@;`Py#6&UNW39kJ2$l>|U6~qR>c7#3 zjApeyITwmcFf$4`yKBR0t|0KNJm}$d6HGOgzx(Hu*SQQg1lemz3RrWETjrC)1N&jm znEE!jw5q9kqirDPZ20Ah+EG6>#^^Cf`!CeB6DH(F{=WUX!-?it(Vg?1v;3sp>}rYM zmnlzCRkTU_)Wq;3@hipuV(m4@tr*35HwEXFb%Y3ce!yx2ANBL0mf{}KB$G7 z3>hnP+uRcA>*ne(2PQjWrYSfA{Wzk;S^bwHU?4Bj$le|dqEUw0w5)+gWXTw-Iqp4; z)1kQ(w(I3@3+I2?Mt+inw`lZWM+<4q2oe)|+}N*sSt9e-LW-dFF z$lL8{oJK{^O9kN*q8h)JgsrNaeAZFrc+QjW8X=x(;O@7r|| zE~gQ}aPzoW^H^C(Qg(TZWpzZ=f~fSOWnPe6(e2hgRNUXKQYYSWM^t~WZrKc-|K%DJ z=~x z;xaiPuW#4mZEFm()CDQ$PXYF)Dt2XKq4U4eb_v z4*u?Bk=*|Q6uf|%Y?5O2h>u2Ca<)|1?jIVl3d38UW*3%cX}c&(KCA17o9~<3_*s#}iXb2WVak&F81=Bug4({CbN zU7ok+H}pngPNIXZ1zcg9OLqb(^qmGTBREa%bXN(j16@7!)gjQt$zShz9{*0DdhdGmHdI;chJ*wh zxnAd49^-6NzJvMKom-*7R4ZlJZ=u-5Fidjp__n6euHxmmFk?+&t)|3v<1DQiWJVp> z3!wgBpYJnMUx=p`hAgrAg-Kpep&d(VCHn%{x3T9vbhJQfioVsIHZeV7qNik#AP1{Y z;UcRNjytaEyd&+nMA#sUtdjl*t48(DIsS*cGHiwy`KNWMEdxuNe3jl-x_fgQ>aIdK zJaa!kL;u_-EDxJU-osFzDb71*X^73>O&IoeCzs70)MP^T7;5^R!Av6;Yh} z4c}rSti941nk=jyyit|@2H5$e7cXQHZ#Vup@l?C+F7>h&O6AARN59QYzKVaZTAcM= zkLR={k{CpJ>X?m?#dHx!tWgP*o0iJU^^Z({vjfB)fDjVh?&&QTfIk=-V57frS|Oc% zDvJPia%Iao^v(x{8wOG?N|c@3?tq7+kE2MG$o5lt4K-dmVvQOVcSLL6iPj>#}3XB+CaWPsKur4so0-B_kiwuCUrd6XoC@urasW4z(}Bmr~_V<@cjIF^7OQ^3!pQm zMpPU*Su*EpxLOJfurFV)Y^alc&Th-?)Bt`ke@kymN(Uf;1jE@4{%BEDm|*TN_y5{I z6h9oUmhqb05JoXX!X|YNS(pCegc_b~74=*nQn#Yq(xg%^k{S5y!U)sdD|tcLSp1G^ zLmJ#44`XjXnF+SvO>Y25C+iNRV}M%%+ePm@R=-p1Cj_c8@!-=ar*OnriV?jjBfWOJ zSdtDkX+l-+1BR%PN!H&b_~$`yW(_!#gfbXP4s647;F$%Zt|06pK=jkNO-~c7JW6`5 zvNSI=-lZ0vSdK4d1-|bkZfeGqM(NhTqe6yRt55F>YW}mGxg|Lu5au0X)%6^|It_40 zJEjrI%VC%vMrd9{@IzK?)6Xe6QzTP$rcsjSQNm<@rY!j{1RlA^bnx}g6Wk@j7}o2= z<$cKyzEQ1(l% ztcv>PKDMb>JK0Of{lopE9cfCWd`Dq>6=QKRfA>JhvBu^3cGgaMnmg-*bHiU9J6(hn zG%|Eos6)>K?6u(wqO95xgAUubG#s3gv!i9(I$r|P=1oq3r}7Nkdcs~aCEX)M+A8c_ z@l{sIrhp=mjDOsq;nU$E|0+|3nW=?L4jfZ{cKbc&ZAHPF3zU zEL$HbLR#%TG|D3ZG#Yj#GIbm!`(W_%N^4#TAq3H0Jw#AQ>Wf7B;Pelg#g-Yjo%M#GSJ&`5y6* zw(N=gSQn>+X zSVP4e$CHA^#kk1Q_3;L^nC_9#Y_Q)P38nJP!d@3gARX$XI6~+2{V@9awjN?4eC+&#Jg4MKFFkdg)~U|bz^qI|>6&BJw3`t}R#r$Pf0UB(l=mh;PV8fy{{mlka%zha&qZCR=n0qPYz(;+X2 z7fs1XZ7wS{`(R-~v7Mx97<#;xunN>0O{`81>7l(l0Hho`*rH%$W>s)qp+5X(P=wH# zHJjv?8WM%VY{}n8HH#bHR@!f`zMM(#gL^G_)iKJbuQWkM~r+SRW537R~_KSG!@c!MZb}6Cz#uwSdF!ih1$rxx_wGa5UTvp zTny7+c{h`lKl+Cp{Y{k+BwS$2XS@%yIIkEhB24+Lc~63E+_h%Ei!fL`biHZ5b|+Ty z@i%hs**E|QkGadP_nOVYl^<{-Js2Vv8tlI5vrbRvZW#DE=V zcRKs5t5GlJp$$KfX_Bs)GRmPzKRyblj?0cdO9L`Weze1UK-xlJPMw@sZu8?2iIo=L z!46zO55Y8}_TX3O>f-DTHib%5l8vxZthRal{ZqK1n7lFR@@k5oKk2K(TF6d0x|IU> z&E{%_uzT_j8Gu{QjA;?{(Si&{UwrT77(-L=R*dceCGPc)biy= z>c4Ne@f#<9>_a=CC~MY&fu4=;m|B-}1*-DF%bmGs+CuN0xVlSj7AGw~<}CaY$3)LM zV+zL2W%tbMzX|DLMJ;e^gUX@3`%G*`KyT*B(a9fY)1KyvuI{Uv5)~7u`&>)|4_zsKN9D4w! zQ?73Q-r4)i$BEa9BFmDIeyfo2_eh@j`(oN)LDmYwnA3P$`SM@-eweN$#lS`n~5Ff$K5np@x}ovo@1W+5%FAfZp#5^ zcU6VcXNx8?J?6*HiZY4`N7T3o90P$r-paOD3MvZVQUb*j8~ttbR7e>HGOqIC{G1bt zQ@&H%Ge%d9MGrwX0L9i!lYc^)SUe3u)&&)O_ne*%t%w2#di>4xJSC_cwKzCpu5VhY zR1_zdcYL4ACyPG2+cf$2I(W7K5!}|%$eO@5G6V0uTHEvY&V|+4Kl<6YiCkqRW;&tk zXPRF0{w;P+FXENrYJ6gn#8ZzNoPTw-;ziPZ$u>%Fyj|x` zGtvtn-G-n?Erh}KMUhUUl9`({OLnIiwwMKf<3$`R_6AX1up?w>rf5Eb{H2voB@;q3-?Zy0FpQh*${vJ6xOM7Ww8 zU`Jp>esO%C!l}Mz4L$(#=2e;*WELjD;*z1SJ(_mh(A#!_Ke{v0&_0!QeZTzxpAAOh zIW+4)*%>x^Pj>nlPh5kx=irp$O#MOPycjTMDSHT}}Aa42^Hmj_&=A?+l(sqJ~Kie=_s#uwjJ_nMr=BQaS`tCak!XT*F(EW0j} zzs1zOaexfW;+*UJp%q^-WLSe_*j!$ucUL-2XFUOUEI zGcLW~t9tc2&pz3yTJudFNV}=Wk`6`p!82+1qIfkA#_#n8I*22vNEt#5?cw*%j zx(E|$w^KF`iae&*!ur|AW@T({y)2*j2?>DGmJ&Gqf%y6XWzjWMd{2*h@+pQi)uFD` zc~H;xo_}3n0L(+V)-A8GNz$ddJ3@aqFRa&E>gCeR&|6N_KSXgJ)_Xk@#YWOm(ylQg zj{>hjSDbpnPjMcJJ=J%UCVwb?fWQlcYb(Uh`S-EHR({AIT!XooTtI&C4+UCRE9a4 zhKfmr3oE*B`np12@H#b^r##sL- z*<6(7QY&F`T}b*P5^rYvdOXF1W=$?#L-_NseRT;dcUa3mx=;`%lk+B5Dokw)u4R4e z*Or7}fHVIWK_M`X2!1dj%lApXM$CE*z76>WL5O6>%a%CHtlb~K$@8A-8h2kRj5k)k z)9zlR)ju*}7X74&Eo=KcuA_(?!MkW0i5WKF`8R$uBBHkFH<{$b7}g9KY`C48i5+&)fc$NT@f zp`_g=C+Ll^txY&gwPOk1Yk$?<(^uTVD5>?J8(+v5A z%r>JPXO$3`#}=vt1KO_#{$%ZfRA zmP9-6>Yr66pYkxSQ;Ty9pRLjf#i#$kwdN6h`GJFQU!A#QJg`Zf-HWh|&-aY!` zpEEiGeXuZxTR%m9SJ|iVKym(uKZF~rzbo$96A@esrn_G^6ng5|eZ%uBjt}^hHtCD_ zt*JZ3W38N{fYtp2;~=)tSxo+B*yriA@)SLf8=G_q9xH^EC8rJcE8~C?3P!HNc}FVr zHg8+KUWgJ6S;yRnQa!`}d0_}VrB_zJxnhc*9?&Zvx`y0R_;ha9oMOD5uDD9!F!wcm z{Zbrft6^9yW-a(s^@p?c0(~F#Ma7g~K$23Rh1bn~K}^nA|%zWKiOh7@_Mf7#%X!gV&5- zk63%gdCl26ek?J5rcrt+-KiW>Em&gOtsd7#@{07n%0^CP`&PpCwdz~nl9>LVD}&s( zJ0-&g2oo&MCRU)b_>-fiY+E-s%7sIJ{_Z$VR1B*@C(T^cNFtp62CqtH2a=d3HiLjZE2qy zctv51%izzPeHn%OStfm&?Zp{B*beiW#Ll(&BgmDpoV5D!^X;W?RnaA12xhRhOqJ#V zPfS$pb`9!FlUWItl=sec7>XN0fxMX5^_2tG)>Lw03N2hdSaw6D3MOPdrmAl-tM)N_ zwTDB&EKuudlgic3^5(~)(j$zZRb5Q<3Tm(}B}1KH_5>da5r4AtcApvq{f!-*Xnigw zXaz_YWF8yzG(aLh{;gwPe0n;^xQ&~|Dl^;bPiJp^BF0kyl2yr-+8xa zODQc{sk40nzd}utD-2)7FoF8G|>oA*R3C<$kHCtle88)jNlE`@Cd{DZbP*t|u&0 zm9VQESO4Ix5oew3kf#2D?cGnSFKI!R1J5aY5e0R_#`mrte+#x?nZ9|oIq7wK;KfPv z3s%%kUG)_Ia};{{jF!0jIxW=dH;ZMbH3_z!mJj{$CYkGd&C}-8oL;Fz;>i=Oz2CMm z2h(``_Ls8QbkBxr2Ax(>b!czCGmnz6-euvJr&f4>^#D`H21jgoX8~=o0V9rOI_Zb0 zx<(i!&+0<#4j}{)R?mJe2)5W&m-fF` zOO9cqT6t7od!YB5+&|(OmW-%eUAE4h|LQ(KEr*+=>19JM7QTJzAiBqA|EBHU3bVx% zq`_*?ioyPH|D%>n4YC%8n5|DIZ~Z}*m94D{_kc45w?*d4Z^W1LzzV-+S9E8QvWBgq z0!OcxsHmU=mV%m0CnqN>Am4otAFiyNGNz3v>7QXIUjHRueX2~ObLk4Uf?2l+ z*8lZWe%eV0lSWvs1GN{^*!yTOY3=rZd%gcrt7HXt7^{uRi-(QOw36bL^2SJ1b>`Js zO^66VeBNzk=2OXXY)gmL5$uW*2izZmca`jrRsd#Y^}bo<#h|)QJhU`Z7QKs}pZqaW zcnD|~L4e}?Fqqx*;CDs-gxm3@*A%McxK_@GMG&ap@owWFxjCeSf1SKfiHyeZ`}4gOwC_51Grt-9 zP!>CfY>A%BpQBC9`|dR&8JLw-B_83a0+r$dYmG#B#NJWftWRlfu~9K?*cj^(^e5+* zvtn5UBPIMce(cAF?ZR4*^Y*|}(uxm!3T&7B%W|@m^5sv&hm)q1kSmkTQx`t5`sppN zRZZx}sJLgXnw?E9RLHVd-RtHT`7tA(#Kitj#|SxR#5Z;v-xJDBI=#G zE;4RZbs&7;pgjGi)MfyHbzcl@VI`A5W~jgv<%La-;4H{emTpQji8 zk`j1KGO&e+o(a5#OEHO`A=E^w!7rS1VqAVk0%%$Hxvti8Rd$JTTtybm2X?o6j;*N=ifd}gyY3tcZ|e(~7p@%frw@juj`yx6AsitZ<5^M9&~KrC6JWfY z6CiW<>N{`^O3ZJP{%ldYVT&nlL^>;K&B0ARsZrAe5Ko+WRklOTc8^j@`BuMwA0+T7 zanIVaXCL<`@}KiH`_7pa;{L2{SRPQnj=M&U*-R{z8aiDBYU*u03~W};(Ioq>BEas6qFl^ZQA!G3q zZ3pTOF5CQ1I@^~K6K5qHEnF{&o%xy3b2uo6ElH%T3T1d+M=_@n)f-e3z^Fm7@jH~F zdKox+xv0>V3LZL1>DiidkGv)B{b- z)?01>q8$;*u(~(9M{6a(#|j+$Jk~o5rOFd^D-Y?I+ zN!ruNa6eM)xIX;r@x>vPqg|ApYGuquTG!Q)v7lOnM)%37y-!CJJ)YIa`zZhPv1jPD zJ|ABU?1$dCtPHkJuggTNVeJq_(Xow`u*1FUXZN7ABs#hYi-c?_c()|`mOuL^k4^YnEBBxe-%@*>5 z%09%?g0hE3vd{U7XSKa_6}z$GW#+eLQa0U6dn{qMtI53X>v`O8oxJiSL(7=cax`mP zN>QN4;iy8#KrD%PEOEGs;*j$-_wflsN=k0n&4+H{Ba3j=K0)biX&*AI^-5#LQ4?do zvIed2F;t`u&C^*9Y37&MI8ImUC445Xy^toG8TDvLN`|Fe=oRai`4w<0Y;tg{^mXfMrA z=f~jSq@2W1t%A_zplJ6K$=9l{rM&a-gH!z(Q@D`E-D+EWo7uY&o9I{|y?bX;cQS&M z{BZ5!^%R2^Q37FQn#h9P^i_Bu$+6<|v0RsTwrIikeE$`eVp>|BE0I^cc|Xdo+4q-E z)pPx}8Gj#FKTC*cmui0E*>%D3b^h4}>n;98!m}YkuI2zf*a6!Q z6OVEAlv5js8d2fIf`#M=UX?#M`nj)vdj<7P3xuTzjVz}mLh#;WgW$|mY%n| z^o+@-d0v&rlZGQ7B*ZM4BPt)H+v)A9dW5Ow-sbJRlf?|EDKS@igC!&cQ}L735&QSR zoVzzH@g%3+bX%Mw?4X*>i?AN!WVZ|{SoaR3eJeN?znFJrwy5qCeBqqmp<%a~r<+qm zZ+gzIK#<-oRT2L2pbuD7&&F7&PAwSmQcPp=TPjzTltb*+?W3RiOgECbt|zDJaR1KU zP50@M=en5HS{C8#-Q=Hx6jYwKy&knj&vg0#WUPWui+9#RSX_TzEeq~)P){A$sO2%e zeBRaUtNm4IpK42Oic|C?+gqh($wsw>ko869GykvS`*PI?s2WuXNtkfimoGB_F}i_j z$3qB@VFcT$$oF|3E}AFrc6lG|-Y$8kNUNAFJdrqq>RRhwYh3w?S7~*JY56P{gx^q| zUw102?{|09lmi`Vaqj)ba;IAK*_GiByM*m}f@liBf%`U$`>Mvn&P1m>xjK=+S|;KFdtc3_SujCq*cTV2g8CA%*15JGp^(;c`h##=v$( z%jYRXl3Nd)(_!7g`J1GRiD6zXH+&vOP(Jj3nLikrY317JI3KDYHVOBZ z{%yCN0$7E9>PpX`L23 zA`1oHEHj7 zzL-_qxJ#xVvH@nD5o*EGEg&;1DO-p18{i_XeVs>C-{&6I&SEH{CKfil=qV^)-Nn^t zUbF6FJ+~zr`ZWvINjZP>xk;LjYM)C~7sev>CV*$58XJ{z{$p9PCbeS1HhuL8^H7`j z#cM|K*~ijD=nDmo1ke$AOuqiF&(fV}hrnB9*7KeCxUH~So!xZQR)ViFA3mLqho#;z zAa3iGi5!`EuGs!k#BEJc-EQDUNGu(diU2GRr*rDYTa)KB(PpQ)H@koF_3us@n{6&# zqSuEuTTnEzRwM6?eZYpSE;?z{x*`QGE*;%tK2&MA%?JMXxnigRFh|sUEkDu~u|79x z4Hl|Cuq4S&f#*0gZ;z~O?8N8K-isawo-)UrlAUV~TBZVJ?LVxi>sEf?^|XULRg{?g zlXruD01*y=OQFZoHT35^rGu}OsP#3Iv+i8Lmu?@Q)Kg00p96ffoshMI`hYCA7I;Ad zXPnP=B%|6_-?Moj($ir>I;z-(&CR(Py53kp0vwN&Kpv;F3cxZO^Nw%9{orfh>jr75 zvo)wk`+7$=z9mHSasq_$qn_i;@i-W$A&!qr7=obZXzB5-ZnicZC??}O&&h5eE4Hh2r#JT$74Zu;TAW-KZV z4#!D98l+|)cj4i}CIS_AD9Es_hXRp_*#&B>0tx3hT-9iiCYmi26x zo)PpNR-E7h#c!9L$4`xG9?laz(z&Lt*G3XykT-89wyKynUwkovXnnt+_u}CFh?-Jq z=-9X9b%}MDq**DUCfBc-b@rXk%Bk$>+ z1WoE4ZO-Z9x(YiHaKUQT{o!`Z^&9w!EMYmn>PJEEq;>RnQUmVv-oEhfRe#E>_YDoo zv*0mos#Q)zz?gC-0OeKBz7U%>3Qs3kN?}a6Ru}fuct?3Vw;%&Qj6$GXgSRm^~K?2<5NCkXX$#W|Ym1 z3`+YSIavyPM?s9Sw4G1R=c^S4@h-YxY|+EiK+z&Kaj^P%ew05*^ICEf*VwXAK&jzg zQcl;tV%~;odb!XzL@mYO^B=d8cY&vS-(J2ld9u53@DYfIYs``9p4GocE~###zxqxj zZ~&#MZG#7g6j-@k-!|UY{^#DC-ck5ukJay@=v^$F8qco*n>cA*fL$4euU3V#8HDAD z$j3hX3u+ycrgPp%HJ^4gWh|7;OH*>z@at%RBb!ml*Y%2C66wD*Jj-^*q)c|TUWZIY z097c(37ve@d?bzD$+Be~P!6-pr3J}k86JxoblU(s+=NIkjQ zA@Y4p0>pR4@HeQ&C$5mz4V2?(x|v(?3&|Jih2*qvWiMEk+yss6B~6XxT@W54DHk{W zOB<)*A8uTkWQ&4M9(ETgxnv6!HcX3U==M>9qw60usVyiM+_zRFEVNbB?b?AIA>&+O zce<~lZPgXzJ~2PLn{(?l_CY#Jge39tFX^=~4xBFsODX-Gj7@?8A63cg^EPbe)Ag$r zoANt1IF}}VBKU$U(;@`R-+8P2@;(YNa#;`7^=J*rZVT&T0V4#e?nzY}LYcL*dc^IZ zWcxbl!xOge*G@-(+FI`_xaIQStcj6go1fiJrY*hYD&8b7*wJU=xuOPod=v-lz4z%d zJr7vTKRS;+_m}VEs-j7y>nL$yrb^u;H2^cDHdXHVaZZHHVUVrp5xh<9Hz^Ecoyu!~ z$!QNhd5{r9ajIiq;oy3&m7fjFLfRs)-TfoiYOOC3 zTxyly_MQYgfc@+Y8tTMZhv0f02cV?OBqerS^uT4VKkhE%M6U90tssG`MXmn6F%fH> zl{dG6u%p(}F)!r<1+Q!%#R2aVtfHWziYD#XcZO;*!-fB{-ZN9N*z^1-KRCl6e&put zYek8TzOMUu?m+Uzz@ALdjLg5U^%HrrrHML%1~MSXc2eg5#z*}EN-3oyNcB<%Gn)Wvo5Cjd% zOEa5P9&Seb7^7IUYDVgK%|)}7(f1g~{y-cD9}p4syRvytm#ePT7cG_r^ypdy%x;EP zI(oct*eAj-H0IkY0<0liT)gLH6tL40(O0w4-#?6!R$s>LC8z7O(jMl`Rk5skxyfq4 z-w5+1DDmpI5&@zfc6t&1#s{Ywea0~Nj8>{1bShHC7YUx2cnXen>V=%wGa1VcGlL^< z9%Wc~&RS@iN^R7lnmP?os7rt%9E+J8R!@DBV-|rx>SBBgrEb)el}LdS3oLm5o?@5j|q$qg7r$E)cFLKgC5-0>~_ussSg0sV*{ND=NJoB?Y=h z^_g!$K0a+bAARy8*T2=qMIEy-n1{tlbQmoprs+m}$qE-UNc2Lz9q}GF!n}16iI+F#4i-TR3*iXK&^_d z35lV8QTuFn&0N5%2e!hfD|}39bV3?UkW&wAHCVivK#Qx%5JM79V&IL zO$75B6HI%;_iXij3`CHK4J64MRrz3UDKPx9xSe@V!{npBivokF>-|4jPmKI_0iuX0 zzG221=93!%8W#V(kguPc1G?WWzoe}Dw1U2kj+KKe3FUR0N;YNq1P9+Vz_TAcXk%O(Dwyr9atq>Ow0&324T5npFF zG0=s%x~`#-C(!9oS}64XABC_k$xarbany*~`QK%~0)P1aME?YtKJtoO#+Svt&8BPb zgqS_bUiQgnhvrSwf^CWb1y3DwZ8B}p`ll4g|1o5Ph6iOU zrkeEBYMWD%`@KL?Eo&F)?v?zLI;bF0e>Q$)E<^X{AS$_#(C+3|`rI0S-5wkP^|u86 zxOm4f%;J)kWG8Fsg^$QgO8_lC&K4}@^}>Cg9KmN!O?hIKRPJ6>Hs+s&zA8=gr#m|b z%u{~lSoejjIkza;%u(MI_X6Qi%Vn#!DfvR-CH!YLIwg>6?J2h|T3vb_p`};T+fZ#= z1RBbe@@YMl61q%_rgC2$+GO+E8|Up9i(P`NG-utz#68e#bTZXK>V!c`qktKh18LmjY@&fMS}-_ zFY((ZyEW^#&X|=(Fw(skht9TpN0?SgtCFZYa9bTDlS$TY&(-2@>&t60WW2;HG))cC zepK}e6LCl3jn23g{7p^quF#n9_{>aQh$zi{*13aOdw2|7#7k%OENzTK@=IBqNcP#X z-?pWn)|8QRi)?zXJzrwO%G$bB>3Fil=9wk$QPb+3lBjG&Kn@rps=mvMX_Cz{F54{K z_GgvB)VJ7QeFa`;J?x7nHs?dna5b-7o4+gfZ>_0_zYok>Fe-+YS4Jsvfi!Fctk1-8 z?8o;!$I;VATfRt0su<^*O~H!Y=C6SKb_sz9R~K3XnU@>j1<@YE-hwxeKLQgLUgg(F zV77|FS}iLG9bb9|5_!D1nn?h(Ip#8P`yu^zMcqnw^4jkt8(%CiOR--TjNQ%`nb@p+ z-kJuraqBU)4I~3Tz~g6aAe-h;PN8JAR6@o>62fEPVMk1Vwo{hI=IzK^tgP<^I*7Wz z1x-N4?y1sZ-%)mH34jf~z*y}MNH>3~TsdLn<=&kp&1a;ri?c-q>Ssl-I!M z8j&kO_xR9QzgRY-U!@pE>_S24Yw%@|`|UopH03-uaFgb-E%|f4rnbiJq+}D(jZ~OV z+dfYQ8^Nu_s)U%AL+;a%6^8C;E#8@DX_qONPw=75#UBYe$8{Ou$rQf^c-5~!oa=om#zj(Q=i8^YSe^~XEv^Dut~Fo77u{+d|LQ%S zbJt%paBV3_vi^-Sh$pbC6L5n`%Gw2?cLPxK#GlHcUPNuxtj;F}4;=yQHPSN2-2SG< zp{u^ch+GmQwRQtVK`ugv!=uLp;}X(Z&khg=k2Hm%Jw14f8eGTOk`65QJ8R?jizvU6 z3UNGn)3B>UU@F;%5T>$YS=6|ZK`ES$niyi^@()lsSIZwN0L;^Wpsn|dm6)ob0QM|q zT#GN`IJ5olO&#I_QT@CTAqKlr9v{#NH8FNBAOoh(mWO}|uP3Yy)?h2F1)#T`{&3pe z)&Np-`d%(3)i4Q!6KnJe@g44(+@(Tm@;wxr3xzM6J@=FCMgip)@x~TxMdWv(j~cM= zp1(fHwgp6r8wXe!Vb2E*TI!^3i&kP32{^?@j3ed~%As%QaNGSGI=~%J|NoO2ZF+ES z_M8M*Sx~o^eW1{DYMUk4)Tm%G zMpf<`nk7BetB(CUf4HAMjjY|uO88s)!A^&v&jJx}2qh%)cAkzr7-s|Ymcs@~On?9` z7!jY#=#$r_$x$A@IcOByVe!)y91ql#D3LOkt@auly8^D2Vd ziXvp*&uzbnSP>CLZ8@_?*H=hTfsdTM>*1#|pLkFCmwgsjb1c2>BJh=XwW{YJukGv| z03%1o69`;cO~#uFk0OcDZTqD0mWbtOwZu_YFzjK!q5}(QBzU>caWQ#2I(ovdav!3W z2FS=mkE5Q$ZP${xDu}>@#Pf5+orYCMEt*~Pi7bPCJ4|?-h&3H!W8#;QD=c~JHHmVA z*bcb;X>K7qrW={L#O~LbF$pwJe-6sV8=~@|LEeM8ILtFsaSm+1Z~4320*CdbM=Cdm z%XoWE?Md1i?6<|bFe=K#&ekQ<>31LR)rJ$+ ziI>$=B4_`0iiKeAr_J^#{agKe-@6GZFd^ev&({=GLpn}aHa{E}{+@(nD0v`}(v9Jg z@$8LCmh`4>gS1%_kFw!>gQjetpp+7!tM6`54>g7F`JVU79sXG`tl?`3&9%t zAdgF~#Mw)3Um?3Mft@Ebj^pwANl_(rdim(X&`(Gd-C zOIPu3GW#L@iOA;p;^1!!eMUQhf!GLsoz{5q*FqFuSjEs@kr_s=tpxf(nStmKgqrK; z&-voyn#pwBZzJUBMdJ_2#SI<)BF#%pIbBjY_(kC0S#JEJU=@`be3i-REIPqoSD&kp zkgIc@pL9??Y1F{and~bpFUYJd7cH(mia9 zeDcR5%#v?Wis3cumhxEN5g;9jtt|PbGpygf^J?h9ic^+UUQVsBDucreR(CA$lVN|( zV4_wVnhc+F%z`%ZYVm#v_RT_@UgYvVdP;VG^#lyUmV`x|kM;pfhk4V+p;DAdVGQzu zFUXQlOD1_AfJv!+1s(#k?z?$<=aCCH`QSJ8ib9e4o=z5=3au~lR__>8Nn8bFbe!cF zLSxbID?LQ?TN7`kf{j5lS@aXJx~X8BNXx)N>2pJwv#?M7JBCi~=%*o<#+$F%EjIn} zodVwlL7D>**TmM8k2uE@ib^Rf(rEl0JM-=Vp+E9vkHgU^#@@ObrT}pU^U9+}W;6BH zZ7ka-+ZXu;k1tJCH-B)1v~o+??UJVdF~$AQ4-i8(fp)Ncw^N~$#B6w-1MVqF=E`o7 zoq}gKpFnYp{sI`|DY}iPhTAa?L$Wj;<^K;>hS*?t{?Y zeX!eluzbq>R7LP+7z}i;8v8)mM`+j#IbX2-fPW9?ur%q++Y z*zU=RxFl@)>~#)}f#|$1RdRvNhG$fjoMn_**ArpO*UI!B{KlsAe1MUDsn1W_s00RJ zFW7V^&X5J>#1Of(D8$4gOR>7hhALp!F?!@UwNefMEdbKqcAa<|K+$5d?Cb#OSv`;A zF*~;u(h*i;6^L)9qgMA_#LWo&#d_Lwo^qH~t_L1!paX>=JgiP09iWKNzHO9y@1Q_3 zp+R~l_r&~WBZ!>%PvjF9OEGZH?e_<#Q zkB{B_ubnXb@g7+pPElpW8xz_?Fq#)2nFRKK=+=h&tos{14xaP;hg(9ep(py1xkEPE z0d9fC`=Q460^9l9)pnwJRaHBk z?PZD#_yw9blAY(J?*S>JX7{#C){N>?>Sm6FqpJi$fOl!G8TB0QX27BYYX{7SM${)P z=1Klu^f2Ig>*L*7jw9yk#_^O%{L^JFgdrQQC1gmv8bCCXuh^w0gl+Cr39mNT@Z^_B zIZ}wbhp(|i`+_BbMcHG>UvBMwTc?I?UzQ)As1KlfL*<2SqFGyIbs3T{uAG>e>}OBOUeq2|OG`f7tUS!=3ho1`rwoKHu^7jKhWeAbWmK+r5-%`NSNq@8Y~Q6Q z4O6RR^EQoNNYz7Etd&Ifz6@O~Q2NmLi}LjMICI8-+WX`j?K)|rEBG$kL;&fjw>jV} zjh+6jEg)Ck`z>RkqY$9m#OnPp@uI9JDMe(dhtCsx=9uOkJXO2Llm3oB+Z$Tn(RRA01Uba?C{5iYYVMpEw$zE13uWK-_^u3Fb;q$=_lQLb(ulg_*#NV0&o(AEu#Jlxz2O5Gapv2qaI9Ef@T_9s*q5_Ef9h4DK`OM^jTT=-ySG{Gq&Ci^ z{yRro$9Nd}B%8`(P9uGMMbsm-#*8FmjNd+uBGl0D8BYSH2d_fUpCx`9_wUbiNQDK4 zvZ{T}9FEoK?l1);#j=Xp*DpjyYZD7;p{r}gej`6V=^15W@@qZD|E76rpLJ#&<{XUD zz3NvE$$4iIt9ZX9xbt*g*kD()w*sukRWjU$-1J&(?$SeSt-zhce^@tc{aqL>jSb1UHls;mb^jc*XKp(W{7op3<3q#wr5^!{ zuVFE<@~+m*IuKem+4w@SW0vIcEc^Rr zsv;xrl%Ij!jVj3wZmM}lD`R1n%c6e%Z*T%zO(N$XHvI(TS8YrpTLLoU|&D=l!-duY8FhE-vTSkvHoaY`M;_-T5tfC%a9zuO<;O zV`b-SYi1C4a1$0iuC;8qhFr#-5AW;APwea&Fdi4)khE zrOXQITwKCFo7p6EFqwNFNctF-$)(?S5n%8w0SKnxZkdl$UN>T`e56W{+68cxo@@En z$7yRFTjy--D@3n0E7S|46+uK6xmi*9qJVU(sr)e%L?uOE?bFCbw9?Pcs+FMh+o)0a zpSj>~ji*B{6u!7GEb&M6Ros&F2l8kE^?|C$ZScy!_xY8ctaO%i_T*bOGc|w5^YQi` zml@(zHq5S61BRE;?~IsY_5x<(OQ1m$rr4MVzxA&>GUv90ME(sK zJX31>#-%*)br<^qE2v0Qv^rn&NFcHB@!b@6@|5w}6oAxs{x>pPJz@Uv=#GXNw!t93 zCc#JN(6iPbojihXVmg~#Tx(1~VCOW4*Y5h9&wH)F#tqj86KcQ(9`|>W@^NOuae$dn! zn$#GDJ0w$uFuNThZ%HzGuv4P>&7O6B5;f6l6@N20uV9 z@=MBQ0XdOYH}Ds*uSc&~#EpDIqbAo*{8x$=KwYJYeCU5JhTMPoSrvIQDB2vV>1v3N z&nhFapA!hbBOX%u@&HAhY)sdY!A`tfbPmi;BlDJ5R;!R1!D3;vmDSg0=;b>TFO8mW zhUbc^VATj4?_tY(t)?25(r@Owjhj~w)>IWVxNmc2&fkN;LvC&Rr>;AwNad}n(D4j$ zdeFRL|0Y8lyD}R;Klfw8Ix7Ika>B1P-M8r&^atw-tZ2%3lv0d%*lNlrTM8qbeR<3qVF$fG8=p(*l%0*bjKLnl^bm_bmx36 zqjE#2^$zmXgS96BI8UT5b(Xe_KUhtsatIMLeThO(T)DSEv6G@wE%BE~0pF39Vb-3? zx%%(=e#|UR`xOMYSlT#>t=(>xpC1>7n^Tp@^o)47-S0oYI%J(GYKV2T!p`{&$zdBd zHzrg2L79nixFqOUdM=IR4MBnw-x;Us)pldg1*zMr#RJK@8XW>zCJjl^j9|)PvOHZr z7S_4TFrOO+hPKpC5}XruFTg;+yOr>JpXWeqS`{#gIUw4R=J)=QZ?8;g!3T^Zb#9=# zgA()SqL%8d{SoEE5V*{?5!-BG)areOcW^sWR ztL$p%Bs!i&n^(wg*E|?eL|)od4f|!lmikXDAZcWh9fZ@H)PY=%ITLuwN;292(_Kzh>U5B?E~^^D$gsY6@XL zmwK$JfT|!serHxQdICfb_(!kE=;QhlK!`a$FT+(8bj`qf?o<-dTzbMnh7lNf7BMDu zQ6p4XpxZv(H(p#XSHibQu*M=>P$b&AvJObF`?F@Z6LiyVuIikQz4rqlmVSkeF{s;b z6*T{`k3i=!A_7#t;bUth^SgZ_!jy@R7t3!_gzy$TsK$SPmMAWCFGI;Pu3R) z8zEX$?@sXz^@Gr8)+ivd$%a9yJmTOs13&8ocy8*V#RmD;4}K=@ zMPk8AVLIw5bPhM;0FbpJ6Pzoj9iYfZ+*-R_R)X{8HUPPSnPLn8l=h#fvK<({CHO(py!RQ)7kr8)U`K$A+(ZAy6~N`E{$FiqKt&)q3j!K6MxM@jN7$#cghVsg_P zxkTa+xryr#d{f)rXOGVVmxmi(PCQj4^#7^1E1L0z0MzO04LLZ)>x3#i!LswHJpuCRDTfP2Xa zp>yRYbg~J=kxZy7GGjpSLW|=Fo6;CU<=AL`(0A+p9_h%LrrIrZf68u_FK|i|C4wu#tW_%{Nqnc60 z0XGZ#1cq?}EegIvK4;8ryrglwCVM&C>-O0ZNsg8>aL%;wc=((2d!jvI$r&11zA0jL zunJINFMzU^+gD&i|m9rtt%(kmVML0|%| z<1iMR?*rrVzejf3w%W0JdHQlU_^dXVDX%67G;i{e*l6BpcV7sRZmRe2hoaCb$@@Cx zrttoy?##I4az(h+_07urc0lG2gx)Nnki25h*??{;NJiDf1be{FkzdMb zAMda5%20=glu5|pEpy94OR)nF%^ z6W{IYR~t#6yuH(<%qY7KK>8|qhl3%<`5j%94_>k%hLHWSn8Ryf zExj|NucE-O!aBAdTizoB3BOz~8VCEk<$*YGT)Pi1%q7>1C3#=_pGS+NBQl6^GlZLE zRB&+Q{DYsmMHvDeLff8pi+EXCx-ykf4zzYxYhl3H;I?eiCF5@&OIs03ci_!;J{=|{ zofESh#LNDZ(e4fpb=>nlk326;)?Crys#=`zbmc`%jal%le9PVoxfwwYn$3j-85V^9 zk#xwP|7ch1^ShHcx~k)8R3F=-wXQr1e^i`VVU7pCRTxm0qPV*w&*yo3hRs5V1vG}b zBxg{#G#Rx8L`n)gH(lrZn=5$s)ZV~aaTVF_y!pQz$>-cy)C8!ppUt7$Fj5{4kFfQy zVjD+=v{oiE36{G80SBSwgf9*GJAnYd_d6LrelRAk64?zL4Gv}QW7F(a?tCaP zcZ$FKK7W`82y0+!;g^kmG8Wlq1Pw(SQq%W33qVII*+o^nYGtdm0HtOpY@>iH$0>zm z$njxlx8d_{Lt}%Ni~9gdvcz!A04yVWOu+=(x0hmDjJI0S?KDzT*SZ3kf42eCdoN8^|(}z(^A~Sf(34iM@N9HUS63PGn~Wdm9@Ja zuCO`vY2J-Xz#{q$jC}rf3DdDnWazX7#e&+NjbHFBl4yiAnTvABRis&&YsxCgen)LG z9yT%_Ky5?GxyQLZKo7TvGjuu{Wf0;wZfM+GKP+$_md4bJWOot{4{q!SF}m?!C$r8J z9mioMeD0fbLNyjVPl^sv<>uT|%S^yd!ThX;J?ThV?*No>Fhm$I9qg1Ew;%jB3XT6RYDG7v{9?25?E_=J9gzx-98OxVJ{QLFy#c4Uxs` zx}6BjZ;{OiY<*N{RD$?8?QTu5RndeS0cT9n+qu1u@zvTbM+~@ZU$<^3b$95jRR!(J z91rE{~?jtjO~ESQ^w zpcTLZD14-wXiud)%~>SgD}Ao#N8_kfbE)*2eEsd{gokUPDAn|%^jfPi7qK<*(Q|7f z{g04(UrW2c58V@68=g(-8kbshFB^6PbT5(}3!!Sm24b~oZtlQyLuCA0y^Jg>iH2r$JM(d{M}6XhrGedo~Y&w{~b%=@8W< zh;?k}_SM?P@Ln3-f?SR}VJkTjEU?4`AMPretIdEz&^!iBeY)K3F4y0`iTIV+1~5Q0E}_EKHePeYzTSqM0RJ{C`(5XvnJg`nGrQ ze!Dv2EbNz0Bsbc5na;0G)CK6HmXwI3sx$O~xp~=S&U&zzS6z3x-WH5+0-kDY3DJMQ z8B%GUt64d=B=nJTNUjzu48LKrz*oFw82l#taRO#>)c}eZz1+pXjO-;Iyi=#xToJ_U zal4Nui4zm}1l}jz(h;N@Z^-Q{j0;0yW!H6t?B4Ph98QH12q@XenF0KTyg59yBQt^E z`kQ2?RVBfMy8gE)6Y2wKEN`6F>sU1Hmea+fd}TsT37G1>P}m4)vcm#<0W^?1!+SAr z5C4vxn(BOO^8#0|H~L_5rI~?=C%; zy}3oM1br2zmPk*@TkZGyaBh6ZiS#XLAmrW6(04$#d0PlCwiK_L&7Q8=Y$-^rc4VZ} zS}&bAarMl9O(}Y1Ma)h8G*pVrY+WA7>>0FW9mgBOoZ5lfxw9b{VU}tM^;5TiPU`ku ze5A5AtXD*>wTFd4GLj6~=sJ<7Qi>&f0$|OO%8DZp*1i(H1y6M#tBvvDFQ`sBz_0ZE zas?e=6qIrw#ckm2FWSUbGw}%N!2`>(10`OwHF1@YNli`*&~YFGBl*g_gwK=}H>~(! z0sJN=hH%Zad+|X!)RP^vYydhtD6?k92v};gnzA!79_h_haO=htW|!^&!QH?b-A4l_ zyEa_VE~yNcYzDg8sFEiXG=jXzvCFE?G%VVLXcu3Mx1Qg0H#%uO);Dh1RueSalOZvc^Rh z?@DWXMkZVby9D{p>}O=Y|0Va#$^ALVnz(9GR@`8r?up zp}q2`ae3<#4MC(AzLeQl0)I{VLb7k|2HOn!#xqo3T1P~$;~!*AEUnlpV0_5jc^A%L zNSkL8(`bheo0hfM!#65ri@V)?{U|rzl{GuT`~Gw%ne^|LCB5%tT+Oud`GrrDtv!Ri z<^0t_zuh#4APs0x)qp2*T>0;tzVKp%XLnBY01E~y!K_I+ZMKcY8GXUtEf9$737n?) zgxs}~hu?Tq6e?Uw8*3I>*Uv4yzkhMN&ZA@Otys0!l+mKVWxS>XkHom}rXnVH)#>kw zk})h&^7z<4d;8cZFdouz@^t}371xo;f7+zj;xc#*aeKfK%6c{nyLblhiHILZP)Wc8 z^^}%?jBEMWPq~74~n^PFdnUP`pxfzhT8y*( z+VnT$;ettRwg0N&4ZJ@!FjODop@)&%<$8LOR??JZk{T9bTN51QV(r)=;ONd{__o=h zD4GAZ()b|uX#r!84k!)=n90;ntndzvwI`3%KJm(F%5^8b`2?kR*O>7LH+1jK z{YA1vc7na$by_VnY>o`U>h`oFPdQvlH3^6jn;x{$)^IdAFJxYT!kt0KciiOvBA;E~ zT$|mr_@O{`eQ>j?oy@gA_m^9#!2-wG-5RDXl-Yh3q)?p#t*Vjx#%ud&&Pmj(?6RmvB}p&s zXsI&VNVjW5##k{mVj#`}zQnh;Ma~4ekcjX^zmT3@%}*W zSS~^o20pd#`ruRj+WP~q&xm+3*sa#IbbC7{aPh!124EA5c1)4EFk6wif{gn5(4&c> zF;7$ZU`A<_X9KaUu)sgVr){&mjH&GS4S*SL@CI-W^^#XpjeaW9tim#mox3hGzX zeZ2iNZLMNMaK)o*n9S%Q(#?7ZlIjQk04ro9H&#~h=>}a8hK*B}NJq-mKl6bOZAVLe zpypcSkzEOS#WD!`V~I#HO`MS_y3MJr|Fj7scSlubBv~Qj;dqdn_8B3Iyt%9GGFRWC z@6BQ$cKCAZG5ki8kgd2lP|zAKX9#7HL(rRcMOz!y70P*Z-DW>I(o2#ZE{{J>PHZ~- z@TDgxmBCYDht1Kap%NXp^m_P+?#Y9VsjY_(1$^Cc_H0#}5k+y8HJk`RZMI)!&~ELq zETOf)aIK0pYCx2x-$<2Jm5Ce~i}UnJDM#tJWmB4xc;vxqx5kOL_2UseW5yD|HpHKA zLt~F0XI<<5{bYZYP*JbH-eWMo{OSAD#djcP4vuLKW$x^givgN0UsSAnT!xM? zi$knQF>;TqdA_#mhgH{QP?zd>usnl(dbt*=6P&wt&9>dXVIp*oy7bD6w#>pHXJ1a6 z`H*eE4*q|!l0&<5+&5%aqQ7F8txc)Wf%cri=e+ppE$RtOK(0}r?$^DiO;U%uDT|I( zH;rV(qRsi&32W!L5<1gthzD!v#r?2Kj~&IDP0~qLHa2v}{S^ritl)5se?UV^{yT*3 z=~E?hU&-BSoK>&~ikp0`yuRh6X~b4OT{zMbVc-D^dLk>e5rq$puAQTJDk$}l+b~VG z<94BAMttQs;=SN%pKeB{!1Wj%xYzc?tH)0m#u;-{|E{~%8*jBGmsM5cRsron$PKmD zWq0ZGGK)u#z~QIU=hw#NNq&BkWx0vg`T*ppy}z{atW(8@Xl>>eR2T0#Zxcg41KE3r z$nxyc!JSgEs+1NPJ!?J3@3k9q;U2}{ikqyLA?G|;tF3bP8N*rW%}zd8YIl5luG2YM zkZMRxG5>mZ3m*4%yDre_Ll}xwL~tD2wC^oA4xKEOR<uS@T#yt{aq75AiTEsv zJlMev35=4{MMkn3&faVMHnh|#n8|F3drH4bC!x}mA#JE`SNqq@UUiHi zM&pBzv_&M8!lM(!Uc2SzGP&zjXlVUJErs0(!f*>c3#x=f zXk3Azm}{YMIbL6U2}|lLp+2kxjPd|{m}U`jzTujYh^2e_$;Xz!qv^@}{2TifK@906 zYuMA+%eC-asgE3~mtl7Q%GLjRk&ZvUXkc;@ybXVV@ArJ#)h3BOA>ZVBq_>;kZD32L zA#KVrnFl~rqSqU^dKwZ)*%l$im>LPGT7+qy-wVr*saDi`3-$@6-&WH)r-l=MFd)G5u9b;*`++CGluK#1LugzM1p$Gc44{@F9benRC&gUjC1zBdrOzyO$qij6?}E z&&=aNH2@2y*FGj5d(}zaFt^KDd}z>r+thRBKKFnV=^zB?_>~vEi>dENdQ~mxn5TYd zt+w`kh~haK>zg^lwTMZ;ULgq>@h;7nmI(bB2&M8_OpAwMiiQ|3;`-(GuV*T9GF=_k z>B_Ki$yzG~1&xZ|71_&lMxT&;BjBI5#np4kxN0Gt1YJp=)416%e~3`0l69_*n4s>- z?H_|57R3qT-4dg#QH8tOC8n$NYW<1i-a+$E%X~6L8GBGD^)3f}dHErwNn)G$q_Z~GjDDe&!1Ec?BOT}%7XdGpC~2AW9C1irVKBEju1XjAxUeaQ zZbX4=s}z;?4qy6fy3Bd;nxON~%SNv4C}5QyC4w}4>VG9U(4ucgCg~Yg&wFMX)Bm~S zPdJv6y`Ob-w<-FziU;T4W|}~kZmgiKK)TM$lc{aOu4sKRx3*UKxlwi4J!GuWQnjOU5R`RjAh0eog8_-!AE5G^K3Z#gtpImHpCtxM&tbwaf9xot$dHVrz0))2ptJ1 zb+cgU*j`y()lIxZj?VM*hha6wdycen;J&>0{PV{%Oz@$T>*V254sDl?hh{ts>Dv-f zOH_fG?Z<_*2O%Qky8W~a=}4dhcT8|7gGmv(y*6Xzm1v@gaP$o#kNep;7W0n8P@g9i?= zp~`TX27@`8x8u5*zeY4z1ZY9$WhGDb9?bpGyR|)T{Qayva0fKT6JfvWLKugy*+9k* zdQ}$3C|jf>{YB&Its|(4Mn`jAjl|mD#vFGDlaahh1YbVUI?%KbyB%rXe_C(QAq`)Xd$(b)aM6MttZv6UJSWg6c=RZts9EDBpD zm$$HUc@#q+yvAq)<>HHzne5(N4_MO|9oc~DjQOIm#|qSOhW`$Phx~niyVN7+huNag zxi;mE<6t)N%YJ|bu)G9FNTSwZa8zD5Bxv=yh6y$9n5Zv|+OL1~I0=N{;sKS*M^gWS8CH>*SoxPj&k7$I+@}j$47mi)luoO$%X<|kFhi|MMUS@8f5#pHK z6HWbF3)h9vWoT@!bpQYZ-0kJ4R%Oz{)CoU?KO)_!l9Cso&MF~;(mFNS3 zTl=CeuW{RPDNmI0YHQ#ePznV8V~7GQWRk<2FKW~LIpQ-NyAANByVN!T_JYkXUP-$p z8|&0Ol@-Q*wA~rrVmM>Q;*hC4V|NlEAT~ic97|~S)FS4JVF&Z+ChnF>H($(@;UnVm%Ao%8jZ zH>uxbN=eV{X>5Olo2%oJ0Sz+P^XxXL|D%V) zfy_gn3lCtF^E?B2S0fWRv>n)bVN`F&|5bcp7veLef+VmKv_P~KXr9&s@=6*vK+~I+ zE7dW;#gefe0(^C@=FbP)tBeTAkC|)tGUX^ve4S~^QdYLWK#~mW6WAIZY&^@0q)6=H z%ao>%`<9$&lus74ga&4cZi%tVi3MccV3HiN;(S_@+kiwtNCk?3r6Jc+QNq)2p8ezG zmkX1+OYIH_lw?hP#8#75fG+#QJxuTVgwz+SuDS@&@VQZFHB*W+M#knLHvg*$DXA0s@cQq zRb*wk!`82OW@=NA*oOG5O*V#ymIwn}gd!!yiaSX~m?f-CUaVY zlJzAs-ZQiRkhgWei!et%T1mQV5t&HWbS^jg6l-X|Q%ht-^1rhHLq8Tq0*^Mtw&zfg zQN_OzbYwF#sf8KPK1QbM5%y_~DcHxUc&fTB&!=W)ZR~;B;DH{pUb>^&8duiw`fz5<-q@>RZTwv?LeV01^gPJVPl84TTv`&EbvZ|K-f~X~|)i~D_<%0Nn zr{lijHi#_^(@BPw*|RdNbiWUtI7zW&4Zzcz?fXrPhlN6nt{lri^oz{zjs# z6XOi?E;k7}Rb&|=%VIgUvjIF9*X4mKp8Nuq03l+tvaMZFoIkg-yZ4L9WaPDbPXjkn zBz6{zQm+GCa3dA;Kv-?GB&8l7&n&&ksveJNuj zIshs{xR(oBC!H11SSamV>zQ4o$^ji+13mG6d3&k@8{zy$spN0>n6&7?En8Nb)yXUi zO*(8(`4}%MeZfC7P2dzkN0Pk{C1f>v(Wex;#a3{>DB9(S*B9ecP~b_K_L4+E8cr zT)cBoksJZpxv8M@{=s-$q$Sbwgd^j|^7B$y$$Nhv>aGLjXYvatZ6l&OiX`>rAB>V) zDX7w>2m2t+FMx; zKIlt(0)MHg>+#`vPtwjS@(yMMy%XDp$o_`x{>(5$9sfP1mJX>aDV#mRHiokS@oK)w z+FB)|@?f9TW(Am8Vl5V^4uz=*e9t$IP1*Tp|MjbE!8_C+OSE6ey4qWK@gi#MmW7|~ z5ZlvMP{m57Tp0}r$XxziS|Zjq%}JhVwuYLrD1xJ~s7RnPg6@8?E5a|KD{fh>Ue!CItrPk=UOUR+Nv5B&OJA5APBBq_SKai@%x_q@z13cNG2z73 z10}QDWV|O!cQ|%LjVTxNdrD6+v^t~XUnJgm=2*3YGzHJAQI}~ir`f1$7H=oWdhJP# zezY;^Hyh+mQ&H=B)o@iEhgS+ZZ6#mp>o5bcJ=WaDkz5P1Qomg8Bzw_Fs=rl3TXMIl zCJM;s88CYd_b0EMKB^wS(O~E_ z`{^}9W{ki;LzZfo-(c1^NM@M-ck!;Z`5d@|eloX?!=9v<#{zi9x?;v_RoKgY`6Lja zO_KTZ5J`NN9^$Zwg4n5C&EH#E8E(6}T&2+f04219AN)Wqgp!Njk+P=Gi^@`4MVlQ-r%)|O|I)f@JnbJD@Slbb9?VO&=-G2Hk+V} z-@gL9zo?==zMOPSn5FR;ZB+#S(s+eKMS`7VL2y>%&d7dCcK;BmRP3ELvk>f9vS!G7$Obj7Kohqn?alA@0 zn-Scn&y!E0)7)3$UuDn-RLIqG)Cas&?B;a~z4Ty|Rgz=h)2kc?{qjsH{G>plXYi&+ zNf|E0;p4=U7G8y?yf+NeFSHmbRK!?}=HL}~uN_cNfsSdiMheCce#^JtO<$}&D6rtr zccl8R@c@-dDZiI#!l`|#_|H4yL+;vO)nynRzc@3MZ@N&!TPH2f{f)A!d)kIIt>LN{ zuj{7yeDdysFgVSfH}W~nO#DpriTB2nh7;^|6VY^Pxnd=_o47W9fuP_~Vv z%3BomCG#gu?e>0#j5W#4%D?HWTW?yf&QV1+?N)Us%s2yaYzzIimUUiQ$PgfLMUw01 zr9MvaFG{6jx*K6B68B`*ucjTU*eaig`dxj`o#`0Rzeigdz*>-f>f(3VO)_EWzNHQ@ zM7U}jf7U$NWSNatK6Bdufm)IaBzIaGLo~r4qWZ8aPI_NrV@;GpyrQ=skVg14YXQ|xQcQui=j9G;D^*@@rIVUmnw9fVkgW^U(218yMitK7eEk2+94(PIbC8s06Rk^=P?t2*dcJ;O2! z+TZk3;{(adFtm5Y(hM+@Q`Q+xgPrALAD?vt>dw;N%I>qRVkWaBN z)AtZYz2g#%I7wOe<>j%E@zgbUJhJLN-Sp*Ut<6hEhev8Tf}3kM&mMOdTjq6(|Gqf= zIld|!eNKMC|9PpKajCWmbMUoveUW<2_|37XWK-FNHPX?GHhcIDrKNFc7|9r94395x z;HV_C3}1dc{%20MsxZ!X8ogW=*Cp6cGN?T4dC-|}^;GC8?@#XDr7w!p$xBOp)&W&V zKdY;g9l8#`dakhcL!)OcG+LL&82(JuXWptAOw}CAyx?<{34X^KQYdWRXT-;lp1^~a z)+*_FRekzmXs;(CbQpV-Z==GU^3K{~+w-2a^;zDrHrKi4?TvNOU^#!~XSBD=Q*Yd= zzS}M6IdQ3zOvnn6Vb#^iZTtDm$f@6yIm&L?n)lzW07vb|k7-juNZw!0wpnMECEd4* zd07cDmp*lbhV{UCJHTp%9nc^Fy@eVEAN z|Hy}wFX$b06hdFpqb|((!fp$bD!SB+74pG*4xQ8CIcI`Um~Ya6{HYYU(Zi6RnPu#f zWfkuq%E$ESv0mV9K^SIfrMY{2F71t~ z=BG2eq}|z+{}0f%G1j#D4jerI8fIl#F6I5ERi0bQo3KO&^((M#5)?KxzT)X_9rX08 zQPY>Y;-^*(HES@^}a4ay%{cu&?L zYBW4FRZf5Cf?Ml!;c<1|@}l6cQeFEt!-|frJ(pt$G$da;(Vc$|`qDfzpJx#rB=~R- zR(0a4m;bJTj5*e$d|^DTg+j+!KU=?bgr-i+&=45YzXJTiyh*>3?=-yjXtirip3}jv zu#EV-8GcD4=#|IXwk!SRFfqfG-{jP``sY9TDk_yyU9a5{Q`opv`?|ARJwXO3r0ZZ+ zf3Wrm3s$R{)xjHm#1kZBc1vKv3ZYQTQ0@uo`D+E=1j!BRDrk-mlZj!$(0jh?HBFvo z%VSnEUQdzDOb|G3Mg)c)0cyazkZ{6siEJU6&8!*~n*Xk-7Uv0c0S&QkFzDH6=sq^- zLSp^l{2(ry3|m*nAI^CO;t1O+FB1Es$GYK7gcArsY(7inm7Xq;7~RnG=k9geVK&rf zYjd&Kl-71c4@|+{EzZHI!IOYm@VC!Y4*%7(mY^@Ko7_|Y;M}>-0VW+ z{fkrY3<_cV*?ja>BS_}`j`mi5+^VWpd7^F4>-|jWwV#cSGz7KJ{gh<)=3AAQZ|0Wv zPK2H_XDz6rEe#8w=CSpkbLdhF(KHf!#H{?fKi&F2{b5R6wF1?w|}Qp-#TG zvS*wwcD$p>J;+dF_(`(@h*wlPfBvYNN#(cr&gUvkjs_X4n$wsYas~a6Xlg*lrtL7c z+@QVL6Iw2+$i!yjoM~tDA#{83)6bF5m4aK+Z|wn|-f9dyp9f9V9ZPt>4%)7Z4=?-O z#w!YSvcNTyzg8+{)pjRZzRF&wYfgU=#m~U-p6A=s30UDil-WF`{PXnd2NrslTYhP% zDtRgJk(+-1?ibbWoQe7;BCRBj>BzcpHL!#`L|rSP$b@sLBrxi+->RxXRwN~O%ypF_ zPv_R$9rkL!k@r<1(T`Yg)X87yUMlUB!?JVXJ*8|Ch|z}?E)BD;4pNmH&oaYm8blzS zObC9M%B!p!>8m?o7+ab&BG}rjsc^D-c|fPI&~unHYrlQ}Jw0^*4NHtHdL%PMat z@;5!oTYJ>`O_yHxRD??4Yj$mdxXPQ#!V^_gua*74e|b1x_>N#Q_EXlWzQ&h2}q2Ku7kN3f%d&Mq4n;64=a6$Dpb{H50?n zNn{Q5klc*CvS|LW%h^=S9=Qj-RJCkhEmYP#U0!pegyXycn4t5BcEJqy+cCi+>Is8& zMeFaT(5$f<@?9x}xHzmzeqZPyY^mkZ>d5?UD7~Z2gyuvn-Z|^E^$Rh+1ZF}6McL;k1T&P-A0GESAwph%rt94Yn<7`hTRJKLxpW*i4iWoY}`e2D4=%()em&(0JZ=~3cuhk^B z+vJrBHudL~#?O!R%o8124dJ*#A;4O8t*n)=cO017!vF#{j!C=PT~eiq;pUkqby#GH7iZr+o3&yp);17J(xZdz>+H&Gzr7l?LUq9 zp&b&6&xf*&OTlx==LFVz58oNVE7JAc-_=HXDQV`e&9?v~my!_9cRwn&A&NJ>4(&>* z1ZkOtq*hrN6>w`&P!$}(&ok1M0@gbYv*7xCb+ z^UhXfKDXFgOa9h0rIm9b5NZ6D%4DfIOeP7LGQqbKLC|&Q^-F05UWgoZrQIMHLE&~n z5j2@+>TqZH(qzXPj?kqLKSe|Bal(|w?d|^qwTtA06=Q0fF`6^H+3& za=&Bwr&Aru5Z4M(HPYDyE-l^l2Lj-|PK{)9N|FC%_+rt{RF){A7IP+PpHhByYzqsYbNstAT*BS;AUNQTH`;1Ltr9x z{zWx%c(h>r;nb6NiX9XEUcz#DrtLYMSugbQj3Ujbog;qtofST%hO}*y2o>w?w}f?t zt49Z@zmU;@nEkb9G`K^RJHqA<`sRlBAN!B+qjuO zc?5kEmNOT2b%_U+9EN#$LLTDrQgk=%j*@D$)}w-}YYQa-3qfc5_5zpg{u6grTR0BM-&5{T_?dJGfjI4|+sDP9hGCPgD9$ zz%Ba%47LrUgD@S`sQD-V;$_5%v9fwXazS!!w#DBhGL=((xqjYEGT-gNDR!^MJN*Qr zRlMxT$6LJ3x2@5m_AacYXVpbRd*NanNuRGdnjSS)%10%9uEA>KiM%O{FJ6H=F)>8D-J-0~CIuPTpXJ%!QyTb6>3BP+%U3}Yi zGfup;s^>xmyCIa5I>T!i%x>*slxY9?rdO(+^m`P>uXckYSVnFxzq-YP`eX6FcK=2# zV>H-rKvHZDmO{vWen0sXEg=H2X2`V-^Tqj8OIJ{ChgXg%S^Mm?U?jLBKwDW^msy2@ z_}a^C&sM+;_bD?-6R?bO2KNdlvidZ?LXlQXU|kI z;n>p4gS4Z)KULw2SRO#lz7$zN`hI9K2&F-$<=*~%5kE>03=RgFAsbrIe@T&?TM2>9 zkSRd;gLQRA942L2m|If!S)`nYJ)77m$GcwG@Y6l?pI}pY_MO| z#J1Y|h1V#TZvQ+*mD=}#>(^zdP=2kv5>Ef+rF0=DsIXR6KwkJu&!FK<#_z&UY8@Lm zb{925I9;4-vvx?5@~#t#pSXQp@H2Ti&kgVa%1^$ElGHHyCUk5QcCCu4dx@Ylwte2) z?<>lBUR?Cf4SeZ!mI?I*zX*QIQ#|(#?T)^%!RXtrzK>_qI`*XQ?u=*DTVG#`j^bkx z8ISzup}uAy)7kVTOXDL!jx*#BEH&z_W&_ehskE4sr^ue(Z+Yx{RmOr@N=orRVHTgzqA9>x!!LjWIXjUoNG-{A##frL5&<1ZV63j+0Uo zy`|hkk_Dzxv2U~{)j>m;}fUy4d_RhoKnF0F2 zzjS)fMHTyIR6;))t|0gD8&?h&sBoHZ`Bz8lyhI=$tH9!jLvd%wnf)+zVH~u9s<(BE zx$$S5`FN<*8NKSkbR-1jU# z_7!V#r$hjn&eS^;W=_N3_=W!YqIx|6vadXyzp2Z^c0b{j^jw@dEL2p6DI6`n`!&<2 zdXqp|Dy-F_P{6G=YgrZD3HXZP*3aR!x9o3^y-}rRw-an+trzhbhb*!ly zl-iH}G}^xoTL0?Ei}=;G5#@8+A`BflGO4EfJOErXtMQY`)wWdo&$*w%pqV9fR!)er2K?8%3A&1ZFwVP;^M&GnAV&2DPc#bcix==^dW2uC_x!8sB|4 zA?0gWvpOl2iic8o>DdVbhS;XetbQ_`_?0iQ((LR@y2 z1C9BJ$igz0Ujdfp?{D^c^@Fy8GDRYy;k=p8hcPqdE3f@u=*z~E96sZ$*ry_RWv(IeYL}|s&O->BhO?dH;PAXd3}*eC*2va!$8hOOgJ=9*QEoeRTgsp z9>cUyoLda@A;`{+`{BRQ``Gp0!eaZFJS)xq7({dKeowje9En{Q#d!W_ZC(G^%lC7D zWAA<-wWg@wQYCuyMn+qB8Ox?oJa=vOC2OoQsetCH9nGqs{Iumf0-I(Na~y|HW7A!*WbFI|k=?1rYX8 zWOhJ+!Z>v2NC`-gZ%O{70qPsqjhMcuIZbWCIy)LvOy6P0BdT5qsn_hiNBTUzwEztP zNK!2J%5vimexF`!4JP#O*1fB%x4vw1?IlKlUg)KAUjlv3`)OBH)$lA>z1;M0W*Wu^ z*>Qo!lO1M^rIzc~N9=Lw+ew8P#9X30BYCuNY)LLH^%mt9Nv>kHa0SnvFW?HZIxdI` zaVQiZD3(Brx0k|6TOk1v8?QLUD~cpT6Yz=0-gz@yW;s=-=wZLsUq@4KbN%x;I-c@o z9@KV%I_MCF^Qd<^=ckc1zd@OVZ>{@gef60{Ec^0_xn{T_Jpya7rH$B~PGmEBMiPu&bxBH7HFe03pImO|achQ_o497yX)Z!&RDOL5k+o;pP`;?G`* z81E)g=Sw&(4$T_@Cv3tE=g7}*;+kX4z7uNNF{trR!S@xiiV!QS@q*9SVxO%peuDX zkb#AJ8BSKp)>}*Rin^T?2%$w8qoIxi60gWaL&SZ z!VuUT_ArJ0nP&-}^k#0gwXAdosN#zWLyYi1CW&n6lig%;=5m(=Z7I>aq$d&(8z!3f z<=cOfr8fV8t=PW7Jp^^zCB?F`h|9jWX=K`e=t@{ZH09Z7!?0g|vPQ(Tpyu@TA+`c1 z%i#OGpAA}{Y++=Y zt2{+i3an(ND7|5b%^%vg=bD{Lvhc{kR&+?8=GbgXZ)D>=pnl7(D+WHw51%pcj@AHZ z7LP1S+qQTR9NTS(rND-8t_PCa<4G3GFqKKi@uT1Ft%Pm+iGHl0#YPC7i3m(7B%Dk0 zh$?&PXEs40h<1m=w-cE~BCh?@GQw`=XZj~1={|5yC7F6DpE|Z(CEz5t_}slVDa;D> z&%@KVo2EnTLJ3}{r+JPY_B`z64bNT}%Pw>zj@bH}9L}~8m6Xb!lfW;{4%e_Yq?b^* zv8*%nWf2uKIp9t(_0GrY(WcP?5VK#{n*BeM$HK;UX0j-1 z`**Fws8%e`!Tpxc6Ysxo2&?qySM`F2E*@-)iYe+T2VP2p;1qPM&t~)tbGu7M<%C{N z=(0K<^*}zzQFVPs`rRspeuyzeU9ry2%lu#z{LN0OXYCWn z(#imT7o@?yS10(hlr3{@BCZDayy>{6jqb_EURc}R#-A+l7v}3vv)AoAKWPAkx*e(9 zDm_qVsh!>0y1$<(HN@6xBT-H@&;>tsV`4APf?21eOPaQH01EZ&bgS6n6Ni|w$?e@@ z$`*}HsrYN-C-$Gid-)IT5}tBj2{w?}m|sYdaqvp@>!FvABU`EFy8$kHwN)F)2-m z{B5bpMV*WM1>xKR=0J@y1Y5*g_ZW)V50v2_eIb|h=&SETmP4$japgmmV@q2=vV=?X zp+8DWv1;Jkbu(8`#sdluiiZ(_Pe1}H0aXjSt5MBv+f$|Wjy+iNKHLs&{`O~Du{Wn|Uk4Lq zLFR-n&(50Z9i?-}!8({c*riT0%FR&ZPAl;xI3vHBU6{MPS-&d@-yBWxIBH<6ZUHRj zTLizGM2|7SGtoAVcG=o17}@GC(N}bzbesu6$Dm|n!;rVd_~<9s38`pf``2ix#+oI- zvT^6~AFafx7Gcvo?Vvlz@9j9^xBbryUItCPZ|)guL~f?bb=(eIJ{q&sOlYp`y%XQ; zwE3oa3X|f+@&v)#RLicBUcD)m^*57Y|HWiMXScLgx%Z!rjB^NN>vOJqT|Ps3i7eY{ zhe9~^zBJ`5|B8z@Ol3l9W?Pv~py#WNWFtTN@E!|~n#0ZSAf+E2`t>$^UH|2?W=|7Z z2gYUmNhHLx2!qjc48FJ*WEZ^RpEL!4BDl{8BX(n4kI*z-hZ zK_1U#ql;!*j@7d*4=~g2eUKNQ8)P__lC;VPKIgQJ6LTDx;mW2H3x4HA?H>c>ct2_d ztr}W_Kx3R~NP%U?Am|gae*nQ7nfvov*W%34SFXiC^8g=mF$@o#hJA^)lWF>hs4G^&4i0x7xNL;Bq20(f_3IW5WD2lQe{W zPq$#Q02TKFBXlBQ=)-B&aez;9*!OE>l+EIF#znTk#Tt0KN4}NEng>b>$MGjm$608i zFW<6IpsHUH(LvSM+gh-$%GIn)*0$tZek8PJDYI{u_?$qU#qKei0EccJ3Id%Fi0-EC z@4H%vfEi^CtNbylaC35etO9j3Z?1vK)bAqu>$afGx~Le1g=%?VizJl-9r?j2d)YDN+_J1xIFG z5Tvww4eXZ!Y`m;m>m0RjBPjlK~}R?D`+ zGBye(-m~c!K|f3cx;A)r7sk2SFY$gg$L5dy%u@o4Z8XETe+6$vfE=bYpO^Fq17eTG zG(3O)Cz8fi?$C+x{@S?KUJ93PeLSd;bC|Ep`>OwC0^e3vw80{a0NT{r+E=~CJOH3_ z9cO09TA-lAi%K*~P1cuPFwQ8-w;VQm|M=ocef?QZI8fD1SyQlHqE5UOVU(mau|+Al zX-j8`nYAp>K=X(-QuXoxmJ5=IWU5Q!_sWQE#VZZ#;E$SnkKuPpA+K|iRf~IAKTLdC zeJ&QsMUtWl_dJzoi?dTz6`9gAgYSa_PtWQd6g#QLtoG4bxl%o#6$8QmLR0S95SE{1 zWaxO}f}bB&@sS6FH-%(TL457go4ilv{%x^i7U_79H^u%tCEpTQn&WdHls3Z&Scv^D z_*?Z6!@js8Q>luu_~_WK_S*Bbm$oBBHLgDln?Mh0bwe~&q9Tyd8cgK~3wGF-i1cLS zc(zrUOgAmEmadH3hJq`q6(^^#zm83FG7s{3AR_ej%T?_>yNYemJUa~GtG`NPfX^ny zQw%G#i$t}URxB7mE5=NOjZ;Rbsl#%t)XMRa1H#Lga{=lHNhwF3FKY(RSlYzPdXIRM zjMQUE)e>mx87!&)&b#02dh9CI7@H!cB}V0l{;%)s9H$x$CY0iri?#k%Qn(e_H;A=@ zCF>g34pe+}lS8{tPMPs})sqcF8#ST23yV~L4e&$S=zV=S-|?#)cdtZQAFB3Hv+p`= zl9nK={~ye%HN|{vtD>or-5my)cCBwxUtj4z{A$I zI(;=5)1gUqwMv7Bs;%5W$KG{E?)as%v)AUSw|0~eX@b9%*oMx5g(eEDJaDqrRz4qO zy{D_K0``x7zo#M?vShP7S(#s@IW3LWZ;FpG$4RdDMe=^I@nUM2AHaqPKUR5}ijc20 zzo#8)3wSg@2EOm7INgZYY5R5ZZ1vq#zv?g>FmW0)w_W=Di)zH3jnEj2w}01-m&fgo zT4hT*S7hGujADp?BD$A{N==oEDLv)!Hzzcmm)8AP-TfEc4na4slQ$2a5QKp{Q(qp} z5lu0nQzOrAXDJz8b6yfiBWFtOFZ6jSfcN&3sbHbDe-%|pJDvVJo%xnI zZDm$Ro@@rdwHR6IGZ4uyOEpE}MdXVfV{Mo2rOn@vh`(_kwl~hbXm!M;i4P_e&eBm$ z(VHY;pu=?1abShzqw>8uT=^Tb-2ySmKl)7>z_>!mH^guPp?z z(VObDAs*3MEt_)N86+!>XyxKHo}k7Krmrf7F9~2vID;M%y;>>{ zU+WxOnYX)Yv(O4Jd&gk!WXc3q)FTcC>hcE6{efA%9GEseiPBC5KwzIh- zFr<)^D(`=Wy47guktyt4(j3P^X7QVAhkDw0j(N#^VB7an_MK6X=}lufe1zq}oE?p~ zema8=F+>!L1`Rwr{&;Z)5sW%>V|R_{%JaR9JD5~GNf~+R>T#UPB30eNBhX?+x&)kPM6Y2ftc`GHh^Ub;jpZ;EQe|qcf;9NtGbB=`9er? zGt?Qqd(+W8n1c+jP~B`yr%pJQinfG=;yQC?MT1sDQ{-@+5QELgqe@;3Q1!EK`-bj8 zQ1vaOR7NS_*@kFPG~`}cPwTzUiKFpSQb!cq85dAMd;KGa*Rc)sNl?p=6Tddg_VN^E zM1IIQCW6)Zhf$p;%Xtwa!p$ozA~QO#th2r)U%2*vkv&EIc&A@;we@tZHFxal&E3qP zIO_fm>Nuhu&9jfE)Q9{DrS*(b&GH{&R;2%6y!8}zdbpgoJ45~Gs^;muX71RQ?SH9T z5$h#cJ!#Q@<2~D8-w*%uQO}D0FQI(Z`JbIID)E0FVS3*Gji6wZ`G*z9?!SNHkiW6> ee?u0IHH$u|XFjROBmQfe)ngr#N7dR+(f<#PzFqwQ literal 0 HcmV?d00001 From 3c46f4dbefa96ac862ddd8806886b0a3b182931b Mon Sep 17 00:00:00 2001 From: Brendan Collins Date: Mon, 6 Apr 2026 05:34:42 -0700 Subject: [PATCH 7/7] Add KDE section to README feature matrix (#1143) --- README.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/README.md b/README.md index 0d0099b4..c04b640f 100644 --- a/README.md +++ b/README.md @@ -517,6 +517,15 @@ Same-CRS tiles skip reprojection entirely and are placed by direct coordinate al -------- +### **Kernel Density Estimation** + +| Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray | +|:-----|:------------|:------:|:------------------:|:-----------------:|:---------------------:|:---------------------:| +| [KDE](xrspatial/kde.py) | Point-to-raster kernel density estimation (Gaussian, Epanechnikov, quartic) | Silverman 1986 | ✅️ | ✅️ | ✅️ | ✅️ | +| [Line Density](xrspatial/kde.py) | Line-segment-to-raster density estimation | Standard | ✅️ | | | | + +-------- + ### **Multivariate** | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |