From b6753569a1a68b499a2a86dae5eb60deed9eb586 Mon Sep 17 00:00:00 2001
From: Panagiotis Tomer Karagiannis
Date: Mon, 19 Dec 2022 22:39:25 +0100
Subject: [PATCH] wip
---
.gitignore | 2 +
drawing_to_fsd_layout/common.py | 5 +
drawing_to_fsd_layout/cone_placement.py | 301 +++++++
drawing_to_fsd_layout/export.py | 250 ++++++
drawing_to_fsd_layout/image_processing.py | 244 ++++++
drawing_to_fsd_layout/math_utils.py | 638 +++++++++++++++
drawing_to_fsd_layout/spline_fit.py | 135 ++++
drawing_to_track.ipynb | 924 ++++++++++++++++++++++
requirements.txt | 6 +
streamlit_app.py | 328 ++++++++
10 files changed, 2833 insertions(+)
create mode 100644 drawing_to_fsd_layout/common.py
create mode 100644 drawing_to_fsd_layout/cone_placement.py
create mode 100644 drawing_to_fsd_layout/export.py
create mode 100644 drawing_to_fsd_layout/image_processing.py
create mode 100644 drawing_to_fsd_layout/math_utils.py
create mode 100644 drawing_to_fsd_layout/spline_fit.py
create mode 100755 drawing_to_track.ipynb
create mode 100644 requirements.txt
create mode 100644 streamlit_app.py
diff --git a/.gitignore b/.gitignore
index b6e4761..f347387 100644
--- a/.gitignore
+++ b/.gitignore
@@ -127,3 +127,5 @@ dmypy.json
# Pyre type checker
.pyre/
+
+.vscode/
\ No newline at end of file
diff --git a/drawing_to_fsd_layout/common.py b/drawing_to_fsd_layout/common.py
new file mode 100644
index 0000000..73b4614
--- /dev/null
+++ b/drawing_to_fsd_layout/common.py
@@ -0,0 +1,5 @@
+import numpy as np
+
+FloatArrayNx2 = np.typing.NDArray[np.float64]
+IntArrayN = np.typing.NDArray[np.int64]
+FloatArray2 = np.typing.NDArray[np.float64]
diff --git a/drawing_to_fsd_layout/cone_placement.py b/drawing_to_fsd_layout/cone_placement.py
new file mode 100644
index 0000000..46a2248
--- /dev/null
+++ b/drawing_to_fsd_layout/cone_placement.py
@@ -0,0 +1,301 @@
+import networkx as nx
+import numpy as np
+import streamlit as st
+from scipy.ndimage import uniform_filter1d
+from scipy.spatial.distance import cdist
+
+from drawing_to_fsd_layout.common import FloatArray2, FloatArrayNx2, IntArrayN
+
+
+def circle_fit(coords: np.ndarray, max_iter: int = 99) -> np.ndarray:
+ """
+ Function taken from: https://github.com/papalotis/ft-fsd-path-planning/blob/d82ba8f93c753a9d0fe0c77fa8c9af88aafad6ea/fsd_path_planning/utils/math_utils.py#L569
+
+ Fit a circle to a set of points. This function is adapted from the hyper_fit
+ function in the circle-fit package (https://pypi.org/project/circle-fit/).
+ The function is a njit version of the original function with some input validation
+ removed. Furthermore, the residuals are not calculated or returned.
+ Args:
+ coords: The coordinates of the points as an [N, 2] array.
+ max_iter: The maximum number of iterations.
+ Returns:
+ An array with 3 elements:
+ - center x
+ - center y
+ - radius
+ """
+
+ X = coords[:, 0]
+ Y = coords[:, 1]
+
+ n = X.shape[0]
+
+ Xi = X - X.mean()
+ Yi = Y - Y.mean()
+ Zi = Xi * Xi + Yi * Yi
+
+ # compute moments
+ Mxy = (Xi * Yi).sum() / n
+ Mxx = (Xi * Xi).sum() / n
+ Myy = (Yi * Yi).sum() / n
+ Mxz = (Xi * Zi).sum() / n
+ Myz = (Yi * Zi).sum() / n
+ Mzz = (Zi * Zi).sum() / n
+
+ # computing the coefficients of characteristic polynomial
+ Mz = Mxx + Myy
+ Cov_xy = Mxx * Myy - Mxy * Mxy
+ Var_z = Mzz - Mz * Mz
+
+ A2 = 4 * Cov_xy - 3 * Mz * Mz - Mzz
+ A1 = Var_z * Mz + 4.0 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz
+ A0 = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy
+ A22 = A2 + A2
+
+ # finding the root of the characteristic polynomial
+ y = A0
+ x = 0.0
+ for _ in range(max_iter):
+ Dy = A1 + x * (A22 + 16.0 * x * x)
+ x_new = x - y / Dy
+ if x_new == x or not np.isfinite(x_new):
+ break
+ y_new = A0 + x_new * (A1 + x_new * (A2 + 4.0 * x_new * x_new))
+ if abs(y_new) >= abs(y):
+ break
+ x, y = x_new, y_new
+
+ det = x * x - x * Mz + Cov_xy
+ X_center = (Mxz * (Myy - x) - Myz * Mxy) / det / 2.0
+ Y_center = (Myz * (Mxx - x) - Mxz * Mxy) / det / 2.0
+
+ x = X_center + X.mean()
+ y = Y_center + Y.mean()
+ r = np.sqrt(abs(X_center**2 + Y_center**2 + Mz))
+
+ return np.array([x, y, r])
+
+
+def create_cyclic_sliding_window_indices(
+ window_size: int, step_size: int, signal_length: int
+) -> IntArrayN:
+ """
+ Function taken from https://github.com/papalotis/ft-fsd-path-planning/blob/main/fsd_path_planning/calculate_path/path_parameterization.py
+ """
+ if window_size % 2 == 0:
+ raise ValueError("Window size must be odd.")
+ half_window_size = window_size // 2
+
+ indexer = (
+ np.arange(-half_window_size, half_window_size + 1)
+ + np.arange(0, signal_length, step_size).reshape(-1, 1)
+ ) % signal_length
+ return indexer
+
+
+def calculate_path_curvature(
+ path: FloatArrayNx2, window_size: int, path_is_closed: bool
+) -> FloatArrayNx2:
+ """
+ Function taken from https://github.com/papalotis/ft-fsd-path-planning/blob/main/fsd_path_planning/calculate_path/path_parameterization.py
+
+ Calculate the curvature of the path.
+ Args:
+ path: The path as a 2D array of points.
+ window_size: The size of the window to use for the curvature calculation.
+ path_is_closed: Whether the path is closed or not.
+ Returns:
+ The curvature of the path.
+ """
+ windows = create_cyclic_sliding_window_indices(
+ window_size=window_size, step_size=1, signal_length=len(path)
+ )
+
+ path_curvature = np.zeros(len(path))
+ for i, window in enumerate(windows):
+ if not path_is_closed:
+ diff = window[1:] - window[:-1]
+ if np.any(diff != 1):
+ idx_cutoff = int(np.argmax(diff != 1) + 1)
+ if i < window_size:
+ window = window[idx_cutoff:]
+ else:
+ window = window[:idx_cutoff]
+
+ points_in_window = path[window]
+
+ _, _, radius = circle_fit(points_in_window)
+ radius = min(
+ max(radius, 1.0), 3000.0
+ ) # np.clip didn't work for some reason (numba bug?)
+ curvature = 1 / radius
+ three_idxs = np.array([0, int(len(points_in_window) / 2), -1])
+ three_points = points_in_window[three_idxs]
+ hom_points = np.column_stack((np.ones(3), three_points))
+ sign = np.linalg.det(hom_points)
+
+ path_curvature[i] = curvature * np.sign(sign)
+
+ mode = "wrap" if path_is_closed else "nearest"
+
+ filtered_curvature = uniform_filter1d(
+ path_curvature, size=window_size // 10, mode=mode
+ )
+ # filtered_curvature = path_curvature
+
+ return filtered_curvature
+
+
+def calculate_edge_curvature(edge: FloatArrayNx2) -> FloatArrayNx2:
+ window_size = len(edge) // 25
+ if window_size % 2 == 0:
+ window_size += 1
+ return calculate_path_curvature(edge, window_size=window_size, path_is_closed=True)
+
+
+def place_cones_for_edge(edge: FloatArrayNx2) -> FloatArrayNx2:
+ """
+ Place cones along an edge.
+
+ Args:
+ edge: The edge to place cones on. Shape (n,2)
+ cone_spacing: The distance between cones.
+
+ Returns:
+ The locations of the cones. Shape (m,2)
+ """
+ raise NotImplementedError
+
+
+def calculate_min_track_width_index(
+ edge_1: FloatArrayNx2, edge_2: FloatArrayNx2
+) -> tuple[int, int]:
+ distances = cdist(edge_1, edge_2)
+ idx_min_distance_to_1 = np.min(distances, axis=1).argmin()
+ idx_min_distance_to_2 = np.min(distances, axis=0).argmin()
+
+ return idx_min_distance_to_1, idx_min_distance_to_2
+
+
+def calculate_min_track_width(edge_1: FloatArrayNx2, edge_2: FloatArrayNx2) -> float:
+ idx_min_distance_to_1, idx_min_distance_to_2 = calculate_min_track_width_index(
+ edge_1, edge_2
+ )
+ edge_1_closest_point = edge_1[idx_min_distance_to_1]
+ edge_2_closest_point = edge_2[idx_min_distance_to_2]
+
+ return np.linalg.norm(edge_1_closest_point - edge_2_closest_point)
+
+
+def estimate_centerline_from_edges(
+ edge_1: FloatArrayNx2, edge_2: FloatArrayNx2
+) -> FloatArrayNx2:
+ shorter, longer = sorted([edge_1, edge_2], key=len)
+ idx_shorter = cdist(shorter, longer).argmin(axis=0)
+ shorter_sampled = shorter[idx_shorter]
+ centerline = (longer + shorter_sampled) / 2
+
+ return centerline
+
+
+def estimate_centerline_length(edge_1: FloatArrayNx2, edge_2: FloatArrayNx2) -> float:
+ edge_1_sampled = edge_1[::10]
+ edge_2_sampled = edge_2[::10]
+ centerline = estimate_centerline_from_edges(edge_1_sampled, edge_2_sampled)
+ return np.sum(np.linalg.norm(np.diff(centerline, axis=0), axis=1))
+
+
+def split_edge_to_straight_and_curve(
+ edge: FloatArrayNx2,
+ curvature_threshold_upper: float = 0.02,
+ curvature_threshold_lower: float = 0.01,
+) -> IntArrayN:
+ """
+ Split an edge into straight and curve sections.
+ """
+ curvature_threshold_lower, curvature_threshold_upper = sorted(
+ (curvature_threshold_lower, curvature_threshold_upper)
+ )
+ curvature = calculate_edge_curvature(edge)
+ start_index = np.argmin(curvature) # start from a straight
+
+ # 0 is straight, 1 is left curve, -1 is right curve
+ out = np.zeros(len(curvature), dtype=int)
+ out[start_index] = 0
+ for offset in range(1, len(curvature) + 1):
+ previous_index = (start_index + offset - 1) % len(curvature)
+ next_index = (start_index + offset) % len(curvature)
+
+ state = out[previous_index]
+ if state == 0 and abs(curvature[next_index]) > curvature_threshold_upper:
+ state = np.sign(curvature[next_index])
+ elif state != 0 and abs(curvature[next_index]) < curvature_threshold_lower:
+ state = 0
+
+ out[next_index] = state
+
+ return out
+
+
+def decide_start_finish_line_position_and_direction(
+ edge: FloatArrayNx2,
+) -> tuple[FloatArray2, FloatArray2]:
+ track_point_types = split_edge_to_straight_and_curve(edge)
+ straights_indices: list[list[int]] = [[]]
+ for i, track_point_type in enumerate(track_point_types):
+ if track_point_type == 0:
+ straights_indices[-1].append(i)
+ else:
+ if len(straights_indices[-1]) > 0:
+ straights_indices.append([])
+
+ longest_straight_indices = max(straights_indices, key=len)
+ start_index = longest_straight_indices[len(longest_straight_indices) // 2]
+ direction_of_straight = np.diff(edge[longest_straight_indices], axis=0).mean(axis=0)
+ return edge[start_index], direction_of_straight / np.linalg.norm(
+ direction_of_straight
+ )
+
+
+def fix_edge_direction(
+ edge: FloatArrayNx2, start_position: FloatArray2, start_direction: FloatArray2
+) -> FloatArrayNx2:
+ idx_edge_start = np.linalg.norm(edge - start_position, axis=1).argmin()
+ edge_rolled = np.roll(edge, -idx_edge_start, axis=0)
+
+ edge_direction = edge_rolled[2] - edge_rolled[0]
+ dot = np.dot(edge_direction, start_direction)
+ if dot < 0:
+ edge_rolled = edge_rolled[::-1]
+
+ return edge_rolled
+
+
+def place_cones(
+ trace: FloatArrayNx2, seed: int, mean: float, variance: float
+) -> FloatArrayNx2:
+ rng = np.random.default_rng(seed)
+
+ idxs_to_keep = [0]
+ next_distance = rng.normal(mean, variance)
+ next_idx = idxs_to_keep[0]
+ while next_idx < len(trace):
+ trace_from_last_in = trace[next_idx:]
+ dist_to_next = np.linalg.norm(np.diff(trace_from_last_in, axis=0), axis=1)
+ cum_dist = np.cumsum(dist_to_next)
+
+ offset_next_idx = np.argmax(cum_dist > next_distance)
+
+ if offset_next_idx == 0:
+ break
+
+ next_idx = offset_next_idx + next_idx
+ idxs_to_keep.append(next_idx)
+ next_distance = rng.normal(mean, variance)
+
+ randomly_placed_cones = trace[idxs_to_keep]
+
+ st.write(np.linalg.norm(np.diff(randomly_placed_cones, axis=0), axis=1).mean())
+
+
+ return randomly_placed_cones
diff --git a/drawing_to_fsd_layout/export.py b/drawing_to_fsd_layout/export.py
new file mode 100644
index 0000000..2f035c5
--- /dev/null
+++ b/drawing_to_fsd_layout/export.py
@@ -0,0 +1,250 @@
+import json
+from enum import Enum
+from struct import Struct
+from typing import Dict, List, Literal, Tuple, cast
+
+import numpy as np
+
+from drawing_to_fsd_layout.common import FloatArrayNx2
+
+
+def export_json_string(edges_left: FloatArrayNx2, edges_right: FloatArrayNx2) -> str:
+ """
+ Export the track edges to a JSON string that can be used in the FSD layout editor.
+
+ Args:
+ edges_left: The left track edge. Shape (n,2)
+ edges_right: The right track edge. Shape (n,2)
+
+ Returns:
+ The JSON string
+ """
+ cones_x = np.concatenate((edges_left[:, 0], edges_right[:, 0])).tolist()
+ cones_y = np.concatenate((edges_left[:, 1], edges_right[:, 1])).tolist()
+ cones_color = (
+ ["orange_big"]
+ + ["blue"] * (len(edges_left) - 1)
+ + ["orange_big"]
+ + ["yellow"] * (len(edges_right) - 1)
+ )
+ obj = {
+ "x": cones_x,
+ "y": cones_y,
+ "color": cones_color,
+ }
+
+ return json.dumps(obj)
+
+
+# --- Export to LYT file ---
+
+
+def angle_from_2d_vector(vecs: np.ndarray) -> np.ndarray:
+ if vecs.shape == (2,):
+ return np.arctan2(vecs[1], vecs[0])
+ if vecs.ndim == 2 and vecs.shape[-1] == 2:
+ return np.arctan2(vecs[:, 1], vecs[:, 0])
+ raise ValueError("vecs can either be a 2d vector or an array of 2d vectors")
+
+
+class ConeTypes(int, Enum):
+ """
+ Enum for all possible cone types
+ """
+
+ UNKNOWN = 0
+ YELLOW = RIGHT = 1
+ BLUE = LEFT = 2
+ ORANGE_SMALL = START_FINISH_AREA = 3
+ ORANGE_BIG = START_FINISH_LINE = 4
+
+
+HEADER_STRUCT = Struct("6sBBhBB")
+BLOCK_STRUCT = Struct("2h4B")
+
+
+ConeTypeToLytObjectIndex: Dict[ConeTypes, int] = {
+ ConeTypes.UNKNOWN: 25,
+ ConeTypes.YELLOW: 29, # 30 also possible
+ ConeTypes.BLUE: 23, # 24 also possible
+ ConeTypes.ORANGE_BIG: 27,
+ ConeTypes.ORANGE_SMALL: 20, # 20 is red, we use it to represent a small orange cone
+}
+
+LytObjectIndexToConeType: Dict[int, ConeTypes] = {
+ 25: ConeTypes.UNKNOWN,
+ 29: ConeTypes.YELLOW,
+ 30: ConeTypes.YELLOW,
+ 23: ConeTypes.BLUE,
+ 24: ConeTypes.BLUE,
+ 27: ConeTypes.ORANGE_BIG,
+ 20: ConeTypes.ORANGE_SMALL,
+}
+
+
+def _to_lyt_heading(heading: np.ndarray, input_is_radians: bool) -> np.ndarray:
+ """
+ Convert real world heading (direction) to lyt heading
+
+ Args:
+ heading (np.ndarray): The heading in the real world
+ input_is_radians (bool, optional): If set to True the input will be converted to
+ degrees as required by LYT.
+
+ Returns:
+ np.ndarray: The heading is required by LYT
+ """
+ if input_is_radians:
+ heading = np.rad2deg(heading)
+ return ((heading + 180) * 256 // 360).astype(np.uint8)
+
+
+def _create_lyt_block_for_cone(
+ x_pos: int, y_pos: int, heading: int, color_index: int
+) -> bytes:
+ z_height = 240 # suggested by documentation (puts element on ground)
+ flags = 0 # these are simple cone objects no flags
+ block = BLOCK_STRUCT.pack(x_pos, y_pos, z_height, flags, color_index, heading)
+ return block
+
+
+def _create_lyt_trace_bytes(trace: np.ndarray, color_idx: int) -> bytes:
+ if len(trace) > 0:
+ trace_looped = cast(np.ndarray, np.vstack((trace, trace[:1])))
+ heading = angle_from_2d_vector(trace_looped[1:] - trace_looped[:-1])
+ # NOTE 2 in format
+ heading = _to_lyt_heading(heading, input_is_radians=True)
+ all_blocks = b"".join(
+ _create_lyt_block_for_cone(x, y, h, color_idx)
+ for (x, y), h in zip(trace, heading)
+ )
+ return all_blocks
+ return b""
+
+
+def _create_start_block(x_pos: int, y_pos: int, heading: float) -> bytes:
+ z_height = 240
+ index = 0
+ flags = 0 # start position has 0 width
+ heading = _to_lyt_heading(heading, input_is_radians=True)
+ block = BLOCK_STRUCT.pack(x_pos, y_pos, z_height, flags, index, heading)
+ return block
+
+
+def _create_finish_block(x_pos: int, y_pos: int, heading: float, width: float) -> bytes:
+ block = bytearray(_create_start_block(x_pos, y_pos, heading))
+ # width of finish object
+ block[5] |= int(width / 2) << 2
+ return bytes(block)
+
+
+def _create_checkpoint_block(
+ x_pos: int,
+ y_pos: int,
+ heading: float,
+ width: float,
+ checkpoint_index: Literal[1, 2, 3],
+) -> bytes:
+ if checkpoint_index not in (1, 2, 3):
+ raise ValueError(
+ f"checkout_index must be either 1, 2 or 3. It is {checkpoint_index}"
+ )
+ block = bytearray(_create_finish_block(x_pos, y_pos, heading, width))
+ block[5] |= checkpoint_index
+ return bytes(block)
+
+
+def _traces_to_lyt_bytes(
+ cones_per_type: List[np.ndarray], offset_in_meters: Tuple[float, float]
+) -> bytes:
+ lfs_scale = 16
+ offset = offset_in_meters * lfs_scale
+
+ cones_in_map = [(c * lfs_scale + offset).astype(int) for c in cones_per_type]
+
+ cones_bytes = [
+ _create_lyt_trace_bytes(cones, ConeTypeToLytObjectIndex[cone_type])
+ for cone_type, cones in zip(ConeTypes, cones_in_map)
+ ]
+
+ right_in_map = cones_in_map[ConeTypes.RIGHT]
+ left_in_map = cones_in_map[ConeTypes.LEFT]
+ start_finish_in_map = cones_in_map[ConeTypes.START_FINISH_LINE]
+
+ pos_x, pos_y = (left_in_map[-2] + right_in_map[-2]) // 2
+ start_heading: float = (
+ angle_from_2d_vector(left_in_map[-2] - left_in_map[-1]) + np.pi / 2
+ )
+ start_block = _create_start_block(pos_x, pos_y, start_heading)
+
+ finish_pos_x, finish_pos_y = (start_finish_in_map[0] + start_finish_in_map[1]) // 2
+
+ # divide by lfs_scale because we need actual width
+ finish_width = 3 * (
+ np.linalg.norm(start_finish_in_map[0] - start_finish_in_map[1]) / lfs_scale
+ )
+ finish_heading = angle_from_2d_vector(left_in_map[-1] - left_in_map[0]) + np.pi / 2
+ finish_block = _create_finish_block(
+ finish_pos_x, finish_pos_y, finish_heading, finish_width
+ )
+
+ # put checkpoint approximately in the middle
+ # this is so that lfs counts lap times
+ half_len = len(left_in_map) // 2
+ left_point_half = left_in_map[half_len]
+ right_point_half_index = np.linalg.norm(
+ left_point_half - right_in_map, axis=1
+ ).argmin(axis=0)
+
+ right_point_half = right_in_map[right_point_half_index]
+ check_pos_x, check_pos_y = (left_point_half + right_point_half) // 2
+ check_width = 5 * (np.linalg.norm(left_point_half - right_point_half) / lfs_scale)
+ check_heading = (
+ angle_from_2d_vector(left_in_map[half_len - 1] - left_in_map[half_len])
+ + np.pi / 2
+ )
+
+ print(check_heading)
+ check_block = _create_checkpoint_block(
+ check_pos_x, check_pos_y, check_heading, check_width, 1
+ )
+
+ final_object_blocks = b"".join(
+ (start_block, finish_block, check_block, *cones_bytes)
+ )
+
+ n_obj = len(final_object_blocks) // BLOCK_STRUCT.size
+ assert len(final_object_blocks) % BLOCK_STRUCT.size == 0
+ header = HEADER_STRUCT.pack(b"LFSLYT", 0, 251, n_obj, 10, 8)
+
+ final_bytes = b"".join((header, final_object_blocks))
+ return final_bytes
+
+
+def cones_to_lyt(
+ world_name: Literal["BL4", "AU1", "AU2", "AU3", "WE3", "LA2"],
+ cones_left: FloatArrayNx2,
+ cones_right: FloatArrayNx2,
+) -> bytes:
+ offset = {
+ "BL4": (-261, 124),
+ "AU1": (-50, -1010),
+ "AU2": (-138, -696),
+ "AU3": (-66, -50),
+ "WE3": (64, -1200),
+ "LA2": (538, 548),
+ }[world_name]
+
+ offset = np.array(offset)
+
+ assert offset is not None
+
+ cones_per_type = [np.empty((0, 2)) for _ in ConeTypes]
+ cones_per_type[ConeTypes.LEFT] = cones_left[1:]
+ cones_per_type[ConeTypes.RIGHT] = cones_right[1:]
+ cones_per_type[ConeTypes.START_FINISH_LINE] = np.row_stack(
+ (cones_left[0], cones_right[0])
+ )
+
+ bytes_to_write = _traces_to_lyt_bytes(cones_per_type, offset)
+ return bytes_to_write
diff --git a/drawing_to_fsd_layout/image_processing.py b/drawing_to_fsd_layout/image_processing.py
new file mode 100644
index 0000000..8aa2e75
--- /dev/null
+++ b/drawing_to_fsd_layout/image_processing.py
@@ -0,0 +1,244 @@
+from pathlib import Path
+from typing import Iterable
+
+import matplotlib.pyplot as plt
+import networkx as nx
+import numpy as np
+import streamlit as st
+from scipy.spatial.distance import cdist
+from skimage import io
+from skimage.color import rgb2gray
+from skimage.feature import canny
+from skimage.filters import unsharp_mask
+from skimage.transform import rescale
+
+from drawing_to_fsd_layout.common import FloatArrayNx2
+
+Image = np.typing.NDArray[np.float64]
+
+
+def rotate(points: np.ndarray, theta: float) -> np.ndarray:
+ """
+ Rotates the points in `points` by angle `theta` around the origin
+
+ Args:
+ points: The points to rotate. Shape (n,2)
+ theta: The angle by which to rotate in radians
+
+ Returns:
+ The points rotated
+ """
+ cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+ rotation_matrix = np.array(((cos_theta, -sin_theta), (sin_theta, cos_theta))).T
+ return np.dot(points, rotation_matrix)
+
+
+def load_raw_image(image_path: Path | str) -> Image:
+ """
+ Load an image. If input is a string, it is assumed to be a path to an image. If
+ input is a Path, it is assumed to be a path to an image. If input is an image, it is
+ returned as is.
+ """
+ if isinstance(image_path, str):
+ image_path = Path(image_path)
+
+ image = io.imread(image_path) if isinstance(image_path, Path) else image_path.copy()
+
+ return image
+
+
+def create_close_point_adjacency_matrix(
+ edge_1: FloatArrayNx2, edge_2: FloatArrayNx2, distance: float
+) -> np.ndarray:
+ """
+ Combine two edges into a single edge.
+
+ Args:
+ edge_1: The first edge. Shape (n,2)
+ edge_2: The second edge. Shape (m,2)
+
+ Returns:
+ The combined edge.
+ """
+ dist = cdist(edge_1, edge_2)
+
+ # find the closest point in edge_2 for each point in edge_1
+ idxs_closest_in_edge_2 = dist.argmin(axis=1)
+
+ # find the closest point in edge_1 for each point in edge_2
+ idxs_closest_in_edge_1 = dist.argmin(axis=0)
+
+ adj_edge_1_to_edge_2 = idxs_closest_in_edge_1 < distance
+ adj_edge_2_to_edge_1 = idxs_closest_in_edge_2 < distance
+
+ combined_length = len(edge_1) + len(edge_2)
+
+ adj = np.zeros((combined_length, combined_length), dtype=bool)
+
+ adj[: len(edge_1), len(edge_1) :] = adj_edge_1_to_edge_2
+ adj[len(edge_1) :, : len(edge_1)] = adj_edge_2_to_edge_1
+
+ return adj
+
+
+def combine_edges(edge_1: FloatArrayNx2, edge_2: FloatArrayNx2) -> FloatArrayNx2:
+ """
+ Combine two edges into a single edge.
+
+ Args:
+ edge_1: The first edge. Shape (n,2)
+ edge_2: The second edge. Shape (m,2)
+
+ Returns:
+ The combined edge.
+ """
+ adj = create_close_point_adjacency_matrix(edge_1, edge_2, distance=2)
+ # remove one directional edges
+ adj = np.logical_and(adj, adj.T)
+
+ # find the index of the closest point in edge_2 for each point in edge_1
+ idx_edge_2_to_edge_1_pair = np.argmax(adj[: len(edge_1), len(edge_1) :], axis=1)
+ # if argmax returns 0, it means that there is no edge between the two points
+ filtered_idxs_mask = idx_edge_2_to_edge_1_pair != 0
+
+ edge_2_keep_idxs = idx_edge_2_to_edge_1_pair[filtered_idxs_mask]
+ edge_1_keep_idxs = np.arange(len(edge_1))[filtered_idxs_mask]
+
+ # combine the edges by element-wise mean
+ edge_1_keep = edge_1[edge_1_keep_idxs]
+ edge_2_keep = edge_2[edge_2_keep_idxs]
+ combined_edge = (edge_1_keep + edge_2_keep) / 2
+
+ return combined_edge
+
+
+@st.cache(show_spinner=False)
+def load_image_and_preprocess(
+ image_path: Path | str | Image, target_size: tuple[int, int] = (1000, 1000)
+) -> Image:
+ """Load an image and scale it to the correct size for FSD layout."""
+ image = load_raw_image(image_path)
+
+ if image.ndim == 3:
+ image = image[:, :, :3]
+
+ image_resolution = np.prod(image.shape)
+ target_resolution = np.prod(target_size)
+
+ rescale_ratio = target_resolution / image_resolution
+ # we need to take the root of the ratio to achieve the correct scaling
+ image_resized = rescale(image, rescale_ratio**0.5, channel_axis=-1)
+
+ # convert to grayscale
+ image_gray = (
+ rgb2gray(image_resized) if image_resized.ndim == 3 else image_resized.copy()
+ )
+ image_one_channel = image_gray[:, :, 0] if image_gray.ndim == 3 else image_gray
+
+ # apply unsharp mask
+ image_unsharp = unsharp_mask(image_one_channel, radius=5, amount=3)
+
+ return image_unsharp
+
+
+def reorder_track_edge(
+ g: nx.Graph, cc_idxs: Iterable[int], all_nodes_positions: FloatArrayNx2
+) -> FloatArrayNx2:
+ cc_idxs_list = list(cc_idxs)
+ start_index = cc_idxs_list[0]
+
+ idxs_ordered = list(nx.depth_first_search.dfs_preorder_nodes(g, start_index))
+
+ cc_positions_ordered = all_nodes_positions[idxs_ordered]
+
+ distances = np.linalg.norm(np.diff(cc_positions_ordered, axis=0), axis=1)
+ mask_remove = distances > 5
+
+ first_index_over_k = np.argmax(mask_remove)
+ cc_positions_ordered = cc_positions_ordered[:first_index_over_k]
+
+ return cc_positions_ordered
+
+
+def extract_track_edges(
+ image: Image,
+ show_steps: bool = False,
+) -> tuple[FloatArrayNx2, FloatArrayNx2]:
+ # apply canny edge detection with increasing sigma until the number of pixels
+ # designated as edges is low enough
+ for s in np.linspace(2, 5, 20):
+ image_canny = canny(image, sigma=s).astype(float)
+ m = image_canny.mean()
+ if m < 0.03: # prior, found by trial and error
+ break
+
+ if show_steps:
+ st.image(image_canny, caption="Canny edge detection")
+
+ # treat each pixel as a node in a graph and connect nodes that are close to each
+ # other
+ edge_pixels = np.argwhere(image_canny > 0.01)
+ dist = cdist(edge_pixels, edge_pixels)
+ adj = dist < 2
+ adj[np.eye(len(adj), len(adj), dtype=bool)] = 0
+ g = nx.from_numpy_matrix(adj)
+
+ # find the connected components of the graph
+ # we expect four connected components, two for the outer edge and two for the inner
+ # edge
+ cc = list(nx.connected_components(g))
+
+ best_ccs = sorted(cc, key=len, reverse=True)[:4]
+
+ best_clusters = [edge_pixels[list(idxs_in_cc)] for idxs_in_cc in best_ccs]
+ # keep the longer version of each track edge
+ if len(best_clusters) == 2:
+ outer, inner = best_clusters
+ cc_outer, cc_inner = best_ccs
+ else:
+ outer, _, inner, _ = best_clusters
+ cc_outer, _, cc_inner, _ = best_ccs
+
+ if show_steps:
+ plt.figure()
+ plt.plot(*outer.T, ".", markersize=1)
+ plt.plot(*inner.T, ".", markersize=1)
+ plt.axis("equal")
+ st.pyplot(plt.gcf())
+
+ outer_ordered = reorder_track_edge(g, cc_outer, edge_pixels)
+ inner_ordered = reorder_track_edge(g, cc_inner, edge_pixels)
+
+ if show_steps:
+ plt.figure()
+ plt.plot(*outer_ordered[::3].T, "-")
+ plt.plot(*inner_ordered[::3].T, "-")
+ plt.axis("equal")
+ st.pyplot(plt.gcf())
+
+ return outer_ordered, inner_ordered
+
+
+def fix_edges_orientation_and_scale_to_unit(
+ edge_a: FloatArrayNx2, edge_b: FloatArrayNx2
+) -> tuple[FloatArrayNx2, FloatArrayNx2]:
+ points = [edge_a, edge_b]
+
+ all_points = np.concatenate(points)
+
+ # peak to peak (max - min)
+ ptp = np.ptp(all_points)
+ min_value = np.min(all_points)
+ # min-max scaling
+ points_scaled = [(points - min_value) / ptp for points in points]
+
+ points_center = np.row_stack(points_scaled).mean(axis=0)
+
+ points_centered = [points - points_center for points in points_scaled]
+
+ # orientation fix
+ points_rotated = [rotate(p, theta=-3.1415 / 2) for p in points_centered]
+
+ return points_rotated[0], points_rotated[1]
+ return points_rotated[0], points_rotated[1]
+ return points_rotated[0], points_rotated[1]
diff --git a/drawing_to_fsd_layout/math_utils.py b/drawing_to_fsd_layout/math_utils.py
new file mode 100644
index 0000000..c1d0560
--- /dev/null
+++ b/drawing_to_fsd_layout/math_utils.py
@@ -0,0 +1,638 @@
+#!/usin_roll/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+Description: A module with common mathematical functions
+
+Taken from ft-as-utils
+
+Project: FaSTTUBe Chabo Common
+"""
+from typing import Tuple, TypeVar, cast
+
+import numpy as np
+from numba import jit
+
+T = TypeVar("T")
+
+
+def my_njit(func: T) -> T:
+ """
+ numba.njit is an untyped decorator. This wrapper helps type checkers keep the
+ type information after applying the decorator. Furthermore, it sets some performance
+ flags
+
+ Args:
+ func (T): The function to jit
+
+ Returns:
+ T: The jitted function
+ """
+ jit_func: T = jit(nopython=True, cache=True, nogil=True, fastmath=True)(func)
+
+ return jit_func
+
+
+@my_njit
+def vec_dot(vecs1: np.ndarray, vecs2: np.ndarray) -> np.ndarray:
+ """
+ Mutliplies vectors in an array elementwise
+
+ Args:
+ vecs1 (np.array): The first "list" of vectors
+ vecs2 (np.array): The second "list" of vectors
+
+ Returns:
+ np.array: The results
+ """
+ return np.sum(vecs1 * vecs2, axis=-1)
+
+
+@my_njit
+def norm_of_last_axis(arr: np.ndarray) -> np.ndarray:
+
+ original_shape = arr.shape
+ arr_row_col = arr.flatten().reshape(-1, arr.shape[-1])
+ result = np.empty(arr_row_col.shape[0])
+ for i in range(arr_row_col.shape[0]):
+ vec = arr_row_col[i]
+ result[i] = np.sqrt(vec_dot(vec, vec))
+
+ result = result.reshape(original_shape[:-1])
+
+ return result
+
+
+@my_njit
+def vec_angle_between(
+ vecs1: np.ndarray, vecs2: np.ndarray, clip_cos_theta: bool = True
+) -> np.ndarray:
+ """
+ Calculates the angle between the vectors of the last dimension
+
+ Args:
+ vecs1 (np.ndarray): An array of shape (...,2)
+ vecs2 (np.ndarray): An array of shape (...,2)
+ clip_cos_theta (bool): Clip the values of the dot products so that they are
+ between -1 and 1. Defaults to True.
+
+ Returns:
+ np.ndarray: A vector, such that each element i contains the angle between
+ vectors vecs1[i] and vecs2[i]
+ """
+ cos_theta = vec_dot(vecs1, vecs2)
+
+ cos_theta /= norm_of_last_axis(vecs1) * norm_of_last_axis(vecs2)
+
+ cos_theta = np.asarray(cos_theta)
+
+ cos_theta_flat = cos_theta.ravel()
+
+ if clip_cos_theta:
+ cos_theta_flat[cos_theta_flat < -1] = -1
+ cos_theta_flat[cos_theta_flat > 1] = 1
+
+ return np.arccos(cos_theta)
+
+
+def rotate(points: np.ndarray, theta: float) -> np.ndarray:
+ """
+ Rotates the points in `points` by angle `theta` around the origin
+
+ Args:
+ points (np.array): The points to rotate. Shape (n,2)
+ theta (float): The angle by which to rotate in radians
+
+ Returns:
+ np.array: The points rotated
+ """
+ cos_theta, sin_theta = np.cos(theta), np.sin(theta)
+ rotation_matrix = np.array(((cos_theta, -sin_theta), (sin_theta, cos_theta))).T
+ return np.dot(points, rotation_matrix)
+
+
+@my_njit
+def my_cdist_sq_euclidean(arr_a: np.ndarray, arr_b: np.ndarray) -> np.ndarray:
+ """
+ Calculates the pairwise square euclidean distances from each point in `X` to each
+ point in `Y`
+
+ Credit:
+ Uses https://stackoverflow.com/a/56084419 which in turn uses
+ https://github.com/droyed/eucl_dist
+
+ Args:
+ arr_a (np.array): A 2d array of shape (m,k)
+ arr_b (np.array): A 2d array of shape (n,k)
+
+ Returns:
+ np.array: A matrix of shape (m,n) containing the square euclidean distance
+ between all the points in `X` and `Y`
+ """
+ n_x, dim = arr_a.shape
+ x_ext = np.empty((n_x, 3 * dim))
+ x_ext[:, :dim] = 1
+ x_ext[:, dim : 2 * dim] = arr_a
+ x_ext[:, 2 * dim :] = np.square(arr_a)
+
+ n_y = arr_b.shape[0]
+ y_ext = np.empty((3 * dim, n_y))
+ y_ext[:dim] = np.square(arr_b).T
+ y_ext[dim : 2 * dim] = -2 * arr_b.T
+ y_ext[2 * dim :] = 1
+
+ return np.dot(x_ext, y_ext)
+
+
+@my_njit
+def calc_pairwise_distances(
+ points: np.ndarray, dist_to_self: float = 0.0
+) -> np.ndarray:
+ """
+ Given a set of points, creates a distance matrix from each point to every point
+
+ Args:
+ points (np.ndarray): The points for which the distance matrix should be
+ calculated dist_to_self (np.ndarray, optional): The distance to set the
+ diagonal. Defaults to 0.0.
+
+ Returns:
+ np.ndarray: The 2d distance matrix
+ """
+ pairwise_distances = my_cdist_sq_euclidean(points, points)
+
+ if dist_to_self != 0:
+ for i in range(len(points)):
+ pairwise_distances[i, i] = dist_to_self
+ return pairwise_distances
+
+
+@my_njit
+def my_in1d(test_values: np.ndarray, source_container: np.ndarray) -> np.ndarray:
+ """
+ Calculate a boolean mask for a 1d array indicating if an element in `test_values` is
+ present in `source container` which is also 1d
+
+ Args:
+ test_values (np.ndarray): The values to test if they are inside the container
+ source_container (np.ndarray): The container
+
+ Returns:
+ np.ndarray: A boolean array with the same length as `test_values`. If
+ `return_value[i]` is `True` then `test_value[i]` is in `source_container`
+ """
+ source_sorted = np.sort(source_container)
+ is_in = np.zeros(test_values.shape[0], dtype=np.bool_)
+ for i, test_val in enumerate(test_values):
+ for source_val in source_sorted:
+
+ if test_val == source_val:
+ is_in[i] = True
+ break
+
+ if source_val > test_val:
+ break
+
+ return is_in
+
+
+def trace_calculate_consecutive_radii(trace: np.ndarray) -> np.ndarray:
+ """
+ Expects a (n,2) array and returns the radius of the circle that passes
+ between all consecutive point triples. The radius between index 0,1,2, then 1,2,3
+ and so on
+
+ Args:
+ trace (np.ndarray): The points for which the radii will be calculated
+
+ Returns:
+ np.ndarray: The radii for each consecutive point triple
+ """
+
+ # TODO: Vectorize this function. Limit is the indexer
+ indexer = np.arange(3)[None, :] + 1 * np.arange(trace.shape[-2] - 2)[:, None]
+
+ points = trace[indexer]
+ radii = calculate_radius_from_points(points)
+ return radii
+
+
+def trace_distance_to_next(trace: np.ndarray) -> np.ndarray:
+ """
+ Calculates the distance of one point in the trace to the next. Obviously the last
+ point doesn't have any distance associated
+
+ Args:
+ trace (np.array): The points of the trace
+
+ Returns:
+ np.array: A vector containing the distances from one point to the next
+ """
+ return np.linalg.norm(np.diff(trace, axis=-2), axis=-1)
+
+
+def trace_angles_between(trace: np.ndarray) -> np.ndarray:
+ """
+ Calculates the angles in a trace from each point to its next
+
+ Args:
+ trace (np.array): The trace containing a series of 2d vectors
+
+ Returns:
+ np.array: The angle from each vector to its next, with `len(return_value) ==
+ len(trace) - 1`
+ """
+ all_to_next = np.diff(trace, axis=-2)
+ from_middle_to_next = all_to_next[..., 1:, :]
+ from_middle_to_prev = -all_to_next[..., :-1, :]
+ angles = vec_angle_between(from_middle_to_next, from_middle_to_prev)
+ return angles
+
+
+@my_njit
+def unit_2d_vector_from_angle(rad: np.ndarray) -> np.ndarray:
+ """
+ Creates unit vectors for each value in the rad array
+
+ Args:
+ rad (np.array): The angles (in radians) for which the vectors should be created
+
+ Returns:
+ np.array: The created unit vectors
+ """
+ rad = np.asarray(rad)
+ new_shape = rad.shape + (2,)
+ res = np.empty(new_shape, dtype=rad.dtype)
+ res[..., 0] = np.cos(rad)
+ res[..., 1] = np.sin(rad)
+ return res
+
+
+# Calculates the angle of each vector in `vecs`
+# TODO: Look into fixing return type when a single vector is provided (return float)
+def angle_from_2d_vector(vecs: np.ndarray) -> np.ndarray:
+ """
+ Calculates the angle of each vector in `vecs`. If `vecs` is just a single 2d vector
+ then one angle is calculated and a scalar is returned
+
+ >>> import numpy as np
+ >>> x = np.array([[1, 0], [1, 1], [0, 1]])
+ >>> angle_from_2d_vector(x)
+ >>> array([0. , 0.78539816, 1.57079633])
+
+ Args:
+ vecs (np.array): The vectors for which the angle is calculated
+
+ Raises:
+ ValueError: If `vecs` has the wrong shape a ValueError is raised
+
+ Returns:
+ np.array: The angle of each vector in `vecs`
+ """
+ if vecs.shape == (2,):
+ return np.arctan2(vecs[1], vecs[0])
+ if vecs.ndim == 2 and vecs.shape[-1] == 2:
+ return np.arctan2(vecs[:, 1], vecs[:, 0])
+ raise ValueError("vecs can either be a 2d vector or an array of 2d vectors")
+
+
+def normalize(vecs: np.ndarray, axis: int = -1) -> np.ndarray:
+ """
+ Returns a normalized version of vecs
+
+ Args:
+ vecs (np.ndarray): The vectors to normalize
+ axis (int, optional): The axis to use for lengths. Defaults to -1.
+
+ Returns:
+ np.ndarray: The normalized vectors
+ """
+ return vecs / np.linalg.norm(vecs, axis=axis, keepdims=True)
+
+
+@my_njit
+def lerp(
+ values_to_lerp: np.ndarray,
+ start1: np.ndarray,
+ stop1: np.ndarray,
+ start2: np.ndarray,
+ stop2: np.ndarray,
+) -> np.ndarray:
+ """
+ Linearly interpolates (lerps) from one sin_pitchace `[start1, stop1]` to another
+ `[start2, stop2]`. `start1 >= stop1` and `start2 >= stop2` are allowed. If ns is a
+ 2d array, then start1, stop1, start2, stop2 must be 1d vectors. This allows for
+ lerping in any n-dim sin_pitchace
+
+ >>> import numpy as np
+ >>> x = np.array([1, 2, 3])
+ >>> lerp(x, 0, 10, 30, 100)
+ >>> array([37., 44., 51.])
+
+ Args:
+ values_to_lerp (np.array): The points to interpolate
+ start1 (np.array): The beginning of the original sin_pitchace
+ stop1 (np.array): The end of the original sin_pitchace
+ start2 (np.array): The beginning of the target sin_pitchace
+ stop2 (np.array): The end of the target sin_pitchace
+
+ Returns:
+ np.array: The interpolated points
+ """
+ return (values_to_lerp - start1) / (stop1 - start1) * (stop2 - start2) + start2
+
+
+def calculate_radius_from_points(points: np.ndarray) -> np.ndarray:
+ """
+ Given a three points this function calculates the radius of the circle that passes
+ through these points
+
+ Based on: https://math.stackexchange.com/questions/133638/
+ how-does-this-equation-to-find-the-radius-from-3-points-actually-work
+
+ Args:
+ points (np.ndarray): The points for which should be used to calculate the radius
+
+ Returns:
+ np.ndarray: The calculated radius
+ """
+ # implements the equation discussed here:
+ #
+ # assert points.shape[-2:] == (3, 2)
+ # get side lengths
+ points_circular = points[..., [0, 1, 2, 0], :]
+ len_sides = trace_distance_to_next(points_circular)
+
+ # calc prod of sides
+ prod_of_sides = np.prod(len_sides, axis=-1, keepdims=True)
+
+ # calc area of triangle
+ # https://www.mathopenref.com/heronsformula.html
+
+ # calc half of perimeter
+ perimeter = np.sum(len_sides, axis=-1, keepdims=True)
+ half_perimeter = perimeter / 2
+ half_perimeter_minus_sides = half_perimeter - len_sides
+ area_sqr = (
+ np.prod(half_perimeter_minus_sides, axis=-1, keepdims=True) * half_perimeter
+ )
+ area = np.sqrt(area_sqr)
+
+ radius = prod_of_sides / (area * 4)
+
+ radius = radius[..., 0]
+ return radius
+
+
+Numeric = TypeVar("Numeric", float, np.ndarray)
+
+
+def linearly_combine_values_over_time(
+ tee: float, delta_time: float, previous_value: Numeric, new_value: Numeric
+) -> Numeric:
+ """
+ Linear combination of two values over time
+ (see https://de.wikipedia.org/wiki/PT1-Glied)
+ Args:
+ tee (float): The parameter selecting how much we keep from the previous value
+ and how much we update from the new
+ delta_time (float): The time difference between the previous and new value
+ previous_value (Numeric): The previous value
+ new_value (Numeric): The next value
+
+ Returns:
+ Numeric: The combined value
+ """
+ tee_star = 1 / (tee / delta_time + 1)
+ combined_value: Numeric = tee_star * (new_value - previous_value) + previous_value
+ return combined_value
+
+
+def odd_square(values: Numeric) -> Numeric:
+ return cast(Numeric, np.sign(values) * np.square(values))
+
+
+def euler_angles_to_quaternion(euler_angles: np.ndarray) -> np.ndarray:
+ """
+ Converts Euler angles to a quaternion representation.
+
+ Args:
+ euler_angles (np.ndarray): Euler angles as an [...,3] array. Order is
+ [roll, pitch, yaw]
+
+ Returns:
+ np.ndarray: The quaternion representation in [..., 4] [x, y, z, w] order
+ """
+ roll_index, pitch_index, yaw_index = 0, 1, 2
+ sin_values = np.sin(euler_angles * 0.5)
+ cos_values = np.cos(euler_angles * 0.5)
+
+ cos_yaw = cos_values[..., yaw_index]
+ sin_yaw = sin_values[..., yaw_index]
+ cos_pitch = cos_values[..., pitch_index]
+ sin_pitch = sin_values[..., pitch_index]
+ cos_roll = cos_values[..., roll_index]
+ sin_roll = sin_values[..., roll_index]
+
+ quaternion_x = sin_roll * cos_pitch * cos_yaw - cos_roll * sin_pitch * sin_yaw
+ quaternion_y = cos_roll * sin_pitch * cos_yaw + sin_roll * cos_pitch * sin_yaw
+ quaternion_z = cos_roll * cos_pitch * sin_yaw - sin_roll * sin_pitch * cos_yaw
+ quaternion_w = cos_roll * cos_pitch * cos_yaw + sin_roll * sin_pitch * sin_yaw
+
+ return_value = np.stack(
+ [quaternion_x, quaternion_y, quaternion_z, quaternion_w], axis=-1
+ )
+ return return_value
+
+
+def quaternion_to_euler_angles(quaternion: np.ndarray) -> np.ndarray:
+ """
+ Converts a quaternion to Euler angles. Based on
+ https://stackoverflow.com/a/37560411.
+
+ Args:
+ quaternion (np.ndarray): The quaternion as an [..., 4] array. Order is
+ [x, y, z, w]
+
+ Returns:
+ np.ndarray: The Euler angles as an [..., 3] array. Order is [roll, pitch, yaw]
+ """
+ x_index, y_index, z_index, w_index = 0, 1, 2, 3
+ x_value = quaternion[..., x_index]
+ y_value = quaternion[..., y_index]
+ z_value = quaternion[..., z_index]
+ w_value = quaternion[..., w_index]
+
+ y_square = y_value * y_value
+ temporary_0 = -2.0 * (y_square + z_value * z_value) + 1.0
+ temporary_1 = +2.0 * (x_value * y_value + w_value * z_value)
+ temporary_2 = -2.0 * (x_value * z_value - w_value * y_value)
+ temporary_3 = +2.0 * (y_value * z_value + w_value * x_value)
+ temporary_4 = -2.0 * (x_value * x_value + y_square) + 1.0
+
+ temporary_2 = np.clip(temporary_2, -1.0, 1.0)
+
+ roll = np.arctan2(temporary_3, temporary_4)
+ pitch = np.arcsin(temporary_2)
+ yaw = np.arctan2(temporary_1, temporary_0)
+
+ return_value = np.stack([roll, pitch, yaw], axis=-1)
+ return return_value
+
+
+def points_inside_ellipse(
+ points: np.ndarray,
+ center: np.ndarray,
+ major_direction: np.ndarray,
+ major_radius: float,
+ minor_radius: float,
+) -> np.ndarray:
+ """
+ Checks if a set of points are inside an ellipse.
+
+ Args:
+ points: The points as an [..., 2] array.
+ center: The center of the ellipse as an [2] array.
+ major_direction: The major direction of the ellipse as an [2] array.
+ major_radius: The major radius of the ellipse.
+ minor_radius: The minor radius of the ellipse.
+
+ Returns:
+ An [...] array of booleans.
+ """
+
+ # Center the points around the center
+ # [..., 2]
+ centered_points = points - center
+ # Calculate angle of the major direction with the x-axis
+ # [1]
+ major_direction_angle = float(angle_from_2d_vector(major_direction))
+ # Rotate the points around the center of the ellipse
+ # [..., 2]
+ rotated_points = rotate(centered_points, -major_direction_angle)
+ # [2]
+ radii_square = np.array([major_radius, minor_radius]) ** 2
+ # [...] [..., 2] [2]
+ criterion_value = (rotated_points**2 / radii_square).sum(axis=-1)
+
+ mask_is_inside = criterion_value < 1
+ return mask_is_inside
+
+
+def center_of_circle_from_3_points(
+ point_1: np.ndarray,
+ point_2: np.ndarray,
+ point_3: np.ndarray,
+ atol: float = 1e-6,
+) -> np.ndarray:
+ """
+ Calculates the center of a circle from three points.
+
+ Adapted from http://paulbourke.net/geometry/circlesphere/Circle.cpp (CalcCircle)
+
+ Args:
+ point_1: The first point as an [2] array.
+ point_2: The second point as an [2] array.
+ point_3: The third point as an [2] array.
+
+ Returns:
+ The center of the circle as an [2] array.
+ """
+ y_delta_1 = point_2[1] - point_1[1]
+ x_delta_1 = point_2[0] - point_1[0]
+ y_delta_2 = point_3[1] - point_2[1]
+ x_delta_2 = point_3[0] - point_2[0]
+
+ if np.isclose(x_delta_1, 0.0, atol=atol) and np.isclose(x_delta_2, 0.0, atol=atol):
+ center_x = (point_2[0] + point_3[0]) / 2
+ center_y = (point_1[1] + point_2[1]) / 2
+ return np.array([center_x, center_y]) # early return
+
+ slope_1 = y_delta_1 / x_delta_1
+ slope_2 = y_delta_2 / x_delta_2
+ if np.isclose(slope_1, slope_2, atol=atol):
+ raise ValueError("Points are colinear")
+
+ center_x = (
+ slope_1 * slope_2 * (point_1[1] - point_3[1])
+ + slope_2 * (point_1[0] + point_2[0])
+ - slope_1 * (point_2[0] + point_3[0])
+ ) / (2 * (slope_2 - slope_1))
+
+ center_y = (
+ -(center_x - (point_1[0] + point_2[0]) / 2) / slope_1
+ + (point_1[1] + point_2[1]) / 2
+ )
+
+ center = np.array([center_x, center_y])
+ return center
+
+
+@my_njit
+def circle_fit(coords: np.ndarray, max_iter: int = 99) -> np.ndarray:
+ """
+ Fit a circle to a set of points. This function is adapted from the hyper_fit function
+ in the circle-fit package (https://pypi.org/project/circle-fit/). The function is
+ a njit version of the original function with some input validation removed. Furthermore,
+ the residuals are not calculated or returned.
+
+ Args:
+ coords: The coordinates of the points as an [N, 2] array.
+ max_iter: The maximum number of iterations.
+
+ Returns:
+ An array with 3 elements:
+ - center x
+ - center y
+ - radius
+ """
+
+ X = coords[:, 0]
+ Y = coords[:, 1]
+
+ n = X.shape[0]
+
+ Xi = X - X.mean()
+ Yi = Y - Y.mean()
+ Zi = Xi * Xi + Yi * Yi
+
+ # compute moments
+ Mxy = (Xi * Yi).sum() / n
+ Mxx = (Xi * Xi).sum() / n
+ Myy = (Yi * Yi).sum() / n
+ Mxz = (Xi * Zi).sum() / n
+ Myz = (Yi * Zi).sum() / n
+ Mzz = (Zi * Zi).sum() / n
+
+ # computing the coefficients of characteristic polynomial
+ Mz = Mxx + Myy
+ Cov_xy = Mxx * Myy - Mxy * Mxy
+ Var_z = Mzz - Mz * Mz
+
+ A2 = 4 * Cov_xy - 3 * Mz * Mz - Mzz
+ A1 = Var_z * Mz + 4.0 * Cov_xy * Mz - Mxz * Mxz - Myz * Myz
+ A0 = Mxz * (Mxz * Myy - Myz * Mxy) + Myz * (Myz * Mxx - Mxz * Mxy) - Var_z * Cov_xy
+ A22 = A2 + A2
+
+ # finding the root of the characteristic polynomial
+ y = A0
+ x = 0.0
+ for _ in range(max_iter):
+ Dy = A1 + x * (A22 + 16.0 * x * x)
+ x_new = x - y / Dy
+ if x_new == x or not np.isfinite(x_new):
+ break
+ y_new = A0 + x_new * (A1 + x_new * (A2 + 4.0 * x_new * x_new))
+ if abs(y_new) >= abs(y):
+ break
+ x, y = x_new, y_new
+
+ det = x * x - x * Mz + Cov_xy
+ X_center = (Mxz * (Myy - x) - Myz * Mxy) / det / 2.0
+ Y_center = (Myz * (Mxx - x) - Mxz * Mxy) / det / 2.0
+
+ x = X_center + X.mean()
+ y = Y_center + Y.mean()
+ r = np.sqrt(abs(X_center**2 + Y_center**2 + Mz))
+
+ return np.array([x, y, r])
+
diff --git a/drawing_to_fsd_layout/spline_fit.py b/drawing_to_fsd_layout/spline_fit.py
new file mode 100644
index 0000000..f6abfa3
--- /dev/null
+++ b/drawing_to_fsd_layout/spline_fit.py
@@ -0,0 +1,135 @@
+"""
+Taken from https://github.com/papalotis/ft-fsd-path-planning/blob/main/fsd_path_planning/utils/math_utils.py
+"""
+from dataclasses import dataclass
+from typing import Any, Optional, Tuple
+
+import numpy as np
+from scipy.interpolate import splev, splprep
+
+
+@dataclass
+class SplineEvaluator:
+ """
+ A class for evaluating a spline.
+ """
+
+ max_u: float
+ tck: Tuple[Any, Any, int]
+ predict_every: float
+
+ def calculate_u_eval(self, max_u: Optional[float] = None) -> np.ndarray:
+ """
+ Calculate the u_eval values for the spline.
+ Args:
+ max_u (Optional[float], optional): The maximum u value. Defaults to None. If
+ None, the maximum u value used during fitting is taken.
+ Returns:
+ np.ndarray: The values for which the spline should be evaluated.
+ """
+
+ if max_u is None:
+ max_u = self.max_u
+ return np.arange(0, max_u, self.predict_every)
+
+ def predict(self, der: int, max_u: Optional[float] = None) -> np.ndarray:
+ """
+ Predict the spline. If der is 0, the function returns the spline. If der is 1,
+ the function returns the first derivative of the spline and so on.
+ Args:
+ der (int): The derivative to predict.
+ max_u (Optional[float], optional): The maximum u value. Defaults to None. If
+ None, the maximum u value used during fitting is taken.
+ Returns:
+ np.ndarray: The predicted spline.
+ """
+
+ u_eval = self.calculate_u_eval(max_u)
+ values = np.array(splev(u_eval, tck=self.tck, der=der)).T
+
+ return values
+
+
+class NullSplineEvaluator(SplineEvaluator):
+ """
+ A dummy spline evaluator used for when an empty list is attempted to be fitted
+ """
+
+ def predict(self, der: int, max_u: Optional[float] = None) -> np.ndarray:
+ points = np.zeros((0, 2))
+ return points
+
+
+class SplineFitterFactory:
+ """
+ Wrapper class for `splev`, `splprep` functions
+ """
+
+ def __init__(self, smoothing: float, predict_every: float, max_deg: int):
+ """
+ Constructor for SplineFitter class
+ Args:
+ smoothing (float): The smoothing factor. 0 means no smoothing
+ predict_every (float): The approximate distance along the fitted trace to
+ calculate a point for
+ max_deg (int): The maximum degree of the fitted splines
+ """
+ self.smoothing = smoothing
+ self.predict_every = predict_every
+ self.max_deg = max_deg
+
+ def fit(self, trace: np.ndarray, periodic: bool = False) -> SplineEvaluator:
+ """
+ Fit a trace and returns a closure that can evaluate the fitted spline at
+ different positions. The maximal spline degree is 2.
+ Args:
+ trace (np.ndarray): The trace to fit
+ Returns:
+ Callable[[int, float]: A closure that when called evalues
+ the fitted spline on the provided positions.
+ """
+ if len(trace) < 2:
+ return NullSplineEvaluator(
+ # dummy values
+ 0,
+ (0, 0, 0),
+ 0,
+ )
+ k = np.clip(len(trace) - 1, 1, self.max_deg)
+ distance_to_next = np.linalg.norm(np.diff(trace, axis=0), axis=1)
+ u_fit = np.concatenate(([0], np.cumsum(distance_to_next)))
+ try:
+ tck, _ = splprep( # pylint: disable=unbalanced-tuple-unpacking
+ trace.T, s=self.smoothing, k=k, u=u_fit, per=periodic
+ )
+ except ValueError:
+ with np.printoptions(threshold=100000):
+ print(self.smoothing, self.predict_every, self.max_deg, repr(trace))
+
+ raise
+
+ max_u = float(u_fit[-1])
+
+ return SplineEvaluator(max_u, tck, self.predict_every)
+
+ def fit_then_evaluate_trace_and_derivative(
+ self, trace: np.ndarray
+ ) -> Tuple[np.ndarray, np.ndarray]:
+ """
+ Fit a provided trace, then evaluates it, and its derivative in `n_predict`
+ evenly spaced positions
+ Args:
+ trace (np.ndarray): The trace to fit
+ Returns:
+ Tuple[np.ndarray, np.ndarray]: A tuple containing the evaluated trace and
+ the evaluated derivative
+ """
+ if len(trace) < 2:
+ return trace.copy(), trace.copy()
+
+ fitted_func = self.fit(trace)
+
+ evaluated_trace = fitted_func.predict(der=0)
+ evaluated_derivative = fitted_func.predict(der=1)
+
+ return (evaluated_trace,)
diff --git a/drawing_to_track.ipynb b/drawing_to_track.ipynb
new file mode 100755
index 0000000..12157ec
--- /dev/null
+++ b/drawing_to_track.ipynb
@@ -0,0 +1,924 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from itertools import count\n",
+ "from typing import Iterable\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "import numpy as np\n",
+ "import thinning\n",
+ "from drawing_to_fsd_layout.math_utils import (trace_distance_to_next,\n",
+ " unit_2d_vector_from_angle)\n",
+ "from drawing_to_fsd_layout.math_utils import normalize\n",
+ "from drawing_to_fsd_layout.math_utils import rotate as rotate_points\n",
+ "\n",
+ "from chabo_common.utils.spline_fit import SplineFitterFactory\n",
+ "from chabo_simulation.lfs.lyt_interface.io.write_lyt import write_traces_as_lyt\n",
+ "from scipy.sparse.csgraph import connected_components\n",
+ "from chabo_common.utils.math_utils import my_cdist_sq_euclidean\n",
+ "from skimage import io\n",
+ "from skimage.exposure import rescale_intensity\n",
+ "from skimage.feature import canny\n",
+ "from skimage.filters import gaussian\n",
+ "from skimage.transform import rescale\n",
+ "from skimage.transform import rotate as rotate_image\n",
+ "from sklearn.metrics.pairwise import pairwise_distances\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 119,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# !python -m pip install --user --upgrade pip"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Load image and basic transforms\n",
+ "\n",
+ "- Change the **path_to_image** variable to the path to the image you want to load.\n",
+ "- The image is loaded as grayscale.\n",
+ "- The image is resized to a (smaller) size to aid with computational speed.\n",
+ "- The image contrast is boosted to improve edge detection"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Rescaling image by factor 0.063\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAD8CAYAAAARze3ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAC+zElEQVR4nO39aZCe13klCJ6b+45MAAkgsREgCZAECO7iIlosihQtyZZMl7yEXFEz6gpNaH54uqq7Z6ItdUdMxfxQhHumw9E9MeOJUbjcI0/bpWLLJVuWbJEUyxQlS9zFBQQBEsQOJBJAYsl9f+dH5vPyfCfPfb+kbZnJnrwRX3zf9753ee5zn+c85977LqkoCqyltbSW1tJayqeGD1uAtbSW1tJaWu1pDSjX0lpaS2upTloDyrW0ltbSWqqT1oByLa2ltbSW6qQ1oFxLa2ktraU6aQ0o19JaWktrqU76hQFlSukzKaUjKaWjKaWv/qLaWUtraS2tpV90Sr+I6yhTSo0A3gHwGIAzAF4C8DtFURz6R29sLa2ltbSWfsHpF8Uo7wVwtCiKY0VRzAD4FoDHf0FtraW1tJbW0i80Nf2C6t0G4DT9PwPgPs6QUvoKgK8AQHt7+907d+5cVklKqeYbAIIBp5RqfvN3URQ1ZVyduVTv/N837z9FWqk8rKf4X6WzlbaxktmJtq3lq8a86jifr1f3P2UqimKZnUZydsvf/LtKb5rceOZ8ReXM6eiDzDz/IXp2ZXN9cfmq6qxX7uzZs7h8+bIV/hcFlK6xGimLovgGgG8AwM0331z80R/9URxHSgmNjY1oaGgoO9nQ0FCejzxsEA0NDWhqakJKCQsLC2UZVnBDQ0N5bH5+vqyDjSfyRDus3DgHAAsLC2V+TlH/wsKCBXXqPxobG8tz/NG+aTsNDQ01fQx5VIZlg0J9VVDituIYl2P9ON1wvg96fGFhAfPz86XcuTHncvE/xoH7xzYAAHNzczXni6KoqVv7yIClDloFejzmWiba1ADPfS6KAvPz8zX1cD8ixfmGhgY0NjbW5GGZXd1NTU019lEUBRYWFspPpMbGRisD94XLq/6j/YWFBTQ2NtbIpj4V/+fn50u/j3qjjvC9+fn5smz8Zp/mfoQc0ZfIH/95PBcWFvD44/lJ7y8KKM8A2EH/twM4l8tcFAXm5ubQ1NRUgp1jETEI0WE2AlYsgwobjjoVn1PH4/M8kOo80SbnVcfWOlm2+Fagc+zKgWKcZ/Dk+kM+bZPr5d+qEzZ+BgOnO9aB6llld8FDg486qJ53AUVlVHBjWbUfPH6RAsBVBm2f6wq5OS+DBX+z/UTAYHt1KcY68jU1NS0DOh5XPhflHFAw6WBA4jFweuT61F+AxUCVY5kMhqzzOMf2wn2em5tbBuScX/1e22Q7YXDPpV8UUL4EYE9KaTeAswC+COBf1CvkjDeO82CE0QHvKzrYmTITZQDK2Fj57CAR2ZQlONnYuBUAuH7+zf81ErO8nDgQsHGonBzFOQ/30bFOBVFmj+68AhTXwwDAZdWp4lx8mpqaljmcBgwGXnVqFyy4Da3HBT7nYBysNLApcESKMVCGm9N7Sqns/8LCQgmAbF+RXPBjO1C2ymyWZ1PAIjDxMZ6tOPDN6dH5p/quY/IsM48lj60e18R1uEDF8jBQrjT9QoCyKIq5lNL/AcCTABoB/HFRFG+tpGwIr8wylKyA5wYFqJ1mqrJyDE/zuPbUQHnwWX7nWG7QuB8KRGqIUX9jY6OdGnF+ZQRRTtmqtqfn3PTXBYVcUKrqjzuvzqTAmgNNJx8HDzdGWs4d47bDceM3g5j2mQMN18vJ/Vc75/5w4ulm1VIL59d2lBVqUpB24xr1OMKhOuf+qV+HHzP7dAHQkZdgkrp0o8yUg6sGww+LUaIoir8G8Ncrzc9OylML7RA7ujq3GkCwzEi6tsfMix1dp6txfqlfZXk9ryxLnVDPcVsOFLhulp0jPxtTY2NjjWHzFJblZZ3yUoXK5Bwk55DOefi8TqPdb9WNsl/WjdORa1vLOnmd3M62oj7Oo+tduWUObZvrYTZUBaxhzwxKyvr4/+zs7LJ10QC1XFLdK1PkNXhO6g/OjrV9HjNdkmH9MzBqwOT/jgApoDLGOFDNpV8YUH7QFMbjBj7AjCODKiwAIo7r9NOlpqamsm3nRDlnUmCJcvWSi/A5NsTtcf9C5hh0Xq9zxsR6Db0oU1ZdattOD5qPjTXHnNhI1RGjD8rAVS4XEF2gdIGUZWG9OVl1qpxjddzfHHBz4jxO726JROVw9hdJ1/zcco1bL9bk2CaDFtuSBo/4jjX9aFODXW79nOVza9pqw7rW6AgK61fBN/JU2fmqAUpgebTTdS51Rt41jvLA8vU3xzLUEZUV6LSNoyjXx0DnIqgyCDYkHjReZHfLCQowrB/uO9etOmhoaKhZA+RzUY6/uT49XsUitbyCOrev62rKJFgnKwEfZRtRjzJnlofHlc/nGFHIzb/VobUN7psL+K4NLlcvaOnSjZMp+pqbqmv7Sh5YJy4Pl1e5dTbExzUgOqBj+2d/VLthvWmfnM+zbVbpZNUApQrp1lTYYTmacD7HFFSJDDA8gAqQCnAMsCqXMhB28sjrgDAAwBmamxKrboI5a5BRPTCzYrn0vOpOWY+yOVc+5FKAVwBQx+ayPG7KEjQ48fgBy6dRWk/0g0GDna6q3yq7OrECGY+ZCyAahNhu1Oa5HNfLG5ncFq+jxu8cG436mQW6vrKMrDsNRg6kVa8K3iqLY4U8zpGif+yv3BbX6UDf9VvTqgFKBrlIajDOkasAEfAOyG1qPbk10Kq2+VwOHNnJIy/Lpn0Ped1ShP52ded0lmOT3Kccc9Ox0Hpyjh061WtXNcBEn1kWriN3zqVwZr0sRx2XZee63XhwfpWfA9VKdOj6UY9duk0It97IG3fxn8GIbT0Awi0FsaxV9uTyq++ofytIar6cH7lvXZqL344kKIHJ6VDTqgFKnYYwU3JOrg6Xo84KXpHUcLhdB9ihTJ3uu3UWlpEHXdmrkyunm6rBdCDCdfNOrfaryhHqgatzagYoZSK6ZqysTYEzZ+AqayQHpI7BuPUyZZaqE2WMjt3mQC+ChK6bO/B2AFNFFEIG3aRU3bJ+1Fdc4HPjofk5hZxuhhbn3RopB1DVJ3+q/Mzpjv3ZsUoe66qgGGnVAKUTWI1cgVGTA8scEOluN08ndPoQ+XIGHf+5vZBFB6feoLipS3z4cgddg3LTUa5PHY7zOIZTxSarWLoz1ioWq2vSmjSIsR3oUgnX75yDE5etxzJYv8FatE2uT6d//K3LAsq2FKhZdh4Dvshc9cZLMRy0dOqtrJRtLY7zcpOzk+gTy6obJtx2zk71v7bDdsXLOspiXR2cX/M5f3BpVQClMsIYnBjA5ubm7JTYschIDFYKxOpMGm2rohXXEWWA5bcR6uaLAxgFztw6I8vv1puUQeTYB6eccdRjUithNgqszulVB8qccutpnNexUK3LOYhu4Oj4B9PR20y5/dAPM16Vj20j6mNZNQhwee5P3IbJfqKsjceB9evWbHNjz6DI635ch/Y3bqNkG1G74c04DXghq84gWY8KkLpkwH3ivKpP1XU9gIy0KoASqDU8VqQDSD7GxupAwUU0x1xiIHVKzE6igK7fbqDrAVe99RGOzCyPTmOrgLwey2ajYoDPAa8LUOrobj2Kf2sAYifKRX/+VkB0THilwc6BePRBZYlvF/SYKQWY6LpgbDxocKj3m+3WLQ9xmp2drQEoZc/RN07sTw5EGCyZgbI+VA9ar7vDSPuq4KaAGu2z36usXF7bUHKhv3NpVQClDlJKqbwzB6hdfGYndtdXsqLdtNCBhzqby89OkxvslTBGbZOTgjivi2qdyqZcH+oZX/yPunXjg+thgNT/kYcB1xm49t8BWI4ROqahfVMWowDMx9TZdQ00zinwunELhqPjpf10AMltxQ0DLJeWd33htphFBtNzm3xcBzPC0OHc3FzNBe4awAL4dG9B5eU2GWxdoOB+6BSdZeF8Ts9Of+oX3C8el1xaFUDJzuYcPTqiDBJYHik08XF1Zk4u+uQcXx1K+8IDHMarRqaDF33SyOwu/cktUHP7VQDF+TkxELvNJmXyLiJrwMgFJl4+cHJWHXf1KutmXagsOfagoOI2luJ3VZBjG3UszrGtaJ/PK3PUjRH1A65Py3Gfc/bD+dzVAtxXbtO1wXpTYI5lCJ2dcR94z4D7G3kCdJW8qLz6Xx9sowQgl1YFUALLFcUKyTl9DlgjORaowBnfLmq6SMfMwzELdzzKal4nJ4OIyu2iqUZtjaq6s6d65naVgXNf1PHUELXfuafluPFSmbgNljOYj7btwNONgZOZv7kfwRDjTihObu00xzxZDnb++K924YBV8/F45cBX21RwcLbe1NRkryfkoJUDRNa9EgSWVZmijhnLmltHBrBsGu8CNTNTXTZgOaqIT6RVA5RALYN0xq4LzcDyAeK6dH0o51gOnDk5R3IRq8rxHPiy8XI+lZN1ooZVxaICXLhfrs8KhrlAwbqN/FG3ey6ibnQoGHMfXCDRKaYLQjkAib4D1UswXL9uMjDYs4w6pi7AcN0xpXYbMdpvZoqsr8ivSwZ86VcuSKteIvE0mBlayKE2HvqJ8eY6nc/oOX0CUM5n1CacT7JsPB4KgDym+pANDc5VaVUBJVB72Y5+q2KdsTkgUoUoKDnn04HQ49yWA0Mnk5ONnUcdLv7Hw1ZdoFDg0ekzG4T2y4EWyxf18P9IVcFEHYx/K2Ng0NW1Zh5bN45ab/x3DITzuMTjoVNTB2YchCLpfx1jdloFiGiPwZnlz21Ccl3KEllebcctQfHY8NIIl4//c3NzyxhiboOK111z48X+yLM7ZZe8e+6A3S1HOJnq6ULTqgJKZjTsLAo0umkTCo0FbL0UgcsCy9dzmP3Ef43sHEXZKBR44rhGNufcIR+DXrSha3gxmDqoahTcL7cxELIoADg54zjLp+uKyl5Z5zw2rn43rmzcet5tYjhmyEmBxAG8A3FlNRpUGZTiE/oBll/LyMyKQUBlUH2ETURySymcX/vF59ke9IYADrA69Vcb5wfQ8LjxmqEClPqhTuGjbmeTmkfHic9V7UG4AMh914DIadUAJQOjLhazMUZSRhjGw4ClQBflIuUGhQ2fyzgn0sFyIMpJI1u0U7ULqLJG3W5gOSrrBlicz7EJF3z4f25XXGV0mznKPHg8VB+5fur6so6TthVJr/Nj2dzGBrfBcrlxU53xODLgqD0pOCp7jHz8hCs+V+XULtgVRVFepK79zQUdLa/6UsBnIpFbT2R/VUDltlhP8dsFq9CJ+jL3g+VnX4iyKaXspUtlvuyZf+Lk1kTYOOuBB5/nyM7n3LcbLC5X7+Pyah0srzoJs51I/G6PnDG7abNjGSvpv7IsbkeNlR/+q3Xx8wo56OVuXdTNCtaLzhq0L9y+Ami0yWPCATja5HZ1ysayqj3yp0rm3D3IDuR0w4TBP1gcz1ai/cbGRjQ1Ndkpu7NVt9anrIpldBtJrBf9r/Xqh/Wr64Q54qHjXQWmueUJrpf75NpwqS5QppT+OKV0IaV0kI6tTyk9nVJ6d+m7j859LaV0NKV0JKX06Xr1l4Jk1haio5HUUSOxwt1GT87xtH6d8rOB5XY7ebC0XQWs3DnOo2tdfM5Nr7g+NXzXR9e+PuRY+6JjkqtbA0OuDgaqOOdkVjkdC8vp0tWh7Ejz8rqbrovlArkG3JAtAFiXlFwdrDtd+9Nx551p7o8CQQRcPufA0/UtB+hOXjemapv8HX3QzT/9xFqotuXyxXdKadklQKrHhoYGzM3N1Q1enFbCKP8/AD4jx74K4JmiKPYAeGbpP1JK+7D4fpz9S2X+MKVU/fTcpcRGqWtgbHAOlHI7b3oPqjqTOxdtcXRTIFUjdrLy01ty7ebYVlVyywlsKGqc/Hslhs/G5NZtHfPPye4CnrLeKrbjmIMaPY+/k0ud18nqlgJYB7lNkqjHrelxoHXgErJpENdNuip2F8cUdKI/jvm7lGN0bn2f24jkAlEOiJ0ec74dutGgyuccFrA+1Se0DdbjP4hRFkXxHIDLcvhxAN9c+v1NAL9Ox79VFMV0URTHARwFcG+9NgCU0weNqtGxAD0XhdiJ1ODYQVwUdmtDOiD8u4pZcX4FV066Xsj9cdOCoihq3jrHsrik7FmjqzMY1oVjTJF4WqjOwWVZRgcWPL7aVg6QVH8qv4Jq6JLZQ85h3dgyaOrYcB63zsV9Z1mZYTrmGnXxLbWqc77/PLcBouOvds5tu6UqJ7cDOJZBx0+XYpSMsB2GfyuwOxDmfuk4hX70tRUaLN1rYqqCyd93M2dzURSDS8IOppQ2LR3fBuB5yndm6VhlyrEBzROJFeTW14Dl16NpWTcgVYyLy+WiouaL4/G9sLBg3zLo+uj6q4akGwdVkVuZoGM5CjQsP9fBAJBLnE93m11QqwJazR+OFjKpjjgps+K6debAuuF8Oq7Rn6o6daODGQ7n08ClgOXW3HRTSHWq46j9c0sOOVbJ47cS+9Rzrm/Oz3LjF+1rIHU+DWCZ/rl/qjen21z6x97Mcd5jkS+l9JWU0ssppZevXLlSExk5ErCjANXrf3yedzr5xei5qOXqcUanTFXBut4gaOTlKbrLqxs7Kr8abBWIufOhGzZANjxlIq6//DvHPDhvLvhoPfE7ZMpNo7lvTveR103DlFkoWPBv3pxyQYk3b3jsQgZeM3SsPOqJfvIGDveRdcr9U3lyAOzWdnXM1J94eq/nlJ2qvnnJhm1EA7KzHWfPGjRDptBJvH5X+6yBTINfVfr7AuVQSmlgScABABeWjp8BsIPybQdwzlVQFMU3iqK4pyiKe/r6+moi7lK92TW+pfIWSIHahxIsLCzY3eOoo4pJ5CJvjuE5B+Z+uHUdni4444+ka3uufu4f90tZmTM8Zyhazu1Suh1W1qdGf+ecuZ1VzqtOqs6nulLGkQNjty6qIKMPltWpvwJ0FbgpgALVbwDMTYnV9tX2mCCoLAoUrFcNmI6FabBxfpBbF2Yik5uF5caD63B9U9liuaooCnsFAv93vsvp7wuU3wXwpaXfXwLwl3T8iyml1pTSbgB7ALy40kqVLrOCIkrkOsTHI6/bYYzv3IaLAqkqUVmQA0+NVlo/M5VILBPn1acoOaBUp88FFpYt8jhG5diSysvTKBccXOBw8isIq0M4MMwBvuqD9a+60DK5YJFbz8v1Kz56h1nojJchGDQYMDjQ86YSy68BywFSLkAqsGowy61BV40r54tjLqCojKq/AEFlhdwGjwHPFmdnZzE3N1eW1fHl+rQP9Rhl3TXKlNK/B/AwgI0ppTMA/i2A3wfwRErpywBOAfitpcbfSik9AeAQgDkAv1sURfULc6kzjo1wJ5QZ6JSRjdKxOI7s3JZ+O9DknUSXXwdd++WMySXeZY2Ue9uk6koNXkFadajrj6obB16cqpiFHtf80U/OW08v2ja36Zy+HlgrGEefGehUBmVYGlQ4X3Nzc+mg3Adds+X1StVPbn1XfcGNe7Sl7Wg+V5/WE3Ip+83pxNmqWz/kurW93Fqw9ovvCNIpPrel69pMwuqlukBZFMXvZE49msn/dQBfr9uyJDW2XBQLBQRoOTanO1661hNJ1xt1kNlYHQA642JjCVkciDuQcIDrANABkG5MxLkAuxzb5PNud94BjQMlnQ3k1nI1wCiYOZ3mxkqTCwjx25135VydGjw1OKuNqGxNTU3l9YDRF94gya2psp2oc2s/c+v6IaPeb60gz+MYbfN/btsFGdUJB4gqu9QArWzdkaPQnc7ccoCv+uJ21C5zadXcwhhJBxCoZY5sMKpkdYycM9WLwk4WNkitz4GctqdA7Aw2B4pOJj3nAIjlc21zOVeXjoMClDoXAyQnt4aqzEL178ZQ+8THXWBVfSgAad91zYtlU6DUfKpzro/1pOCnl7CwfA6sFDhzd+SoDqP/yuYcYEe+CJ58XH1Qx0hBivPqlFhTbsx5qu3GihOPNTNSHie3blsvrQqgVANj4FN2xoahA8UDxE9o5sHWp0+zobh6nePnnJTL5xb8o4wDQE1hHPw4fy2jjqYMyLUX/Y6nEmkkVtlyQcitA7r+KtC4AKQg7pgW4DdKqsA/Prrj7NgMs3IH3C7IKnAo6KkjVk2FeQYD1D7EgstWMSaWLa784P7m1u2dnnP9ZUbGfdAxZR0z0KludaydDXPf1FbckkL4QtVaOverHliuCqAE/PV6kVKqvfRC13dUyTy47q13DChA7RqHYwZRLhwuB548YAHsuaTTeQVS1xfHlNgp3FS5oaGhxuG0b8oe+Lwa4EqWC9hAnU6VhbL+FIwUSF3S9vmbk4KgY1mqdx2THKirzhkgNHEgVTmVNakMqj8HMM5m3Tqi1sV6Y52orpWJsy+pnNx2+EK9h0/oeLuxVXvQ6br2wwE5t6XM36VVA5RAniUqmOi53GK4DkisF1aBXTC4qEeTY5u8SK3l+E1+LLNjIs7pVRcODFl37Pyqv5wDOofQNlQ+p1eVxV1Wk3PCSG7NK47zuLjpYI4BqZ60Hc3PwdgttbCOVHZmYlVLC1yfsx32BQVQ1y8t4+SNPMoCuWwOPAEsW+eskiHy6C3JvOzBtsFrq1p/LhDmAjUHfy7PdTI5WklaFUCpbCL+O6Nz17E5A3KRM0BL1370EpmcMQGwMrlI5YyJDV/r5Tw8yAqw3B8+r33jfuR2jKO8Y92qX9ZJzoCVfebYUPx2BqxO4wIe52M9qB25mQnrhZNj/1qHYzXaXm7ts55Dcn9YBzzFzS1FaOCItvmZmDoWHAh0huV8x/lHThdxXPWjMjArZhKkOmUCo+yYZVB70Odmsr5SSjW3ttZLqwIogeXOq+uGLqpoRGAlu3cZ67qLA6RI+uh+BW42xKgr5ww5NlFlWACWGWMVA1S2psChBhv1uwXvnIx6zDGFemtN2mfVn+ox57xah9OxOkmOCamdqe54fNxUTusJfeQCvm4mRDm+QiJeH+HYk9Mft68BmcdIxyxkdEGF267awFH9VrFNp3PVDQd/xxZz48X5uA96mRCPi07bc2lVAKUyAH1JvBobK5fzuaS7jFyHAqfWxWt3zGKd0ypQuYHTHUq9g4iBsSpIRFLn4U0wBUUtozLGfw4Q2kcHejxG/NvJtBLQ4nxu6sR6cM6mjCTnWOrwzAZzs5YqxpULSGHP7OT6Dh0dYw2UjmnFeZXDXerixpfBT4OV6obrZnteCcBw3fzNsnHbOiNy4+SW2nLj4AKT02U9ZrkqgFKdKEBEHSI6zm/GcwvV6lCq/PgORfGTmTnasKGyUVetb/DxlNKyd4s4pqzgzX1xIK6G62TgvA6oua3QOZdl3ao+Hfi4qZ/KxGCheuK6om0GW+6PA8wcgPP/qDf6y46vrErrYEd2IKmg4fTB4+2u3nB94P7lgibLoWQgx15ZppX8Z72o3bA9KWg53XD9rg+aR0GW83D9/M0vu9NgqOC4EsBfFUAJYJkRxjFg+RqWpqpdPWZqXCeAmk2dXL3s1Dlmk4vmOTDKMTYGU5dP9eXAxkXGqFOdktmvYyJVBqQgFm1rpNZ6FGi4H+4xZQo0ynqiHnbUKkBVOaIPbkxcINLx5W9newzIORDK1alrljndcJ08FnHMbQYpKXD1rQRAnP0xq2N51TbVHt01lrwkE0kfNhJJAzsTBF0H1fbr9XXVACXwvrB8y14YMrNMBZxcyg1QlNVrKuN4lI1vjcL1WIfmcQas9XF/2SmU2aoRsvw5kHBvtnRg5eTN6cAZGYO9C3gMTlWAxmPmZNXf6pSqbxfsWG8sV65ffKyK0TDA8vKKtq8M3t0a2NDQUN6f7JaQVL9OdvYRHVcOlNon1a0LfNofzrewsFBzxYcbw/jPSx5cf9V6sANdBWkOEtx/DbYr2dBZFUC5UidUZTNwVNWnSZliJAVLPa6/VeYcY+B2ec2EQVmdNdfnnINUye6ipm7g5NgJjwfnc+uHatQ8tdFdVpatihlxXQ5kta9VLI3/c326K8xgp3W4eoMNuSCaA9ucbLw0oMCjulH7Z8akYJfTl/6vshUFXB5HBjeWUe9zZ1/JBZzcfyVPmkf76kDZ1cs3I+TSqgBK7UAoVdfUNOUMQI0t8mpStsHGnJseu/LaRhxTA3ORKxw11xeuV4GYwYj7y4YUSYOLMjYHbo71uXVLl5f77wA595/ZmNuRVPakAOfYUy6guODGQM5TXu5DACMf075G/5kN5pifu1JA+8DjyKCs48HjrmCrQcEBjgNeHVu1ucircjqZ+Zyza9df1UtuTyKOuScHaV26yVMvrQqg5JRbL1KHDCXHOiMbONcT10opYLldTTZMBi+WR9eKFNwYJON/LuU2CByIKOAwQ+MXMKmcMX1jPXLS9lnPrEtXTmVzDsVjo5t0XEZ1ouMeyYG/W19l2bnuHAhpkATyd5FEuw4kFABZD5yPmZbbZHPAoH1yBEJ9RaftXIaBM/6z7XBQ4L5pgHb6VHkiOVCv6kP81w1Zp1P1Tz7GMoYtqi9XpVUBlDlG56KJsjRlAJFYgc6RlHE4gHDA4nY9cyCuifNXXRoUg1h1L7bqQEFH6+S+quGo48Sx3EXY7moA3ZnkAKNtxVULyqRccOSx1eUCVzfXp/81wHGfWW+qBwbj3BRYj7mkbIcTgwM/aagqsR50NhHnlQECtUCr+fVbfUqZrLJ+BlsFZfXB6KdjsyyzewIXs1S1CbZ9BkQXGLTdXFoVQAnUAhCvIbJhqhKA/GOh9LY6VgRf4a9O5pSozqJAlfvWxIPqBj4MU18exY7t2ufkNkH0PMtZFUC0LS3H4Ka7yNpnBfXoJ+sjxkbXNp1xK/tT/a4EPLmP3G+9N76KPbryqsfQu5v5cH+i3dyYs8xcL9+R5RiaSxqYuD4FHZZBZVN5FBzVj7ivPOtxgSnkUd1rwGLmzQGV8ywsLCwLQFq2Kq0aoORIFf+VVTDLUSfWaVJuSh7neBBzUdMZk2Nk2g9Orn2gFuB5SqpG5gCZ9cN53bRDgw+3rbI6tsb5HKjyOe0Tl3EOyeOdewc1j3Mc47U/LqPsRZ1a++H0xW3osVxfnAwAygvLuc+ceFlE64n8ufHgAKOMzOk6zuvsStfHGVTUL/i8bhryu52ibzzLUDm5XtYLt6FjyAxZ2SEvMbAecgxex96th3Oq+yqIlNKOlNLfppTeTim9lVL6N0vH16eUnk4pvbv03UdlvpZSOppSOpJS+nS9NkJBPBislKU6l7GHKBf/1UhDAcDy3WBlqlEXO5VzMGZDLLuCKw+Ctq1Gn3ttABuBi9rKHlRmZuacdA1M+xjRXo2bp1qqfwbiCGjK8nM75c7BQ8fcBx0nJ58GB6d7Z2fatmNmCjCspzjGnxzTUudW+eLbyefGjFkX95vJh26SMaPjMVJgVh90fQ+fAFDDipVFVumU+8k6ZH1rfh3fkCvnH7k2Va8ureSdOXMA/o9FUdwC4H4Av5tS2gfgqwCeKYpiD4Bnlv5j6dwXAewH8BkAf5hSyj9vTJIahzq8YwA8GFpX7i2HyrDY+aNOHaxI6ggsF5/T9TuV0d2qye1r/x0oquNzG1xW23fMQ/sYid9XxLI5WRXM1FFZt+ocytrU0eK46tz1TZP2ux4I8nEHeioTB5g4rk4bKd5hH8lNmxmAXNBVfwgdq6/k9ObGX8dQ31PFAYh9Q21cp8WhwwhA+i4c9hkGZtZ/1OtsXcebgVyTC/z/WK+CGAQQ7/AeTSm9jcV3dT+OxXfpAMA3ATwL4PeWjn+rKIppAMdTSkcB3AvgZ3XaWWbMoUQeIGYC6nwa2dhoNAJVOZMCkjqJtqt5nYFHYtDmnU/uoxo564O/WWb+zq2BqSNUgaXqTB1N+6nMgfurxp9jcq4+Ls+gwAbPgKzB0wGuBhenB7dOqPpTfelSgANlHn9l3FFHbr2MfYR16oK+lnPBSP1L+xIg4m7M4LLK9pSlRnLPReV6HcBG4jZ4LPR/yOWCttbFeetNvT/QGmVKaReAOwG8AGDzEoiiKIrBlNKmpWzbADxPxc4sHatMbg0tkk4rGERdVFeQcoag/6vYSI5BKPioc4YsVRtLju3qwKqxch5ncC6p4bHMGp3rsSfXDz2vfcjtjmtZLccpB1z83/WR9R9g7NZp1e60fu2/yqp2xGCTC9DO1hUsc0FI+50LQK5POYBj8GbdsV9xHtWBBiR3WyrrywUwXQbI6Z/XtZ3tKPDmGHC9V1QAHwAoU0pdAP4cwH9RFMVIzpgBuBPLRi2l9BUAXwGALVu2LGYSoFnKZxkER40qJ809pDfq0Lp13aOKEfB/HUjH9LhtIL8zr33i4ODqc4CRy8usQiO3Gq3KzHXz/bNV8vACu3PIKgN3gSf+T09PY2Jionw96ejoKGZnZ9HR0VHK1tjYiM7OzvKYezWvtuPsyW36uADpdK4sk5Nb49PkgMWNsRtDBgetMxfcdUkk5GSQ5LVALqsyunHm/vA4K2MEameFLIvmc/pyMxi2fe1fvbQioEwpNWMRJP+0KIr/uHR4KKU0UCyyyQEAF5aOnwGwg4pvB3DOdOgbAL4BAPv27St7rAvcDjwY4Pi3Rm1lL01NTcsGfUmWmnKOKeRk0RSGpGxKBzgcWTcOmA1XyaODzP1S51VDzumaGQIDdORxActFc8fgFJD5OOAD2sLCAmZmZjA7O4uzZ8/i4sWL6Ovrw8LCAi5fvoyrV6+iKApMTU1hbGwMCwsL6OjoKINdU1MTmpub0d3djZQSWltbsWXLFrS0tKChoQHr1q1De3v7Mh1FYuBQW1AbrVoPyy0J6Njyg3a57pAvLm9RRsxj4dgd/3eAzyCiwBtyVbFVHk+Wl1ml6lHl58RLUsywNSCx/1fJl2OLuY0nl1byXu8E4N8BeLsoij+gU98F8CUsvuP7SwD+ko7/WUrpDwBsBbAHwIv12uGkTs+dYYDj38D7TukSGwP1zYIQg52bxueiOTsRM+HIzwvgOhV1a3OcXGRWINXImdMDG5U6COuf23Dso17AUMNTAK9adyyKAiMjIxgaGsLo6CgmJiZw+vRpHD16FI2NjRgfH8fMzAxaWlpK4GtqaiqdcnZ2tgSXM2fOYG5uDk1NTejo6EBzczMaGhqwefNm7Ny5EwMDA2hvb7ebZvz8xtAN94flDb0FsOj46vjx8YaGhmWXtDndKpByoOGZgFtacDMYtZ0qpsb9ZHDk86oTltcFbLdBE7rQTRwH4kqK4jzrwcmm+szZcaSVMMoHAfxvALyZUnpt6dh/g0WAfCKl9GUApwD81lLDb6WUngBwCIs75r9bFEXdbSUHRpF4nS+nrCpgiBQsTss44MsxWQYabkOn7NwnXWONc2w0LtqzDBo5nd7UOXLRn8uyEdXLq7+r2o7fCrROn8D760RjY2MYGxvD1NQUzpw5g8HBQczPz2PTpk3o7e3FyMgIBgYGsGHDBrS3t6OnpwdtbW0oiqKcDs7Pz2N2drZ0tqtXr+LKlSu4cOECRkZGMDs7i4mJCVy5cgVnzpzBxo0bcf/996O3t3dZnzXY5VgbB1geY8cSOb/qXYFVHT0HZKxTDphurFWG3Nod95HZoQt2GkgcoOt/JSF8nm88UF07HagtO+DPlXe2r2klu94/gV93BIBHM2W+DuDr9eqWMuVvjozunDJHx66Urud2w3OAqQ4RdSgDc3JwfbpEoA6iEVbrdIakzsuRM8673VMGc42yTs96TPWe2ymsis6OZV+5cgWDg4OYnJzEtWvXMDIyUk65+/r6sGvXLmzfvr0ExKampvL2TpaDQTfaSilhx44d5TR+fHwcY2NjeOutt3D48GFcu3YNg4ODmJ6exr59+7B9+3a0t7dnWRXbDIOXPiVHgyAfq1p7Y92yPnUcGMTdOGj7rj9ujdnJHccjPzNflZX7wwCcIyI5MNPE9TqykeuHa88x3H8wUP5TJ2U3GrX0O/KwESvQ6PQjx/yifQWESAqSOXDj/JEvDEvL5Jhx5HfR0O02ah94OsFlFShVD+q4Luqzk+aicQ4sgzlOTExgcnISQ0NDOH36NEZHRzE9PY2ZmRls2LABt9xyC7Zu3Yqenp6auz3Y6dzaFc8+1F5aW1vR3NyM3t5ebN68GXv37sXPf/5zvPfeezh8+DBOnz6NO+64A/v27UNvb691OtYxAybrJ2TVQKhBu55zVjlxTr8KjBysc36j7akMrOMqG9IZl9anyzkKkiyTA1+VWQMjt6FldINQ3yhZL60aoGQjKoqixtgjVU314lvXYZjC8wDn1tAcEOr0McprfpciAuqaq0ZAHkyO1DkQd0xQDZDPaV41WK3HtatOxzrh9l0qigLXrl3D6OgopqamcOHCBQwNDZUbMjt37kRXVxe2bNmCvr4+dHZ2lmtVWq8yOn0njTobO2iAa3NzM3bv3o3NmzfjnXfewUsvvYTLly/jtddew+DgIHbt2oVdu3Zh/fr1y/Tq7ITtjmc9CjRRzrFK1VdujHMzHdZPtK9sU8eOZWR98f/IG9Nutj/Wc9irBlO1IdWdgpj2W4GXWaWOswIst6kXocc5t9ygadUAZaSVsBQHci6vY0O8iaKMzIGyA5wqwFtJ/zg/s1RlQlUyOWaqjurKMnDrzmNu44o3HDgxWOX6ubCwgImJCYyOjuLUqVM4f/485ufnMTk5ic7OTtx8880YGBjA5s2ba9i+gk0VoLAMuQ0i1h0DVEdHB2677Tbs2rULzz//PA4fPowTJ07g0qVLOHv2LG644QbcfPPNds1MZzrad13mWVhYqJm2Mmg63enYsR4ULFQf3FcFsipWx7JyWQYZBpeqIFK1ieP6rK94UB/jseV+xG9mibzUpkDM9Tm9ubRqgNKtnbnIEr+5DNfBF2crAOUcOgaVH7ihF8qqw+ogOZlUZgYBBm3Oy6Cn7M+dVxm4z+GYql8GZ3UQPlZviSLyBPvgdqempjAxMYGxsTGcP3++vJynubkZW7duRX9/P/r7+9HW1lazU8x9qNqNdOzZjatbM4zETKSnpwef/OQnsWfPHrz44os4deoURkdHcfbsWUxMTOCuu+5aBtZFUXvJjkvMrtyaXiTd9NEAUMXItP8c4NSOHKjqsdnZ2WVyOiYbbFXt2gXqqDu+nQ6UEbrNH5XbsWAODAsLiw/40NfIVG00ubRqgDJSTsnOOHTNkI1M2aLWrQPvpkEKjFqXAw4uw8ymymE/CDvNGTfLlGNSqpucPnQtKSdv1Md9m5ubw8jICM6dO4fBwUGMj49jamoK3d3duPnmm7F9+3asW7euxrmc/lV/2keWOzfDYPDW/jKQRGpubsbOnTuxfft2HD58GM8++ywGBwfx7LPPYm5uDnfccQdaWlqy46SBRwGKgVsds+opQ1GH073zE2ZTrm0XiKvAi+9hd8sIumaubJM32xwrdGPHoBYyst+rbnRDj591qa+JVp1+5ICSd5V1MCKxknNJ1wWVeVU5fXyzLJrX7Uo79seDG+XUWLhtZZhRP99Gpu3n+uGmKO6b63M6D2btdKltTk1NYWhoCMePH8eVK1cwNTWFvr4+7N27t9xN5rpy8iqgqbyqY+6nAyp9mrhO+1imeKjw7bffjs2bN+PJJ5/EhQsX8MILL2BsbAx33HEHNmzYsMy5Y0aS0vvrj5pU3zlAzAUMBSUHNPw7BzbRdpTPXZCt667a3xxZiN+OaIQP6G25obecPlR/PIvR/uZsIQIHA7kCbC6tCqDMRbRQJjuWskBdS+E0NzdXXlysEZSP8XqRKiw3UPwf8Nd9OTCPfuquNtfn5GDm5NaHHGBru1GG8zrw1LLuYQacZ2ZmBufOncPZs2cxNDSEiYkJdHZ24r777sOWLVvQ2tpa8y52ZqFu7BQcuG1mNFVgpHpVh2IQYebJOtm8eTMef/xxvPTSS3jzzTdx6NAhDA8P47bbbsMtt9ySdS4HhNyPeNYjA1FOBzxuuTW6yK8gmrM/1juDeoAP2zfX5Z4ZqjviymAVaCPwuutLdZy4XyFDMFw9zzpx9eY2pkJ+3bzStCqAEvDTXGVY6hi8zueMk4+5wVPgdPIoi+Jz/Jvb4La0Tm4rZ2CcmKGo4+g3O75uinBe/c0OmLvawIHX5OQkLl++jLNnz+LUqVNIKaGtrQ2bNm3CTTfdhN7e3prlEJVJGQ8nZrJs/CpD1UZWbjMv8uWYOTvNunXr8PDDD6OrqwvPP/88zpw5g/n5efT09GDz5s01G2MKOvw7gn4OGJ2uc0xRl0ZYj2onblmCba2pqanmyd+sE8e8HShVsTEFtygf7LsqUOrY8zLA3NxczRgz6FVtjnE+Hot/tIdi/KKTo+cMTjyAHGEc+9I1iVx043Y1Sirjy02LndHyVI4H2j2DkBP3KZy1KIqazSiuU1kDG52CvG4yqR50DEIP3E/Wz+joKA4fPowrV65gcnISzc3NuOmmm7Bt27bylsJgDjmAjjZ5jPm+YA1qCrTOAbgNLqs2w/kYABRgi6JAc3Mz7r77bvT19eHZZ5/FyZMn8dxzz+HWW2/FjTfeiNbW1mX6Z1kbGxvLh3fEOd1sY3l5PBVYncx8XAMH8P6T1t1tjjmQU9vjdhQgc+PhdBvfbhmJbUNfaBZlXABxYK1khr+1nCMomlYNUAK16zw59qagw9HX7dLlnCQXgSMxQPJ/BSAFE61DdwRZvtwAKQjmAFL7o33WJQXHMlkewF+8yzKOj4/j0qVLuHbtGi5duoS5uTns2LEDu3btQnd397IL+l0w0XZ16seGrTpQ/SmQMht0eXiare2pzbGOWlpasHfvXnR0dOA73/kOXn/9dVy4cAFjY2O4++670dzcXNMWM2K213grqAZOHVuW1QGBBn8d1zjH/XTjrmCZC54KSCyfyq+6jMS2FQw7t06rjC+Akpko69IlJ3t8Yurt5HRpVQClTjfYkIG8EeUUzclFuWgzUhW702NRjwKd5uWpF/eRZdH6tE5lpQx2KpMCSgCzm7Ir63QRP1JTUxOKYnE3e3R0FO+++y4uXryIyclJbNq0CbfcckvNRdmuvdwacuTn9jVIqp505zbaCOYT9WngjDJhN64eN64cEIuiwJYtW/Arv/Ir+N73vofh4WH8/Oc/x+TkJO68806sW7euhilqQNXx1+cpan7VIycNBq6v+jvyuKmpBitlcvWICB9T8IljwWzjt+rfLa1xHmbLoSe3bhptxjioXfMS1UrTqgBKoHaqos7v1qEADzaR4j/XG8f5vzNGzqvvodak7TBTVSalDuHkVpBzg+z6ovWxIatend5zU7vx8XFMTk7i0qVLGBsbw+XLl9HY2IhbbrkFO3fuRHd3d019Cphct7Ybhh5TbheYNEiqLmKtKqVUPggjHnixsLCAzs5OtLa2or29Ha2trWhoaKjZWFK96MaVrm+nlLBr1y48/vjj+OEPf4gLFy7g9ddfx/z8PB544IHyjiI3ttx3ZU2cdHYQdRVFUQYu1k2OLcY5tgG+1EfHmuULXXCgY2BSEsDLJcpsGcxya5Lcb2V+mlRel5RputkDt8WA6tKqAEqOGM6hc+DmACy3Fufy5Bij+3aApoxH69K+OMPX9bDI4/rh6ue63CaOAm3OQDRiT0xM4NSpUxgaGsLMzEzpMHv37kV/fz86Ojps/5zDOr1HHl5/rdIlO9Hc3BzGx8cxMTGBwcHB8rbI8fFxFEWBsbExXL16tZS7o6MD3d3d2LFjB/r6+rB///6aTUBlNaEX7VukxsZG7NixA7/2a7+GH/zgBzh58iReeeUVbN68Gfv27aspw4HIkQAFBh7PSLz5oJceqe6Y4WmQ4SCea4t/88aeY5o6lmqDDNQcNHmJRNtm0GTmp31MKZVBUsmHAm7k0/4BtXfuVKVVAZRMxVeSNCrlnNW9QsE5Mhs08H7ki48DRWfMumPowEzBhetwjJZl5jVcLcv5mNWpoarMLM/8/DxmZmYwMjKCd999F0ePHsXExAQ2b96M22+/Hd3d3ejr6ysvuXIsj+VQ4FOWzMdZHzqe4TSjo6MYHh7G8PAwrl69irNnz2Jubq581uS6devQ39+PrVu3oqOjo3zy+fnz53Hs2DG888472L17dwmazCxV3zyG0Vd12P7+fnz+85/H008/jffeew9vvPEGUkrYu3cvmpubS7nd0k60lwMK1mM8cJqBIweWLL/TJfdT2+K1wNhwChbpgDhk0E1PlsttGOl4s19U+RnLpstFfB96tBEgqLMlPe90r2lVASULm1PcShIr1EVubVuZUA4MI7kpsZ4HUN465RbdlUko8Ki8zH5y7XH/3QaW9o37yM9nPHfuHEZHR9Ha2oquri7s27cPN9xwwzLZnH6VdbBT8SuEc33RIDI2NoaRkZHy2ZTnz59HURRob29HX18fbrrpJvT19aGpqal8OpDqtigKjI6O4s0338RLL72Ep556Crfffju2b9+O3t5edHR0LBtLLhuAoWviKSX09fXh05/+NH74wx/i4sWLOHz4MGZnZ3HLLbegtbW1dEjXNx0z/s2MN8oq++YgqGDC46LTap2CatvM3PWYs8Oc3oDa932zjGGXMTVXe+Q8zEbZprQtBnvWhV5twMdVP7m0KoAS8FMRZViRHEtiVsrGrYPDv3VwVI44xzIqyHFeVXgYlG6quCkyt6W6UDlUJ/WcjuXhcnNzc5iYmMDU1BQuXryIEydOYHx8HAAwMDCA/fv3o6OjA11dXXX7r2OmzqkG6YyTdT43N4fh4eHy8WdjY2Po6urCnj17sHv3bqxfv768T1ztwwWGlpYWfOITn8Du3bvx/PPP4yc/+Qk2bdqEvXv34sCBA+UTzqtsj5kJg0VPTw8eeughvP766xgcHMTbb7+N9vZ27N692z5nIDY1nB06cOCx1LVATrm1zfjmGzh0Q8PZjPOFHGHJzXKiLX3gb7TN4Ot2sOO8bl6x/TFz1wDJYMybbFwXgH/4GmVKqQ3AcwBal/J/uyiKf5tSWg/gPwDYBeAEgN8uiuLKUpmvAfgygHkA/7ooiifrtSNtlgrQy3Pitw6aGgLXxdEQqL1VMv6z0TK7yxkeR/tIbu3TRT11SHUON9gsZ+gmxxYde2Y9Tk9PY3R0FNeuXcOFCxdw5swZjIyMoLe3F1u2bMGWLVuwadMmdHV11dwZpbpwifvEOubzqiOW8dq1azh9+jSuXbuGo0ePYmZmBlu3bsXdd9+NHTt2oLm5ucyr93DnZGGH2b59O77whS/gxz/+MZ588klMT09jYWEB99xzT81GSbTBY+Iu9YrjGzduxIMPPojjx4/jZz/7GV599VVMT0/jlltuyQboGMecHuObAy07vdO52lv4hauflxIcy2RQccsPucDIdaic/FtlVGbIa7O8Y6/5uO+6RMFr8FHePXi4yqaBlTHKaQCPFEUxlhZfMvaTlNLfAPgCgGeKovj9lNJXAXwVwO+llPYB+CKA/Vh8Z84PU0p7izqvg1CjBmrXQNylAprcupmLHg48eeof+ViZWp8yTccQc/1jx3VJ69A2XZDQ8trfYGmTk5M4deoUpqencfbsWYyMjGB6ehq9vb2488470d/fX95y6F6bofI44GP9VemOy8zPz+PKlSt47bXX8Prrr6O9vR39/f34xCc+gYGBgXI9MRcAVGcaaDgILiwsvpxs27ZtuHjxIt566y1cf/312LRpU42sCrI5+SNPS0sLdu/ejfHxcbz11lt47bXX0NjYiOuuu65kvwxKCpzcngKY5tWdbse8eNmDgUb146bwbty4XgVmN8YaEDiffhRoGRAd680tO6g8YfdaB4PpPwpQFos1jC39bV76FAAeB/Dw0vFvAngWwO8tHf9WURTTAI6nlI4CuBfAz+q0swxMdJoQSZ3SsZMcYDFwxkefY8eDzazIAZCm3IMrlGFF/1QHKzFWdwmGGofKd+XKFVy9ehVDQ0MYHBzE3Nwctm3bht27dwNYfGXwunXrkFKqebUrs9/cFQXcrgKMgleUYXmnpqZw6NAhnD9/HufPn0dLSwt++Zd/GQMDAzWMlkFLndONNQdOXSu99957MT4+jsOHD+Pdd9/FoUOH0N3dXT4dSFmjAwxnSwBw5513Yv369fjJT36C559/HufPn8f+/fvR39+/jK05AHZtFsX7jz/TsXVjz+eiHXeJDB9T8qDj5UBI24ljeseR5p2bm6tho7ynECCZu62Qx8fVrazS6Sn0ESDqdMdppa+rbQTwCoAbAfw/i6J4IaW0uSiKwaWGB1NKm5aybwPwPBU/s3RsRYmdTdci6yXddSb5rQExY+XjUYYH0gEfp8hTb0dejVTzuPajfm1TL21RXYXBjI+P4+233y6f6NPZ2Yk9e/Zg27ZtNWtovNOpbbngVMUQ47guaygrmpqawosvvoi3334bDQ0NuPnmm3HzzTdjw4YNpUxuJzX6r2wgApAbcwa2bdu2oSgWn6y+adMm/PSnP8XCwgLuvfdetLe3l7qLwKHLJerUOjXetWsX2tvb8dJLL+HYsWOYmJjAgQMHMDAwUPZJnZPLu91YtUl9FqYCANsJAwLXpXpSlqjnWS7H5jhP/Ga75PVEBijHtnUtVYkGf5iBK/Bx2TjOz9wM4P4HA2WxOG2+I6XUC+A7KaVbK7I7VFsmQUrpKwC+AixuHETnl87VfLPCVuLAvDvG00euQxmKka/md24XmdvW3U3dxNFNDT6mhpdjrvzbPbQgysWbBuOp4hMTE+ju7satt95aPsyB+66gpKzKjQsfV0YZdbCTcL7p6WkMDQ3h6NGjuHTpEg4cOIBt27aV4K2OHIAVMrknV3MbKofKGI7Z0tJSPuXoRz/6EVpaWnDvvfcu04Xan7JCBSEA2LRpEx566CH8p//0n3DkyJHSYbdt27ZsM4ZtRnXngnKUCVCMHXkFLwYOBhoHeAo8epwvF1JZWPd6PMq7YM95mDHzcX2uJo85M1gGS2buLBtfMuQ2GXPpA+16F0VxNaX0LIDPABhKKQ0ssckBABeWsp0BsIOKbQdwztT1DQDfAID9+/cX6lw8LXVrLRUyZhmROh7X73bJoxwPfhVQMkvlSOgAKPqjUV1/5wbSOX84xPj4eHmZz/DwMPr7+3HXXXdh06ZN5eUzDszCkGJTQw3TMWvVr7IR7fPc3BwuXryII0eO4OTJk+jt7cWDDz6IzZs3l3U4EHGXWDl7yI2/0yk7ys6dO/H5z38eTz75JK5evYq77roLra2t6OvrWxYoHeNytgMA3d3deOSRR9De3o6jR4+iKAoMDQ1hz5495VKH1qG/dby1zxE4nE54DFSHysajfLBRDUTMUvm/7g3omGs/FMy4TwzMqsuQl9l49D3yBvsPhsjT8Kgjpv38RKMqNgmsbNe7H8DsEki2A/gUgP8OwHcBfAmL7/f+EoC/XCryXQB/llL6Ayxu5uwB8GK9drhT4Yy5i9AV6CKvM2IeQN3lchsBeowVmAMJ1xbLxXmc7Cqb66Mrx8YUGzVDQ0M4efIkhoeH0dbWhjvvvBM7d+60F4lr31VGB3ScdP3WBaD4PzMzg0uXLuHMmTM4ceIEOjo6cN9992HXrl3lTnboi4GM6+LpWbTJbC63eZTTN4957Fw/9thj+Nu//Vu89tpraG9vxx133IHe3t4a5q/Oy8eDFXGejo4O/NIv/RJaWlrwyiuv4OLFi2hvb0dnZ+ey53RG/xRoHIjqWObW8rkuvQxGwYzzsp6VCKgO4rh72ITqTEGTmZ8LFFyW76WPD4//7OxsjQ61PPfZ3cOeSythlAMAvpkW1ykbADxRFMX3Uko/A/BESunLAE4B+K2ljr6VUnoCwCEAcwB+t6iz4x1JL6/hOwQ4MfDxIGn0ZYYS5XRKoEYXhq+OoYZZb5PJtcEA7i71qeoH/486OKJevHgRR48exfDwMBYWFnD99ddj7969WLdu3TId61UErCuOrrn1xYjILLPKFzLOzc3h/PnzeO+99zA0NISOjo6a92fHxhH3O6fL0JXqXPXjAhozEc7PU9aFhQX09/fjV3/1V/Hcc8/h4sWLuHjxYvlUJDd15fp1qYUduLW1FR/72MewsLBQviK3paUF/f39JbPkmVO0wZfoOABVUuDGQO1MWZz6nTI+JjBsmyynriVqvQq6To4ow3atvqdtMcDGcWW+rDu2bw1MVaxyJbvebwC40xwfBvBopszXAXy9Xt2m3LJFc8BvpmjSe2CjPvefnzSjTJDbzEXuOA8sv02QF/l1HZDrD1niFjxg8aJozsP5FAzm5xffZDg5OYmLFy/i7NmzGB4exsDAQHm3SktLyzKn0d+6bqvOpk7KBqUgy0A7Pz9fvljs7bffRmtrK2655Rbs2rULbW1tNUDLzuwCEjusS1Xj5PJEO7lA0NbWhnvvvRdPPfUUXnzxRXR0dGDbtm3LWJ0DRNZP2Eb0s729HQ888ADa29vx6quvYnx8HJs3b8add95Zvkdc62fZNMDpx9lXjlU5m3KBmAMAsHxdPJK2l2OEDLjMbnVXXkFUgVJnEbxurYDPzNMFC5XFpVVxZ45jTrpOFYlZZEq171OOtUbtMA+0OpMCJZfNyRDnom7973aPnTFy+3Nzc5idnS3XB3WazFF0enoaY2Nj5S19165dQ1dXFz72sY9hYGAAzc3NaGpqWuZs0R7r2K39KYtWcOTyqt/5+XlcuHAB7777bnmpz0033YTdu3fXrI9q/zQ4qN5jrHNPk9FxjDHJre/yN/c7fvf29uKxxx7D008/jR//+Me46667sGfPnhowVIdmMMsBU1NTE2699VZ0dHTg2LFjOHnyJNrb23HDDTdg48aNpS51nZGdPHbG9RIcHuuVAKTqgEFG8zqQCVnUZzhv1BEyByipzh1Ict91PZHzx/FIzETddaN8joGyKq0KoHQGHscjKT3nsuy0AGqmSTmw4wilIKqOlJOvam1Dy0ZbmhoaGkpgU4BgoAIW2efVq1dx7Nix8qGxLS0tuOGGG3DjjTeiq6sr22f3nx/ckWO9jkG6fhZFgcnJSbz77rs4ePAgUlp8HNm+fftKuRQUc+zVtcXRP/6rzCobO5UD0DiuthDOv379evzar/0ann/+eTz99NMoigJ79uwp8/POuz5ph/Wq09HW1lbs27cPmzZtwt/8zd/g0KFDuHbtGvbt24eNGzeWT4h3lwlpYObHzCkw6DjW0w0DoI4v/2dwytksMz59GHTog+8l53q5324qz+uguaDEQMiyRp+jbZUhl1YFULrk2BewnOUpK+B8uaTgyAw1juWcOcq75MCG//MLlaqYrfavKBY3CWKj5vTp05ibmyuf6hM7s6wfZ7zanr6mwRkKGyI7P7OCkZERnD9/HqdPn8bIyAh27dqFG264Ab29vfblcI7hqgzKpGMtUx+goCBbtXbsxoSngcrIGhoa0NHRgYceegjz8/N49tlnURQFdu/ejZaWFvsgF12nrJodbdiwAQ899BBeeOEFnD59GpcvX8Z1112H/fv3o7Ozs6afPDbuwnM+7/yGjysDVDaVW8ZiPSoAsX0pKDGoRV26Lhhy8VPLeWxCLtYHj7MCJrNOt5apMlSRHmAVAaUb8FxyhqAdzb2rRY3FRVMHKjk54jzvuvO3MhaOng4oo16+dGViYgLnz5/HW2+9hYsXL6Knpwd33HEHduzYUTqsmwYxa8xtmFSxWF17VX0MDg5iaGgIw8PDKIoC27Ztw1133YWurq5lLF0dz+lUNwz4GOvHTZOYqcRvfQACt81TZHZgtcMYi0984hNobGzEd77zHTzyyCO46667ll3REHI4IFdwinrjioQXXngBw8PD+MlPfoKxsTEcOHCg3G1nR+e6XJDjpQHXLw18zA51LKK+XFuaFJxywMk6Uzn5E4yQp+1RLnTCuuZ1SgeaLJ8C5UdijRLIX86hjs2RGli+0cP3BHMZbkMVwnU6xqNgwQvMmnKgymUcoKmM8/PzmJ2dxYULF/DOO+/g8uXLABbfN33zzTejra0NwPu72FW6iz7oRpl+q6wqG7A4Vbly5QpOnjyJS5cuobGxEdu2bcP27dvR2dn5gZ4r6nSnIOPWQR0LYvBzU0LXVkNDQ80OKQcz1WNzczMeeughtLe346c//SkaGhpw++23W7bIsisb1HFuaGjAwMAAHn74Ybz00ksYGRnBa6+9hhMnTuCRRx7Bpk2b0NTUVDMVZzvVNnTqraxU5dHk1hhD/zwGju0y4KiPKhArkDFQ6fpr1BF5+D1AwX4ZcBUU43humu3a07QqgbLqOO8m80DEGp8rz0nXpSI5ZhnJbdxwNGTjccAT9TswdgGiKBafwxgbIlevXkVvby/uvvvu8sENuaUFBfrY1NG22HDdBpjKNTc3Vz5p6MqVK2hqasL111+PLVu2oL29fZlOc4HPjQ3nc2yp6r+TXfOr3nW9W/8r6+dnGN53332Ym5vDs88+i9bWVuzdu3cZW6+aHUX7aic9PT34+Mc/jr6+Pvz85z/H+fPn8cwzz2D//v3lS81cvxgIcuvyKwFrRyA4uDObZz0yG+SkDJJ1yKxf7VKvjeQxUkYY7fBDL2L9UdmmbgZxXfXYJLCKgBLI73bm8gK1a3JsfM5hXVKD4qmGOkyO6brLTBRIXTssW+QNQDp9+jTOnj2LhYUFHDhwADfeeGP5zEReKM/pS51VmYe2q1Msdp4rV67g2LFjOH36NPr7+3HDDTegv7+/vA6S2Y06XNX41WPkyio58bjzGOmuc66N6LfWz4GEHR1439kfeOABLCws4Ec/+hGKosBNN9207OlTPDZcD/ebr+AAgObmZuzfvx+9vb146aWXcPToUZw7dw6vv/467rvvvvLhyRqcuS/K1PQ/BwoX1HUdMG6PZPvOTW0diEc5BirdZFGSonXyKx/iE4DtAHB+fr7m3m0eS2aOfAeSCxScVhVQRuIozxsg6jA5VqXGGGklgMkD59iFskt95WjUpcbLAKbMqaFh8Y6CK1eu4PTp0zh16hSuXr2Knp4e3H///di0aZPVA9fNIK8AnGO5bBgMvgAwOTlZXgd59uxZNDY2Ys+ePbjuuuvKXVldl1V9Vul7JSnHCrmf6mhalh1Pp5KRj3Xk2JfT53333Yf29nY8+eSTuHTpEu6//360tLTUAADwvv3FRhQ/JV3ZfTzabseOHVi/fj16e3vxxhtv4K233sL58+fxyU9+Erfeemt5vS2DZI4x6/qf6kZZl+orjqsvuPKst9x94Tytj4+7QJzbZTl0eu3KOL3of12n/EhNvcNRdT0jpWTXHuNbo5tjeEA+ciibUjCLevmYrtVwufjtjEvLLCwsYGxsDOfOncOZM2dw9uxZFEWB66+/Hrfffjs6OzuXMTYFIGWtrEtdl3Sgos4STxafnJxER0cHbrrpJmzcuLFktKwr/l0FMJwY8PW8u/SKgYT7p6mKFajeGBijDS2rjJWBuq2tDXfffTcmJibwzDPPoKWlBffcc499Vzf/5g0mbk/XrePi9GCXFy5cwE9+8hNcu3YNDz74YM3NBLzZoWPt7Fj1ESDR2NhYTmOVefKUW/Wl4OPGwoGyXnnCjI/rVP1xP5RJqlzcP5aFL0CvB5LAKgLKqqQOqExNI3OuvAKITtuccbExxifa0sf8q6OxnI7tLCwsYGZmBocPHy6f5r1lyxbs27cPGzZsWPa09ijr5FO2GfXzDmsc0+ABLBrTxMQETp8+jYsXL6Kvrw+7d+/Ghg0bSgaT27iIT6yHxnEHkg6Q9LybOeQ2GHJyRRnHCNVptT+qG74CQdfk7rvvPly9ehXPPPMMtmzZgh07dmSnoMpugNqAw+uts7OzaG5uLq+3fPbZZ/Hee+/hzTffBADcdttt2LhxY40tOPt3Swjaz+ib3tkSAMybKK4dJg/RltbF64hKIrhOB+TMRLlPzDR5qs1grw/HyAXZena5qoDSCRrTFQUHdUQ39dPfVaDF/3V6UG+Kn2Mybooc5YpiccPm8OHDOHnyJBoaGnDTTTdh//795dOw3TqctslG6uRwa0EsQzxIY3BwEFevXkVbWxuuv/56bN26ddlrINz6LbP3ePKQ0+tKkmOT8V+XUJy+GZQY5BSodIqqgBF1cd95aYLZW1NTEz7zmc9gYmICf/M3f4PPfe5zGBgYsHbhGC0/wJb7ELI3NTVh06ZN+OxnP4unn34aJ0+exHvvvQcAOHDgANavX19jX0oK+BzrVfNWLUXoOHBdbpbG+oz6NfjwWKkMOqYAagCbQY/zMOjymKutMIBqO7m0KoCSnV2P6YDyIPBA8HpGJDdgPMVRJ+C23NOIgOWXzjjGpwaqzHV2dracZsctbPfcc0/NcyJV9lx7mrg9Bvvof3zm5+dx/vz5coe1ra0N+/btw549e9DR0VHDZvXyI+6L6tkZJS+nqCwuUFWtSwO1j/eql6I+fV4h24KCiwYofoSZBrxIn/3sZ/H9738f3//+9/HYY49h165dy9gk20Mk3XxSxh+yrFu3Do8++ih+9KMf4a233sLMzAyGh4exbdu28vmZaqvsVwwM6lN88wDLFMc0gHD9CoCcJ+RxTE7tgqfCUR9faK5riioLBx5ms8pCeU3yIweUHOWd01U5Yhyv2gHW/7xT65xVy3BErsdOnTMA7xvN7OwsTp8+jUOHDuHKlSvo6+vDfffdhw0bNgCovS5SmVvUo+DHbWge1uHc3BxmZmYwPj6Od955p3yY7P79+3HLLbfU3GrIOnXgyEmn9yqP/tf1OM2jDFbPu+AU9cYYuTXknGwcVHQN1rESBVcAaGtrw6/+6q/iu9/9Lv7iL/4Cv/7rv44dO/ixrH5DSgGGdcz3thdFgZ6eHtx77724fPkyhoeHy2Db0NCA2267Dc3NzZbZcZ1BKHLLVC4wc8Bie2TQiv+6u81g6TaWQiZmgfE7AiIzSQ06nJ+n/O6SH87DgYn1n0urAig1scED+Q0TPpaL2JGYXbl6c/XF/8jPbeh0VjdU9PfMzAx+/vOf48SJE5ifn8e+fftw4403oru7215CpP1VOVz/XH8XFhYwOzuLEydO4NixYzh+/DiampqwY8cO7N27F9u3by8valYwVHashso6c3oE3r/2TvOp7lj3CwsLNQ8QjnPhiG5TSceBncExIt18cOu4DMIMqFFO15E//elP48///M/xV3/1V/j85z+P6667zm7guOsRdTbAKVjQhg0b8Mgjj+Dpp5/GtWvXMDc3h6NHj6KpqQn79+9ftpnE9em0mfPFtwOhHIt0xCV8113So20pm+TxCmaoMursj0GV69Yll0jMUnns6zHK+nOXf6LERqmDk5u2KYPic5pP15u4Pk455TlFVrEWBpWpqSlcuXIFb7/9No4dO4axsTHs2bMHN910E9rb25eV0/LaXzUsZX1qlLOzs3j33Xfx8ssv49ixY9ixYwc++9nP4uGHH8aOHTuWvUxMQVudWJ3F6SiMmBfTHcsOPSqDrZouhkzOiWNMlKXk5FR9u0V9ZsBVrDmlhJ6eHvz6r/86AOCv//qvce7cuWUBmp2bL47WtVS206ampnL9d/PmzXj44YexdetWtLS04NKlS3jhhRdw6NChUia2ldwF1cz2WF/8X8lCyBUbKJqX73LJAa1eIM5tzs7OLmN8ygRZT8w8I79bz4ync7FtaP8/coxSDZuTGjxHfl3jZIDk427NLQxAzwH+DoWVyL6wsIDJyUm89dZbGBoawqVLl9DR0YFPfOIT2Lp1a/nUoEg6JarqqwMQLb+wsPha1oMHD+LkyZNobW3FI488UoJjXAvJMrAsCmQqE/dXZeWprJtCKxBrHnakkMGlhYWFGhbFZVg+dVoXdHV5IxcUuKxeBdHQ0ID+/n78zu/8Dp544gn8xV/8BT75yU9iz549ls3mpn5qz2yDjY2N2LlzJzo7O/Hcc8/h9OnTmJ+fx8GDB9He3l4+A0DZlGNwLjC52Qn3W+Xl6byuIWo7vH7o/EsDvVtL5OO6pqnTfZU3AF7rcwGS06oBSnYsnbKys6kxsaJ4+sjl3VSSB0XX2NSZXLuaNwYo5J+fn8fU1BSOHj2KQ4cOYXZ2tnwyzLp162re7Bf15aJ+yKrrRJwY8CcnJ3H16lUMDg7i8OHDWLduHR588EF0dHRg69aty6awOT1pvTp91G/dmc6V0cAWDpZbk1SHcsk5RsihAKlOodNqB4zK1ll2LbOwsIANGzbgC1/4Av74j/8Y3/rWt/Dbv/3buOGGG2ra1KlhlZ3Geb5gfePGjXjggQdw6dIlTE5O4vLly/jRj36EW265pbwGV/2BmaObtbj+OrtkWRUkXRvA+wCtlxullMpXODCT1Pb1NwO7AmC0peu8rF/WcxVIAh8AKNPiqyBeBnC2KIrPpZTWA/gPAHYBOAHgt4uiuLKU92sAvgxgHsC/LoriyZW0oUbHTE9ZgU7VojxQOyXmy4vY8NhZ3I656X9NW7zupClAMt4XnVLC7bffjt27d6Ojo6ME/pCTn3TEjqxgzKChwB3tjoyM4PDhwzhx4gQ6Oztx++23l69OVX1xG67OHINm8ORAw7JoOQ5qkXL/ua5/CEjm5FawY71yPleGHTnq0HGJMhs3bsS/+lf/Ck888QS+853v4NFHH8Wdd95ZlndLTBHEHBPXAF0UBQYGBvCpT30Kf/u3f4tLly7h6tWrGBkZwcjICP7ZP/tnaG9vr2FhAJZdbhdAo/XztDr6yrbGiV+Hq+uIfF7zMLBqvQpsPA1nW9Hx4fNatzLRlNKy9VCXPsga5b8B8Db9/yqAZ4qi2APgmaX/SCntA/BFAPux+LbGP1wC2cqkxhqdAGojfhxX9qLTFB1cjdZ6SYQq3gFFjnVFOwsLi5smo6OjOHLkCN5++23Mz8/j5ptvxp49e9DV1VXephb1uGdC6rRQ9cR9juv6zp8/jxdeeAE/+MEPMDg4iLvvvhuf+tSncMstt5TgzP3V9nJTZAUMPhZTd9bz3zepU6gz5cooSLKjODaiAUinoW4KxrMWtQ+VVdOGDRvwL/7Fv8DmzZvxve99DwcPHrQBSYN2tOcCgM66rrvuOhw4cKAEu9nZWRw7dgynTp2qedlZjlHpbYaO2cU6n06FFbQYEJ0fcpuxRsssUz/R7szMTA2DjLIaBHS2EudZdvaxqjVcTitilCml7QB+FYvvwfmvlg4/DuDhpd/fBPAsgN9bOv6toiimARxPKR0FcC+An9VpY1kkVwB0j/DiHc4qp4o6c0wBqDVSBaRI7CScwqAuXLiA48eP49SpU2hsbMRtt92GXbt21UwzNABEvW5XUtkbl5+bm8PIyAjeeecdnD59unwdxLZt20pA5kDBQcEx92iP++SCg5bLMU8dC7er7sZG5XDtKlg5INfjPOa5oMdjomtzbormwIfLFkWB7u5u/MZv/Ab+6q/+Ct///vdx7do13H333eWrMRSUdJdc+6XHFxYWcNNNN+H8+fN4/fXX0dnZiXXr1qG5ubkESh1/XTvkqTLLElNXDRbqcwxs7EeqC66LgTXO8TS+3gxE9a6skwGU1++dndSbkax06v0/APivAXTTsc1FUQwuNTSYUtq0dHwbgOcp35mlYzUppfQVAF8BgIGBgVJgx+LiXA6kqs4Bfkqu7CIGntfKcnkUuCLFmuTZs2fR1NSEO++8s+YOF+4Py5V7E2HU767tm5iYwKlTp3DkyBGklEpAjgdWRH6nEzdNjDtE2MB1GsY6YBlXoi/tWxX75GChMjkd8tiF3K5drjPHHLWsAm3kU6DWe5dZNmDxOsvPfe5z+NGPfoSnn34aV69exUMPPYTW1tayjzENVLLAsnP7rI+Ojg488MADGB0dxaVLl8p6eE2TxyjKK6PmfKpzLct6VWALgGYiw2DL64buAnFtQ+2OL0PiumZnZ5eBPsuoYM36rCJZK3mv9+cAXCiK4pWU0sP18gNwHrBMgqIovgHgGwCwf//+wgFPTCs55XaxRWYLZjmHdcCrdXH7mgcAJiYm8M477+DcuXPlK1l37txZyscOyM7PUxvHWFSOhYUFnDt3DgcPHsSlS5dw44034s4770Rra2tZZw7YtT7XX9aJRlmnJ8eKNQ+vfymjVvYcZap0wSm3AafBlsHRAY2CqbMRLqe60L6og6eU0NbWhocffhgdHR149tlnkdL7TyHSpSTWQ9TvAJwBu7e3F48++ih++MMf4vLly3jjjTfQ19dXvlZC7Zg3VwDUsMdol6e0qjteU+W6c6xTy3IfuX4NSGFDfD6Wubgut0SisjLr1f9VaSWM8kEAv5ZS+hUAbQB6Ukr/M4ChlNLAEpscAHBhKf8ZAHxLwnYA5+o14piDo8PqCDkW4Xa/nfFGu8qYHBvTc1H32NgYDh48iDNnzqC7uxt79+4t7/fNvS4gEretAMfl5ubmcOnSJRw7dgxnz57F+vXr8eijj2LDhg3lbr/K66Iyn2NdOkbh2KMe540HFwCK4v2HZSiYuMSXqbADOQaaA3sHsOqYqhd2MrdeqzqsB+IOcIHFZ04+8MADAICnnnoKU1NTNZsu2kYuaChIxlpbd3c3PvnJT+KFF17A+fPn8eMf/xj3339/+X53rTtYJ9uC2qSu3zOo6DWTzCzdumHk4yk21+fkq1pvrvowRuhUnAOfI1CaVvJe768B+NpS5Q8D+D8VRfEvU0r/NwBfAvD7S99/uVTkuwD+LKX0BwC2AtgD4MV67YSgDkyis7pWqEbPjuUcx0VtPh+p6hIjjZJzc3N49913cerUKWzcuBG33HJL+eQf7UuOPcV/dlguMzY2hvfeew9Hjhwpb3mM+8JVP87JXLRWnbu1TBdguDxHZdVvMJYc0DH71WlQjmlqcoadY5YhA09FlcW7IKhOlNOhgklO7rCL+++/H1NTU/jhD3+IlBIeffTR8k4klkvHkXXiQCWlhL6+Pnz84x/HX/3VX+HVV1/F8PAwPvGJT2D79u01AKF2zeNRxejZThU8Y8yZvdfbnGL9aj5+fQMDKq9txn8G1dx0mkEz+lIv6AH/sOsofx/AEymlLwM4BeC3ljr0VkrpCQCHAMwB+N2iKOo/8A1+OqOddYYd3+pcGsHYKNxAVsnknH16ehonT57E2bNnsXnzZuzfvx9dXV01myhVdQPLn6bCxjs6OoozZ87g5MmTpXMNDAwsi/gMdDnm5Jw7dzz6rUFH+6Jrk7rWpXXror7qmJPrS05WBSad+nFfnPxRVsEwB1I5EOfxc0DHfWloaMBDDz2E2dlZ/PCHP8T8/Dw++clPor29HXNzc+UVBQ6EFZwU+BYWFtDT04NHHnkETz75JI4dO4b29na0traWj2Zz5XTstX4FMx53lU/XCZlt8u4zjzWDmC6P8JRbAVHXPrkeN24sZxyfmZmpnO18IKAsiuJZLO5uoyiKYQCPZvJ9HYs75CtO3CG9vq0o3t8RVoNgEFQQ5em0a8OBSpVsPPBzc3MYHBzEkSNH0NHRgb1795YgyQDJO5hq+OyoU1NTGB8fR1tbG5qamsq6R0dHsW3bNtx6663lBcTBjJTl8u+q/migiXJ8zakGJ10WyIFPjqm7qK2OyrJwoNP+5fqgQOX6kdMT16fv4eFHoYVMAWbstLoOrYnzxxg+9NBDmJubw9NPP42ZmRl85jOfQXNzM4DamxgYHONYsGNdy4+x6u/vx6c+9Sk8+eSTOH78OBoaGvDxj38cvb29Nb4RY84bIwpaPEUO3eqjyhioQpcKcO42xMgb592yRZTJ3X7I8ukY8BqqgqReJphLq+LOnHAKXScC/C6wGj/nc9Mmpv6qEI78DMbcptY9NzeHY8eO4e23Fy8rvf7667Fx48YawNYpB0foSDHgU1NTePPNN3Hq1Cns2LEDc3NzOH/+PFpbW3HgwAFs3769fD+NRn0NAo5J5fSooMYPZs0BIv/X8jljqwI5bVuZSw74c/0G/P3YOYDNMXBlRpF0TU5ZteuDkyVAIZjl+Pg4fvrTn2LXrl24/fbblwEV903ZVhwPQA0Zm5qasG3bNjz22GP467/+axw6dAjt7e2444470N3dXdNfXf5QIIopMBMUzu/KRd16naYCKSfHBouiWHY9KDPThYWFZaDNcrn2dez+MTZz/kkSKyw6z9cexnFN7DAcedU5uW4HIPUYGJ+PB+6Ojo7irrvuKjduFGgVXFx9MzMzOH78OI4fP46RkRFMTEygubkZGzduxB133IHe3l40NzfXXGKkEdMxO21HQSCcjzdiop7clEx1zvXlrpGsAmgHbPXGixMDiDp9FShzH6Msr0vHMQUkZtXuagxlJ255gesN0GtubsYv//Iv4/z583jqqafQ09OD6667ruxjlGN5GZQ4OKvOAWDHjh34zGc+g6eeego/+clPMDIygl/6pV9Cb29vDfAxc+S+hy/yfwY8p6tgu45xcpssN4+P6pdBNYBRgVl1pYDNgSfGR0lSLq2apwcBngUqY2Cj4LysHMc4OVLl2lYZHDOYmZnBe++9h/Pnz6Ovrw/bt29fdiG8Tv8cUwFQbgQdPHgQwOLAzczM4Oabb8Ydd9yBjRs31oBkVcBgPYQM7PxV7FLXfaIs53Xlq1gk64K/ta6VlOexrWpDk4KiC4qOMUdS5wm7U4YSbXF5B9COvUagamtrw2/+5m+ivb0d3//+9zE4OFhTX9X0UVma9qcoFm91fOSRR7Bx40a8/fbb+Lu/+zsMDQ2VZWPanbu+0TFu1Y/TF/ed7Z7HRFkhM0a95lFtN+RkPblbErmfavPahkurCigjuWmWMpsccLBBunVNx4bYGHShOVIo9PLlyzh58iTa2tqwf//+8mngWq9+60DMz89jaGgIr7/+Os6fP4/R0VF0dnbiwIED2LNnD/r6+soL0cOZmP1xqmKOCnAppWXsVHWu+tS2ct8uyGjQyjmScwA15KoApnUycChYOidVOXnM+GEr7j1GWhcHZT7GtsYpxrWvrw+/8Ru/AQD49re/jXfffTerq5yuq3S1detWfPrTn0Zvby/eeecdPPPMMzh//vyywOdmKqyTsCElHwyePHZxjm8XZHkD3LQeDgR8jgGdp996n3ic09sU9Vi9YA2sIqB01wHqoPHDLYDlRuMck5MrG21zcsAGLEbCd955B1euXMGNN96IzZs317Sjcii4x+/5+cWHV/z0pz/FwYMHMTs7ixtvvBEf//jHceDAgXI9UmXSNtxakQMcZ9CqD464Ou1Sp+Q2nU7VCB2wOwDmdSftc3y7+3JZJicPXxepjsbBSG2Iy/HvlFJ5iyjrktd4VV98XzPLr0Fww4YN+Of//J9jenoaf/Inf4LXX3992Rog95OZVzArp4OoY+vWrfjlX/5ltLW14Z133sHf/d3fYWRkpKybx82NE7PP+O82cKK/fIyn9vGJ509Gn3izjGWP+71nZmbKJw3pTEjbyh1j+1zJfd7AKgJKwANWjkWoEbOB6LobG3hu+srgFnUq6F64cAHnzp0rLyqPabFrB/CXI8zNzWF8fByvvvoq3njjDWzduhWPPvoobrvtNuzYsaO8w8bJygbJIMy6UYB305Ycu1bZA2Q4Dxtzrp7ctNoxIT2n9cQ3T8OUqXBe1k+Moy7NcNBVfXAdqhPuP//PTY+1Hg4Gcasd1x/fGzZswO/8zu+gt7cX3/72t/Hmm29idnZ22XucNNA7eeMTAaGpqQmbN2/GY489ho6ODrz99tt46aWXMDIyUo636poBSQFRx4TvrearPFxdbnrM0/0IBHycQdatbWo7vBGlOuEAVC+tCqDkAXevxFTAYsNSJqCGrMbD9XB9DLzRPjvb9PQ03n33XYyOjmLHjh3o7u4u22An03Z44GZmZjA6OoqXX34Zr732Gvbs2YMvfOELOHDgQM2mTcjjpgTsKKEvDRSsMwUINhLnbE5HrMsoG0asET3yujFguVz+qiUW97I3B9S5pRV1BmcnuaDsmIlOr3kKGMe5L6EzHVcFjLDnvr4+fPGLX8TGjRvxzDPP4OzZs3YcGJx07Kv8Ztu2beV1nM8//zx++tOfYmxsrEZPLKuuISqxUAAPneishfXIQBZPB1Kmqq+5VXl4c0gDzuzsbAmSmt/ZSVVaFUAJLKf0fGmNdoTBrYod5ZKyjSijDsgOPDg4iBMnTmBgYAA333zzsnbUEPn8wsJCySDee+89vPbaa+jt7cXDDz+MgYGBZReou11gdjy3TMHOqXpgp1JQYBbFunR1q0wc8VXnjjVwW7pWx86gY8LTQbc2zTbCTpgDjpCDmYkDZw0oym6iLt3c4UeAhQwO6FlHDEKh/3Xr1uE3f/M3AQB//ud/Xr59MdrRfiuAOvmjPADs27cPt912G4aHh/Hmm2/inXfeqZkKcztVs4/QJR/jqW30TWcnCpJues553U44v+KBjztdOGKgesylVQGUjjnpeWe8cU4ZhbIjNxVRoM3lDeMaHh5GT08Pbr311poLy7V8lOWIubCweM3W+Pg4Tp48iZ6eHnzsYx8rn5rk2Jxu3ii4xXHnyGzYbBgsq7JnllMvc1LwVubEwKZsKxiSGqMDDWeoWn98q3yRcjuqHGDU0Zlx6rKNMngGbrUTBWcFeNaR/g45OKWUsH79evzGb/wGJicn8ad/+qd47bXXMDk5WaMPnqbqhknUq2AVvx988EHceOONuHTpEl5++WWcOXMGU1NTy9Yq3ePaom7ut44720Sun6F3nqGwHuM/A7KmICPBIl3gZ+DUAP+RYZQ8jWTDdwxCI2gkHsCqXfFcHTlWODU1hYmJCfT19WHTpk01ctZjXyHL1NQUjh07hpmZmXJnO86pEbHBOD05lqntszMrG9CA4zY0VB8OGLivvKarTqpgonpShqUyqyxRnvvjfuu4Oz0pWOemexwAXCBm2cM2WE+6YcH9VqCOsWhtbUVLSwu2b9+Oj33sYxgeHsZzzz2H48ePL7MNZmROf6wz1kl7ezsee+wxbNy4ESdPnsSPf/xjDA0NLVun5CCha4TqSwpEPHbKnKPPamdVywiaqoKsBm8+x+OmrFTTqgJKdRRnkBohXGIFaXk+r2X4m/OPjIzg8uXL6OvrQ2traxZY3bQxpcXX1B49erS88+b2229HS0uLfalX6KJK1sjjNn0c+3RR2gUP1YNr29WhvxWYFVD5uAssCqLqWAwwKr+OB7MoBRJ2ah1TZUjaroKCAoLK7exVdaDslcs98MADeOyxx3Dq1Ck89dRTuHDhggWssAtml8zIFLRTSujv78ev/MqvoL29HSdOnMChQ4cwMTFRUzcHbx4Lnmq7sWR2y3UpOdDNOfcGThfcdKyjjmg3lkH0ioNYMlkpq1w1QBlJnd9NlRiU3LQ5EivX7YyzwpXFxidAbmxsDJs3by4vC6mXgoXMzs7i6NGjeOWVVzAyMoI9e/bUvL+G225oaCiv2Qv5c5HSrauErhyYsRO6aaKez63bsL4c41SQccHDjRWPZ9UMgtsOQFCg0TLKDmNsnK5cQGGQ4N8KkhoUHNjGf9WpBnb939raikcffRR33XUXjh49ih/84Ac4c+ZMWV+0kXv/C9cVcvFMYvfu3filX/ql8l1PR44cKdcr+TIeDhYMeBo4eOOE+8PjEzpUu2O9Kji6ax81YPGsSmcebllCd9BdWjW3MLrOR3IRWZ0ijN+tZ7nIrgDA5xiQh4aGcOrUKWzevBm9vb01DEXldP25cuUKfv7zn2N8fBz33XdfzdNbIul0iafjCuTKrpyuWD8OwLgelsEFGtab04+L5twmj4kbw+ib1uXycXKBlPvCwZHlYvDOgRPLrAGD19KUwUR/XbBwLApYvJhdmTU/pCLajKn45z//eYyNjeG5557D9PQ0vvCFL6Crq6uUPR7gwbpUGfg/9+nuu+/G4OAgXnnlFbz44ovlUhMHF9a/Anwu0OpUPPrJGz5MXHSZQi9IV0CNtnW81CaZWXN/cj7MadUwSsdcNOprZHYOxHmY6nN+B7LKAubn5zE9PV2+M/m6664rjVrBg+WLSBr5hoeHMTIyghtvvBE33XSTBUl26oh6K9WPW3fJgYgCEDsXO004KjuWsgedysUxBidmkzwOrGe3UcTfOfDW8eP+aB4u54BY1yWBxUtLZmZmLNth2RgodIrLx/W6RB0fF8TiPNtTe3s7PvOZz2Dnzp148cUX8fLLL2NmZqZmVzg2NEJWBqHoFwNPjFlTUxM++clPYteuXTh37hxeffVVXLlyxU6rGdTc+qgCprJJPcb6iTwMkMFqY5dbN5s4ACi710DEOl3pheerBijZuTS56KDHdJAYNHiKVsWGFAgmJiYwNjaG/v5+9Pf3L2tf5dNp18jISPlA3wMHDqC7u3vZYOTYEtfjdr61XaD2+Z050HDHGegjj+qIjzt9azldX1bZ3XFlCvpbA5QDPHZ+rd85bEq1t3RGm1WzGAVPBTvWqa7z8a606lTHk89zANu0aRM+97nPYePGjXjjjTdw+vTpGvkVdLldrVv109XVhU984hMAgIMHD+LgwYMYGxurCSRaXxU7i0DMcgTgOXvWXfu4fMiNYYwHLxGwbHqNKQOuA+6qtKqAkqcbzhncFC8HmEXx/m4uA4FOf1RB3Nbk5CQmJyfR29tbvvJVI7w6R/yfmprCwYMHMTQ0hBtvvLGcwkRiw3fsLRfhlF1zNNSpjQM6dWrOwzvXrCfXLh9zO/GsTwckfN6xlaq6HENQsOV+c2KHjbx8qVXIye2wTeUu12I52CldYK0KPC646Vp8Y2Mj9u/fj0cffRSXLl3Cq6++iunpaTQ2NqK5ubncJNTrRIM1xh0zfI51vXPnTnzsYx/D7Ows3nrrLRw7dqzUw9TUFKanp2uufQz2qrMpDjqRj0FNdeVYohs/HudIuess2Ue4jJO1CixXBJQppRMppTdTSq+llF5eOrY+pfR0Sundpe8+yv+1lNLRlNKRlNKnV9JGCMkX5rqOx7c6LBsSMxpWvObR/9rW0NAQxsfH0dPTkwVtVi47YOweNjY2Yvv27eWrSXVtjDeGuLzKwgaixsJ1uADiznGfdaqurKRKR8qMlA0z4MSx3D37ji0qgFb1Scsr6GgelkFlccxLQaWqbQf8Mf56bWmcU7vVdtSZ77nnHtx+++04deoULl68uGycooyyNxew+HhjYyPuuusurF+/HpOTkzhx4kTNtZUcQPipQ7mpriMV3CcFMWXdXAf/D3YYsuurNFgHjqXXwwBOH4RRfrIoijuKorhn6f9XATxTFMUeAM8s/UdKaR+ALwLYD+AzAP4wpVR3m9hNR9z6VhU70t+OmfA3G20cC4AeHx/H1atX0d3djU2bNtWwCjbiMBY2homJCRw+fBhjY2Po7e0tH5IabXEf9LKI6Lc6jL7wnfXjAKSK1bFT6oX8Vbp0EVcdWB2J21T91Qt0eikR611l0B1zduQqthBjwPncelZMF3WslWnzt+owEi+VcGDJyaeOHsdbWlrKJ5a/8soruHLlSk0erp/ZmoK5yrywsIC+vj7s2bMHDQ0NGB4exsTERFYuvbWR2aUDO2XlubrYjqJeBV1d32Z/4rFRnbA8Lihr+odMvR8H8M2l398E8Ot0/FtFUUwXRXEcwFEA99arTC9oVSWpQTngiuRYhVNO5I23BPLxiYkJLCwsYNu2bejo6Kg5r0yJI+Ps7Gz5srHGxkbccMMNaG9vr5GNDTKS65uL/ApCurTA5diI2Kl1KqcGo0zLARuX02+WN/dYslzg07ZcnQ74nB0oeHO7gN9AcOyDx5nlyl2DqXJpfj7nNkoiULCt8FiwvP39/fj4xz+OsbExvPLKKxgbG1umGwUH7rebPgOLb4vcsmVLGRymp6eXERd9p3auj6x79mMGcAZWPs/lFWAZWLleHi+VWcdRbSKXVnp5UAHgqZRSAeD/XSy+k3tzURSDS40NppQ2LeXdBuB5Kntm6VhNSil9BcBXAGBgYKAmympyEcHdFePylx0gpWkZdtpQ4uTkJFJavBi3tbW1pn0tGwYCLF6cfujQIczMzGDfvn3Ytm1bjbw5Vubk4rx6B4+yLT0WSdlcyFw1zXBg4/IoS19YWFhmmAy8Gmy072y8UZ77obKrsedmEzE2zFL5EWCqKw5GLHvUw0AFvP8uJAfUep2gMngNGPGtu+s5/S8sLGDXrl0YGxvDyy+/jPXr1+PWW2+tyaekgm2W7Uqnvtu3b0dHRwempqYwPDyM/v5+20eVkx9EwecbGxvLjRkFrarLrlwgZ926pYyQJ3QeY6k+xI9sq0orZZQPFkVxF4DPAvjdlNJDFXmdBy7ztqIovlEUxT1FUdyzfv36ZZeJ1AhpGEnuYbnKUNyakWNgDKLz8/MYHh7G2NjYsovDuW52wHCmEydO4NSpU+jq6ip3uqNeB1ravu6+shHHBlD8V4DMMS1Obhqseqsq546rXhQg9Tjnjf+6zuUukXLtan0si7JDtQWWQTfE4rwyvhgX3QVmRqUbFsD7IM0sKB57Vu9OLDed1DHau3cvNmzYgNdee61cr1S5OaCHzeVAdG5uDi0tLVi3bh1mZmbKdcqQQa+NZYbHeuKdaD6nY8RgyP85XyQdVzfrCL2xHTnCsNK0IqAsiuLc0vcFAN/B4lR6KKU0sCTAAIALS9nPANhBxbcDOFevDZ0WcoTlgYk8ubtyWLlcNgZUo5GbUs3OzmJkZAQdHR3o6OiocfzIo0wAWLx+7dKlSyiKAjfeeCO2bNmSXX/j/rq1WHUwnsLGbr7WyVGZk8qvoO0MSAONTlHjW43OsUV3baWyJB1TJ5tj5jzDAPxGQCRmd8rs2CY4hZ7ctabRb7e7zP1iXTBr5XVDzpubIrKe2B4aGhbv6Nq/fz+mpqbwxhtvYGJiYlkg5MDA4MNjzfpuamrCrl27MDU1hYsXL+Lq1as1Y63vtGIw1h1/B4SsvyqbU791ZXgMeTbA8mo5B6K5VBcoU0qdKaXu+A3glwEcBPBdAF9ayvYlAH+59Pu7AL6YUmpNKe0GsAfAi/XaCWF5YNWYuUM59qQMSc9zVM2VnZycxMLCAtavX1+zW80Aw/+jXKznrF+/Hrt370ZLSwuA5WuC0V70Q6cgrA/+razDPXxWF7fVAZRhqc5dW5FySxoKjgwCyi6B2vuClTlrv3N913b0v7t6QIGI23ZBgMElQMEFbD6uzs6/oy889rxBpDaeC+iufykl7Ny5E7fddlv5sjqdUvO6JwORtsGBYPv27eju7sbY2BiGh4eX7XzrmFWth6qNc11ubB0Qat81mLCPcluOlUZybFTTStYoNwP4zlLnmwD8WVEUP0gpvQTgiZTSlwGcAvBbS42+lVJ6AsAhAHMAfrcoiup5FPzOtYv+/M0d1WmeSzrF4ba47MzMDKanpzE9Pb1MNpUv/qeUcO3aNVy6dAmbN28u38xYL1JxYkfRKK0P0KgHJu47gNSxJj6W03EYG7Mix/xUt7kxYRDldSQtC9Q+NEVlY6dtbGxcxtSiHgYMroeno05elZEBV2cE2he1Tdav6kSBTYEk1hRZnsgf/2+99VZcuHABP/vZz9DQ0IAbb7xxWXvRZ2VVbC/xe+PGjVi3bh1OnDiBN954A1u2bClvnODptAYKHisF5uiLkhYut9Jj/DoOBUxuh9uvtw/iUl2gLIriGIDbzfFhAI9mynwdwNfr1c1JnTM39dKIoQ7OBhmDGRHSsSRlrHE+XvbFeYDazYQw2qIoMD09jXfeeQdjY2PYtWtXzQYQy+9SQ0NDzT267jzw/lpVGAjrR/uVO6dTQe03A2n0NcckXf08JmqUCm6cuC23Js3yKwPUcde2VEbdyIk61Wk1iCo4cVtOnzGuDNCRh8GZnTrqcRdgB8DptDzkDbu499578b3vfQ8vvvgiNm3ahJ6enhr2yPk5hf7jE3ravXs3jh8/XrLK9vb20rdCLjf1Zdly78/h1zRUzTpySddvmXUqM43EAWmlRGbV3JnDnVGjj2MMUJE4gikN198aVXOgMTExgdnZ2fL9New0vIjNMk1PT2N0dBRNTU3YsmXLMhk5qZEzI8ixYZXRgYLThX7c9apORgc+TocrYZLh9LqWF/1xjMABAbMTdiCdpuklSaqDHHNkFs9OH3kdCABYBl4BSBGk2fn1KTxRNwOZXibj2o97tvndO9GHdevWYe/evRgcHCyn4CxXlFE95Ka6O3fuRGdnJ6ampnD58uXy3vJoj8dY7ZSn89wHp3u2T5VDx1FtiIFXy4RvcdDnsk7fmlYNUNabFiqL4PN6DKjdFNJpgeZXY45rxmI3XtmrW/cI42ttbcX09HTpJE6Ov49uNKlcTidVOlJ53HQ6ynD9OvWtF5jYKfix/1UgzeXclCtknZ9ffA86s1FeFwt5c87HtqV1c15+94q7PMz13+l+fn5+mW6VILh+cl3sB9rXsLf5+Xns3r0bTU1NePvtt8uNHWcvGhyc/nt7e7F161aMjo7i5MmTNQ+5raort6mmY8B9cXpjfUS9kUd1yGV0L4LrYHuvB5LAKgPK+GYjZcNway1uUwXwF6RXMa+oLwasvb0dra2t1pA5b8g5NzeH0dFRnD17Fu+9917NBbpsDCvtPyf9HxfIu2UJzs8RnOtW5gLUMmXOzx9dqNeozX3h/Kx/rUPr07FwU2puL/Lpe6HZUdyyCwOoysp91qlbbkqf0xmXVSKg46SOzUCYc272Ec63bt06bNy4EefOncPQ0FB5nuVXhqtBI/6ntHg98cLCAoaGhnDhwoWyDrUp1z/WvWN8rD8XPEI/elG6BlJuV0Ffx0dnBlW+CawioGRluzUFnT64nWtmSew8rDS9FMcplS8gdsDF7UU9MzMz5d08Y2NjNpq533ypDxumayfyaz0OWKtAwgUKBozoE8viwNJFYtVn/NYNlpw+ou2V5OXXoeYAz8mrwKIBk8cibIZnHKpntVe2C50WRzl12NyaGoNogIySg9w7wwcGBsqLxXNPDHe7/gqaALBlyxa0t7djfn6+fFUEl+NZlu7IK+FR8IpvXZfMBQWdZcZHZeAH7LANcrCMVG+2t2oe3BvJrS256wx1lxRYfnmLPnwh6syBXzhWURQlo4zpEhsm73Ry3XGrV2tra1nW9YfLORCIQdc1VNWJAirX7Qaedcl5WR9spMretB09xvXoxgVf78j9rBofZrdOh6EnDarRL8eMlDXl+qX5WEbeUGJ7zOlFkzo6y67jULXBo3VFCsCJN3wODg6Wbw7VdljnCtq8UdLT04OGhsUX5MUzKvmhwzpO2ga3lbvtMM6zDthPNOg4ZsjnglRwcHP5tW6XVg1QKjtiR+VoEUaqZXSxluvg35yHN2p4UBoaGtDc3Axg+UW1DuB0ytbb25u9bEHlAmoZM7eRa0+nEXpeB56fyKTTWQWpcBD3hB/tLycGSGVpGiwiD1/yw/W649Hvubm5ZRfbq55WApwMgDz2DogC4JkB8XjHMQ04nBT843esKeqU0rEex/pUPyxz2OHIyEh5qx7Lraw0F7QDFLu6ujA+Po6pqama3WrNy8GJL8Vy/lA1/VVQy12dELp3O+C6y6/rxKy/qrRqgBKo3aFixfL1a3p5Ru4dMwo2zikYnJglBEiy00Y9ChwpJczOzqKxsRG9vb3lE4f0vTpqzHFMp3/qfMzEcpE7F3VZr9x/zqeAztMgDkAOiPiYjpm2EzpWMMk5px5j51CmrkFH5XFyKVvhbxdIXHDmj05BNa+b5uc2naIenj5ynQzaWj/3aXp6GlNTUzX3V3N+Hmu91EeDVHNzM1paWpBSQktLC+bm5pbJ4Zg6kwDWqQO7+G5qalq23hllHOtWvejsIYcFUc4t42laVUCpTIsdgtfx2HmZJapxK6N07XG+AKS4/s2t6wC1dzjw7nZDQwO6u7uxYcOGZUxN2+RzznkVGN0Lzbhe1pXWE/0C/LQr5NL1PMe8c8k5KjMK7a+WjW9dU3WAy2W0fQYgZaaqM/7vHDy+c4Ei+qtTO51S6jRZ18p4F5cBQwGcZWY2qH0MgL127RoWFhbQ0tKC2dlZew2o9kkBKmyc9dDa2loCpQY+nW1wf5UNsq7ivD6flevmvQkeX2bF/NrlHBvnFLOn3A59pFULlCmlGnCM7zAmdWpXB5fJ5WOj4+kXTxnUGNiJoq5I7Owc0bQP8ds5gTpKOKMuP8QAc9uu3w7g3Bpr/A+ZmcFrf3Wswpl4g0Vl0mkP3xbI7MYl1bWOMQcVp2OnZ2XjfM6NP/fL1e+AiAHTtRP95ymiToNVRu0r2wlPL4eGhsq3h0Ybjmm72Y62Hb9nZ2cxPj5ug6f6iG7OxG/d7OHEeuD+RXm1Jx4fneazXWlAi8T37FcRgVUHlPybp73MUJj9AcsX9dloI+WihU7bg/qHwbrEThaD1NzcjM7OTjQ1NWF8fNyCtQNxNQi3gROJwcddaqLTisgbfQq98vTKXRbkdBNJAYKNk41O63GGqJfnRFl2fG5XQcoFKtUng4vWywGMN5zY2Tjg5YKagg/rXvPoWiTrIXchuga10JebIrNNnjlzBk1NTeUTrKINTWz7OiVleTs7O1EUi48gnJqaqnkHFPc1Nz3mtqLfLLvm5eSCoJ7jZZcAdWW4zgbdcU2rCihVCcpY9NFqcU5ZBjMAd16jMKdoQ6eyQC1wsFPFdKS7uxvT09O4ePGiXUdz7IXlccbudKPMpYo9Rt0xLeF6ee1XHYiB2AUddgRm4jH1cdMevcwmF7wYyAMUckEm6uDLY/RRdI49RTvMcKJs6ENtSJcXdFrtZGNHdE7J01IXfBjAmFmprrS9ubk5XL58GY2NjeXaovY3N57xHddZxrHm5mZ0dHSgr6+vZtlJ5QvdcDvKKkPXRVHUvM6B7TPGPm4qUJvLBVX1bw0Q2lfOl0urBijZABXgGJjcgHIZdSCOxq5+NX7ekeWkjsrHwhgbGhpq3jCnDI/bdXLm8rGxKQix/jQpsGoengq5KK2ArKw4x064TQeY7OzuagXWjQNlpxf37QJjbrrtGCszn6pAF8eU2TjnVIePfAyK3Fdl6coiNZhGunbtGq5du4b29na0tbWVfVDW5gCNwTT0EkGws7MTe/fuRUtLyzKZnCwK9E5vfJmRjhkf0985gNeA6Np1l+5VscpVA5RALfMA/JRPF3A1X9X1bDpYuTW9MKrcJTIsQ2NjY7k219PTg8bGxvIOAjYmjoKcHNvNye4G0kXU3HHVb5zPRVIHlrpuzElBvKGhoWbNkh1GQdfJzUGMg4vb/MgFSmUWqg8GBM4bbEbPuaCq8vN5ZVhONp5G65JO9NddFqfAzLq4cOECpqencf3116Ozs3NZWe4n64NBXtlrUSxeX9zc3GwDIbNmTbHBpHYZZRj8uT/Rbg6Ec0lnBCFr6FX1WVUXsMqA0iXHLiNp9Mkx0ly9/O02LdzDFVz9semxbt069Pb2YmRkBKOjo+W7drS8piqHCzB2juryB7ipLrgt7iPgQYT7peBeDyS5jQgkCiAMRK5eB2qOKcZ/Z/TqdPqGTwZcZXRVIM72xmVzOmJwV5DjIOICt858OB8zXl6fm56exs9//nPMzMzguuuuqwluDsiUnfHaZJwP/cXdOVGXMjztEwMut8/rulX2lQN4DYzaH9VPtMPLTZE3bLEqrRqgVGbFx100YSPlsqo4ZgUarTkPD/js7Gz5kAuNQuqI0X5DQwO6urqwefNmnD59GidPnsT69evLjaGQXeXOyVulF2dI7nxVoOClDJ7S5MAr5HWOxrvxkcflDx1y4rsndGyjf7kgyY6srE915BxTd+gVMHW5h/vhAI/bj2UYBiEOGCorg50DYDcOru2FhQW8/PLLOHz4MG666Sbs3LmzxoaUqStIuyl0/J+amsLU1FSN/aqtuGk4jxmPtfPtqM9d0O7qit96sTnXF31gkNSAWI9Vrpp7vWOaxs6rzE0BQPM4NqDRSvOwoiICxsW1ce+2Dqrb5Ekpobu7G9dddx3m5+dx6NAhnDtX+waMGIyqAXHyan41fM2ba0OneFwX94vrU0agKcdyVJ4IJurUKiuPvZtSKrPSpQ1mCgoAOb1GGb7I3oGoyuw2K2KNuoop6wZOJLV17SePRfxWVjw4OIg33ngD27dvx/33318+FzVmB2rrIYt7N3zk4Xd3KwNkvfCT2jmf6l6JiQsGPE5u80mv7dQrDwIQ9YqB3I0gVT4JrBAoU0q9KaVvp5QOp5TeTik9kFJan1J6OqX07tJ3H+X/WkrpaErpSErp0ytpI4QNQ1UwCmWwMWk0cGzUAa1jHJwCsOOSGnbC+GZZUkpobm5GU1MT+vv70d/fj7Nnz+KNN97AlStXslM4o+caI2L53JqRBgp2BDYM1ifLwNMgN8VWoMslBg69y0PzqA6ibzw24XAppWU79QxG7BjKCtySiQJLnIt+s+7Z+Xh9jYGG6+F2uL+uX1GOGZ3TszI0Bkmt+9q1a3j99dfx1FNPYWJiAg888AD6+vpq2mFwdQGIWa0GA7YznjUooOu1k/w6XLYvvXhcA5Dq0/1nPXP7+jAcLevsUkmHppVOvf9HAD8oiuI3U0otADoA/DcAnimK4vdTSl8F8FUAv5dS2gfgiwD2A9gK4Icppb1FnddBKFtko1FG56JwDJ67lITr5frclCvKx1qMRlp3/WIca25uLl8cf+7cORw+fBhdXV2444470NPTs0x2HXiOkPFf+6+szJVzQKt916kLsPwCcKcXl3jscsYc5xXow0ECrOfnF58HGreFNjc31zg3B69wKJ7e8+VdDIAsS24qzv3NMcLIx5fwcF8iMWDoTQHcjral4MU6CgCIOkIvo6OjeOmll3DkyBGMj4/jtttuKy8yZz+KunScGBCZBeolXLxTrExS7dHpj/WhbFoDi5uhOPti/SmgcwDgxCRMcSeX6gJlSqkHwEMA/rOlymcAzKSUHgfw8FK2bwJ4FsDvAXgcwLeKopgGcDyldBSLb238WZ12rLCqZM3rQERZVaSqdSpWbDz5J9aY3KVCDE48OI2Njdi2bRt2796Nt99+G6+++irm5+dxzz33oLOzcxkL5jpdv5UxqyFo8FDdVDl86IONUh2L2+WxYEbBDhyAGk9SmpqaKkEwjk1OTpZ3eASwNTc3oygW31d0+fJlnD9/Hl1dXeWO7dTUFACgs7MT3d3dmJqawuTkZBnMYjziAuu2tjY0NTWhra0Nvb29aG5uLndri6IoL0Zm1si6iieINzQ0oKWlBU1NTeWzG2OcmZWxXTiwY/264zxePK0MfbOO5+fnMTk5iUuXLuHy5cs4evQoBgcHMTs7i/7+fuzfv9+CpMqgY8rH3CxsZmamHCeV3dl0FesOu8mtAyuRYX0rMHNSUqXtRr0cwOpt5AArY5TXA7gI4H9KKd0O4BUA/wbA5qIoBpeEGEwpbVrKvw3A81T+zNKxmpRS+gqArwDA1q1bazrHLDHeJaODV0Wf3TSMI7EbOGZVMV0dGxvD1NRU+e4c4P2NAccsY0B7enpw4MABjI+P45133sErr7yCoigwMDCA3t5erFu3Dq2trZVPJlJg5GPcD2Vp2i9+tiYDWzzOf3x8vAS0YHHMyuMtlOzA4+Pj5VNkGAT5vc9FUZRPeuc2dZE++hxlZ2dnMTExgatXr2JiYqI8zwwTWLxOMN49HU/JmZubq3lIytzcHNra2tDW1oa+vj5s374dfX19aGxsLC9ybm1tLe1samqqBNf5+fnyYRItLS1obGzE8PAwAJR1tre311yREEy4tbUVMzMzGBsbK/UWsjc3N9ew46gnNn9iXGJdcGJionxaT+hwYmICZ8+exfnz5zEzM4OZmRlMTU1h+/bt+NjHPlYTkJlRrQQk+T+DYeijtbW1fHC0MtDcemN8NAC45R72VZWD8yhQ6+yBz8dMJY7z7FMJTy6tBCibANwF4D8viuKFlNL/iMVpdi45DrtMgqIovgHgGwBw6623FnS8RmiOjHppQU2jZmqtx6su9eEyLS0tJYiMjIyUt25FfmZOOoWNi883bdqEu+66C0VR4PTp03juuefQ2tqKgYEBHDhwAAMDA+jp6UFHR0fNdEYHjJ2YBz+cJu6cCHYUxxsbG9Hd3V0yodixnJycxPj4eMnGGJiCQTU1NZUsio2yqamp3MUPPURgCZAJXS0sLGDTpk01wMA6CpljbTfaDXZ44cIFXL16tbxdLoAm2j937hzOnj2LqakpfOxjH8ODDz5YEwimp6cxPDyMkZGREizilQjx5PqGhoYSaEZHR8vXFMcYzMzMAFh0rsuXL2N2drZG3s7OTjQ3N5czj+npaSwsLKC9vR2jo6M1T+4pigJtbW3o6Oio0UVnZyc2bNiApqamcuxCf+fPn8fVq1fLnWYGvHiCfly7OzAwgL1795bAH3JHfgc4kSeO8wvAGFji3NTUFHp7e8txd3UwICp75HrZH3lJgO8EUrDmWZ/2xYGvyqVlHOPNpZUA5RkAZ4qieGHp/7exCJRDKaWBJTY5AOAC5d9B5bcDqN3+NSkUpQvqHJGqKLKCIzNMp+zIx9+RGhsbS8MdHh7Gxo0ba14Xy+tq4bjM9AIwdu3ahXXr1uH48eM4deoUDh06hMuXL2NwcBCbN29Gf38/du/ejdbWVvT09KC7u3sZKw7QGB0dLf+PjY3VAB1PSRYWFu/rbm5uLkEHQMkMQw9xb3oA2bp160qwa29vR09PT3m3kQJa6CyORRsBruEUwTx42uSeSs7jEYa+devW0rF0BgEsrpdNTU3h4sWLSCmhr6+vJoA2NDRgx44dNRsXzFB4+hWBIxgts+CFhcXH7p04cQILCwtlH0NWACWjC/Cdn58vA2AwMGZTIyMjJUOcn5/HpUuXAABXr17F5cuXy3ENFhfyBjg3NTVh+/bt6Orqwrp169DR0YGurq6al+EpSKqv6dTebRZyGh0dxcjICHbv3l3qP+rXet3F3uyHUVbX/1muuPIg8ussKtrjstFOjL/6fSTGEdVFLq3kdbXnU0qnU0o3FUVxBIuvqD209PkSgN9f+v7LpSLfBfBnKaU/wOJmzh4AL9ZrR9dE4reCWm7tRTtaj07reQbXhoYGbNiwAVNTU7h27RpGRkawYcOGGnl0SqwyNDY2oq2tDZs3b8a6deuwc+dObNy4EYcPH8aVK1dw8OBBNDc3480330RbW1tp6OGY7e3tAGqfuDM+Pl6yPu5DS0sL2tvbsWXLlvIi99gZjqetd3d3o7m5uZzy8xKDMg6+VEvZezgCT8s5kMUUsiiKZY9M07FiMOLEhhyBSB0lGFpvb++ysQ+Z3SUvOsUEUKNrzhvH5ubm0N/fXwbQyBPfsYMbutFr+gI81bYZoObm5mqm2byMEfl5SSEAWwMYj2GMkwIgMzv2AcfAIt/g4CDm5+fLGYOOK3/inLvmlmeFunPOY6fj7WZauqEW9UfdsQaqtuFSDicirXTX+z8H8Kdpccf7GIB/hcVLi55IKX0ZwCkAv7XU4FsppSewCKRzAH63qLPjHYmV4RSoSlOWyWUVUDW/ghsbdVNTE/r6+jA/P4+rV69icHAQ3d3daGlpqalf21T5Ajg6OjrQ2tqK3t5e3HLLLRgaGsLx48cBAF1dXZicnMTIyAgmJibQ3NyMrVu3AsCy15EGsw3WGnppbW3F5s2by51OjuJ8G6YCShhklNGIzmU0ILChMrsKJslj4KK2LndEX3KMhsu5b5aF+8cAovbCMgQTj36GQwcQtba2LrNJnc4xe1YQ4Kkz34TAwa6lpaVm6s1gFuyT+6HkguvUNbz4zZtWborKM5nIOzMzg5MnT5Zr6zoz0DFTvfAYxKMBdUahfujsgMc89MmzBO4j258GeidXvbQioCyK4jUA95hTj2byfx3A11dStya3/ih1W5qshsfHIjmgUMDkqXNMe0+fPo2Ojg5cf/31NdckKqiH/NqfqLe7uxsdHR3YuHEjbrjhhhLEpqeny/Wwjo4OtLW1LauDDdglZogKBqofNaA4xkDFeuRlBadzTWzsIbNj4prUuRVE1amijehvyK99cIxVp56sLw6ael7zhbzqrHoJigaeaD9YdbB4fk8TT4dDZn30n5NdZYk2eZmGQZzr4GWJSKdOncKlS5dw8803l5tp3CYv8VTJqEDNtsABWxk5y6wXjDvW6myLA7vTXb20qm5h1Kksg6ayNc2v7EEZh5aJYypDtMtR/sKFCzh79iw2btyInp6eZe2xTC46aj/CIULWtra2ZYwvx+hUdnZCXtdxfVMnVSboZOW+8jeXccDGjCGSKxd1qsOzHnS9TcGSGWCUUUBgnSrL0nbVqTS/BozIG2wpZ6P8JCDVHU+ZY+oYa5s8xQw2qmDCwYEZPINh1KuBi9cONbgVxeJmZGdnJ3bs2FHDurk81+kASQGZ5dKApuDIdevaqLMdtmV3re3fJ62aWxirABF434AiaT6uI37zIMYxjZZaX5yL3ev+/n50dHRgaGgIb731Vrnor6AQdbjNKAax2GTRdcLoY6yD6SPtA2C1DNetTLFKL7n/OWBjvWkQ4noUmFgPkXRKzqxQkwNJlccxJAUiBUmVJ7deq/3W4KyAr0+c0uki65mZvCYXJNW+uF8cMEJ+vlNK1xAV7DkvM8KhoSFcunQJ/f39ZUDnNUC+bdEtfSib5MRjoOyOQdIBLOvJYYNrz0279XcurSpGyb8ZjBxLckxBWYxjAhpl2fiU8cQtiRMTEzhy5AguXLiAkydP4rrrrkN7e7tlUpFU+Wycem2jW38FalmGA1/eiee6OHG/XfTmxNMf1Z/WqWPA+o8+Rp256SG3y2Pn1p1DJl5/UlmiLpZf+8RysDxcB08BdbbAZaM9vUaPGaD20zl6jCXbrrI+F8xy9fP0lq8a0MDndq6Zuc3Pz+PIkSMoigLbt29fdksvt+HGiglE6JV1oGOoMqitMEPkssr83a47l2HZ3LuoXFo1QBkpx26YjqvRugjP5Rks1BA1+rOzNjQs3pUxMDCA2dlZvPPOO3jjjTcwPz+PXbt21VxwrExOp1+ujwwKburhWEg9PeWmsGo0cUzBQXUaSY1agdaBu7ahcihLdH3T/8z61AY0vwKDyukCnLKzOKbn2eH5wbN6h0/IzPnV9qqCiksBMm5aqUCp5IFlDB0wI2T5zp49izNnztQsOcV6pLvXXvvAO9vchhtLHpvQiV5/C6BmMyjqYGzgYOH0pnpd6ZR81QClAoMaPitUF4MdADiQ1fPMVHWKFP/jfcbbt2/H9PQ0zpw5g6NHj+Lq1avYu3cv+vv7s46kLM4BV5WcCvR8XAE6B4p8XoNKOJzqXOtygYdl0fEIZ+T+xVqbBiY3fpFXQVHXUxVgWUadMjtWxkGCy4SjuXbVqZ3enB4j6bJA9FEv8I6UuxKAN650jPgZoC4QcV8YWDjf1NQUDh48iJQSrrvuurKMA/Ng33oJVPTHya+6Cebv+h93UfESFV82F+1zfzQpOWEW/JGaekcn3bTLOYOLQkzxHeiwU/Px3KYMgPJynO7ubtx0001oaWnB6dOncfbs2dII+vr6ykuHdD2U21ZgqzegTgfK7BQEnF5ygKIgy4GoHrNREMm1kwtwenWDY//19OXYtuovQFsDQi4QM6tipuJ0q0ws2g2g4nZ1nU0DKdfL+XTzJRIfY7BlsOIyOabFd+OwXg4fPoxr166hr68P69evL9deOb9jpeyHcUmTbmItLCy/8SBYebQft5nqDMUFBqfbKKOzNB3LnB9qWlVAqd/sMJFy06X4zrEGd9w5F7fJU6MAwt27d6OlpQXHjh3DuXPnMDk5iT179mDLli01t49xH9gguF7Xf8eq1eC1rBqL1qHAE6lqeuz0yfIwsKgs3LabiuamvVxH5GNA5st6csatwS6+mYFoUkblNg8U8HQqzWtn7JzuN9fLZeO4Anp8dErNunBgz/oGUMNalUEy6J49exanTp0qH/ASV2iwjNEG76TredU564jrUB8Jhhh1cTDmcQ6dshw58sHBmdtcCUgCqwgogeWXhSijcNPJ+O+AJ9iiMhnOw8e5HmU8sV7Z2NiI6667Dh0dHTh+/DiGh4fL+6XjtjKtC1ju9C5ppFRA57pYX66c06nWmwNH1aEzZjcOmhzbilS1A5oDZ66Dv916M1C7+B/yFEVh15WjDbYVd94FB+0z913zRnkGL9Wn5s+V4YCjgMy3/aksUWewQ2aA586dw6uvvorx8XHs2LEDmzZtKnXpwJhZsQJv/I6ZYjBtZZ+sUw2SzA5ZNxq0os0gNmpf0Y7bHV9JWjVAqQ6g5/Q4P7VF8yojUnDQunLrdDr4Mf3o6OjAli1b0NbWhvfeew/nzp3Da6+9hnPnzmHnzp3YsmULurq6KvvC9bqNFDY+B145/VXpQfXsgpAaphozOwLXG9MpBq6oh6d9vCzyQQyVZXW7ysww2Imc47pv9+guLqO3N+oMIfJwHU4GvlyH6+N+OJCNPMoENZ/qhc9r3THG09PTuHr1Kt59910cP34cExMT2LRpE6677rrykiBl0Ny/qJv7yu2pnbFNReJ1x3qzhty4cOBwG2ixa89A6nzQpVUDlJHUgTiS8DHOD/g1Lj1er10GB8eCwknj09jYiI6ODqxfvx5nzpzBuXPncOHCBQwMDGD79u3o7+9HZ2envfnfXcbD7agT6m/Om8vjBl7751iMywcs6p3vzeYNBWXRwVZCV3oJDbD8Uh/VBYOiA1cNIPqJPOqsCnbMfBlM+F3oIbsCHLM8Bls+H5sR/CxHB9xRD/eXgc2BJ0+jgcXbXuPxdlFfPLCjra2tXPuLx7VduHAB586dw9DQEEZGRtDY2Ij169fj5ptvRldXV1ZvPP6qc+0/P9wikjJJPla1G8795uUYtWGuV6+z1JsT+FGOubRqgDLn8A4A9Tc7XDiXDq4DDR4QB5JugOJ/3F3T2NiI3bt3Y+PGjRgcHMTRo0fx3nvv4ezZsxgYGMDOnTvR19eHjo6O8uEFPNCuvVybOi1hHSnrVn1yUnai+TkP16dBiA2Uf8eTb7QfatC8oK9t828nl7sGMGSMC6HjsWcR1IDlbyyMMZmfny8fRxf33c/NzaG7u7t8XmTUHbecBkMuisVXMYyNjdXs7MeNA/Got6grACtAjuXj/kVQCvCL510yOE9OTpbyRB5+7mc8KHlqaqpsf2ZmppR3fHy81ElXVxc2bdqE3bt3o6urq+Yi/BhjZa88di6A6XokL3vEeQeWPEahp2iT22fb4roUQFUmtlkOhrm0aoASqGWMbk2iynn1GJB/spCuj3BeXcOomuY2NDSUd9q0tLSgs7MTfX19ePfddzE4OIhTp05haGgIXV1d6O3txYYNG9DZ2YmWlpbyaT78oFwHChowXN+dDqr0G781mHCfFdgU8DgqB3hMTEyUDhhMgg2+KIryoQ/xbMZ4Cg7w/tIGMxT3dJy5uTlMTk6WgBDGHs/eDMY0NjYGAOVmRABNgE1Ki6/86O3txczMDAYHB3Ht2jVcvny51A0/MzSlhJmZmfI5pfFg4HiYBT/EJIJpXEgewNXS0lLeyx/9Y0cNB4+nCAGLl+rE80Ln5uYwOztb2l7oI8YjgoQ+7BcAhoeHyzGI+8v7+vrQ29tbPmVq3bp16OrqWvZkKO4Dg5AyQxc41WYjEORYoQvEXJeyapVJWS8n/b/SGeeqAUplK6wgnj5VMa0qwHDTg9xV/gw+ejmOMtTIk9Li2mVzczN6enpw6dIlnD9/HkePHsW1a9cwODhYAmpHR0f5pPN4oXzUH0+qCSdjIFGZdG2VQailpaXGiHVqwcYa00NeZwz9TE1N1awdzc7OlkwNQAlI/PTteG1D1MHTz0jhmNGPMPRw+mBkwbZ4+jQ/P18+kzOlhNbWVnR0dJSsLPQV9+rHtXihi6mpKYyNjeHatWuYn59HS0tLCb7xeLp4evvw8HDJKFnWlBafBRBsDUDNssT8/PtPSOfHs129erWG9fHF6sDiQ1F4DXN2drZkjbEUEMEj9LNp0ya0tbWhKApMTk5ieHi4fExb6DRSe3t7uX6+fv16bNmypbyiI+xTZw38rb/ZJ6I/6msakB0o8a51Dvw4r2OTKp/m45mX9qFeWjVAqcmxPh48BSrHNDU5Gs5lIo+2zYaj8kReduSuri60tbVh48aN2Lp1a/lA1tHRUUxPT2N0dBSjo6M4c+ZMjVHwfcLxBKOenp6aAed7xYPB8HVqPT092LRpU/kA3Zha6TRzYmKiBLzoR0w/w+HDcMPReQ2J1ytD5s7OTnR0dJQPmw3WUxRFOZ3lJQdmibyTy/e6x9iGXAxY69evL197EK9U4AeMRMCJ6evVq1cxMjJSvvoiAGt0dLTUYdQfVzj09vaio6Oj5v1JMU5xMwIzXnbuYH/A++ux58+fx7Vr19DZ2VkTiPg+fr04f+PGjUgpleMTutm8eXP5bqAAqPHxcXR0dGDDhg2YmJgogTzGcWBgoHzAb4C9Ix9s97yG62Y5YZtsh+pjDJI80+P/UVZBzE3N+b9O3yPx8g/bm6bcxjCnVQGUDuz0vCoU8HTZgaADSDetdlN03ZHmss4YYsoVjtTR0VHeAsnToomJifLJ2jGdiylrOF68uiGMPBwknDheERBTsJaWFoyMjJSXLAVQ8VPHo1/ct5Bpfn6+hq3EuIQzBuOId78AKGXVdTZl87rTG3qKc7ng5ZYCeDy1baD2/uZw9JCf33cT65HxWoV4EHBMRWMMQ0bgfSDlqSOAckMgZA2w082WnTt3luyQ1yh1DS4CUuiGN86CITL7j/GNB0AHM49lCJ6xxNjl7gaqWu7h3wycAeTugnfOw0SE2+ZZo06tq4iPa0uXBmLcGGNYp0qGXFoVQMnJrUnwf452mnhzpqrTuQjK53htTfPyArXWGwMe5SJ/gBmzKGZxsVOpd3PEoj5vXsT0K9aRYsobC/oNDQ3l5hGw6BBtbW01L4ZiQA/HbW5uLqewYcz8Cgng/ZdxscFzpOblAWYQygKUUUZZBXHWtRqzTtMY0NgBA0RbW1vR3NyMdevWYfv27SXoKMMJHTPQ6Zoc9y3AjM/z5iDrIoKNrqfxFFevgWS2DaAEOl1SCcYdYxYBj5d2WI9RNgIWT9HVfxRUmGWynDzt5eM81qxDBk0NMDrmfI7tiXWl7FPzsI2yrfADZlxayetqbwLwH+jQ9QD+zwD+ZOn4LgAnAPx2URRXlsp8DcCXAcwD+NdFUTxZ1QYbHkcA7iwbnA60OmHUafpS02YcU/B1jsz1smNxckbCjJTrYQAOgFLH5F1Hnea6vrKxR/uceNOI5YqyLBPnYQPkPDw91OUQoHZXmwOPc0IeA10PZFblmAwnnkIrO3NjFmuE6uzRtk6pIwWDY3nZ8dme3FSRbSLaUzldoGY9Vz1lPGYRwbarSAG35aawzHi5T+qzOoXmfrsLwLlttjnN49Ycc7qNpGSJ2aROwVfyBKGVvDPnCIA7lgRpBHAWwHew+IKxZ4qi+P2U0leX/v9eSmkfgC8C2I/Fd+b8MKW0t6jzOgh1Io7ODIy59Y8l+ep22DmZgqQyB65bDZgTMwctE+fV0ID3Nx+UyfAF3OFgOi1h3QXo6nmWG6g1aJ6+KZsPB2HQUBbATuv0rFMbpzc3bhrAcpdw6DjxWOoTaJT5ajCIMgEsrDOVk+2HGVAktxHBfWIgibxs41wPy6A2wAFIj3N5loFnJyxPzq9UZ3oumKAjBCqPMmxNLB+X4766IOMAnnWmwSv8K2Zx9dIHnXo/CuC9oihOppQeB/Dw0vFvAngWwO8BeBzAt4qimAZwPKV0FMC9AH5WVbFjSpEUEOK3giiDioKVAzCNOrk22XjZ2FVWBVN+4rU6nWNE2s84x0CqUZXBj9/F4ozVBQMGLxetmcUzg3IgHHqqYvg5nStoOLlVZzrN0nLKolxA5PK5GYTLz20os2QGrSDnbIDLsY3rvd38rYDlbEL16oCZd5r1ulrXhvaJ6+C82n7k1WPcrvoqT7Vd3Zr0gnbua/zm805PufRBgfKLAP790u/NRVEMLjUymFLatHR8G4DnqcyZpWM1KaX0FQBfAYCBgQEAy51bo7oDpfjNTpTrtObRfDmGqe3xfwVQbS+YlvZLwYn7wfWpTDmAVUCLlGsnvpXhMRNyzFGdTY1N9erYlIKW6pWnii4gVgVBbtetPaoNOEDQtll/2h8+Fvmq7C8+zEDjPx/n8+6hEyxH5OU1btV/5M3dA+2O6TKO6se14wIM53GsWXUXbC+SbpixnLoUFPVwnQq2OburSit+FURafAPjrwH4X+plNceWWU5RFN8oiuKeoijuWb9+fekgdL7GSWPwdMDY+HNKiLymT8t+V7EGlsWVc4vJLhoyW1TnjymugoAyaS6roBoG5Ka5TgcsJ7OGKvblgFwBTfWl48PMy6UqkGQ96MYRj0uuv5qf++mAi9tTO8z1SUFRrxOM9nizJr5jk48vYFdgjaso+IoGRwIC9OLaTZYjR0b0SUM6U+N6FcRdgHKBhW2b/ZfXRDlYMGCy/pVta+LzSjYUc3LpgzDKzwJ4tSiKoaX/QymlgWKRTQ4AuLB0/AyAHVRuO4BzK2mAGYybYgO1TK4qUjgG6BzPTWV459mlHIBE4gHUvAo4kV/X/xgcObmpMLfrZNNoq8buZFcZncwrYWKOkeYARhmKY8GRl/Wk4J5SqtnB1cCiNqXy66VAOiVm/SvYu+PO5lyw5//uzpWQjet365LBHDmPslJNOV1o3dG2LiOw7Ay27AO8xsgfvl6Xy+dmQVF3jhmyHTFIullXFThy+iAvF/sdvD/tBoDvAvjS0u8vAfhLOv7FlFJrSmk3gD0AXqxXuSoTWD6lUiBUB9B8zul5pxbwTw1RgFKKX8WA+JtZnWNOVXU5w1UjyIEY182grWyLj+nOrfbFMYJc36vGwzEfZ7yOFSn4OEbq+sy6YxBy9elGQjg6BzK3KcRBTo9zv1ge1z8FGwVEPu42ITjAMEPNAXIwVtWByqYbVk4mHW+WQe3X+bYGoNixB2p9Uq9+AN7fhFNwduXVPkOeqrQiRplS6gDwGID/PR3+fQBPpJS+DOAUgN9aavCtlNITAA4BmAPwu0WdHe8QVKO8o/uspKp1B2WMfCwHvtqOA1oXiXLMkmXTgVXQ0LsIWHZ1vtwalLYH5N97rJFe+5lLDoBXCpwONHTXne1AH5rhAoELnqwjBxA6jgoU2g/dJMgFcw0qORvmcgraUb6hoaHm9kYFf9YH53FAFnXq2ORYqp5XUHMpdK6BNupRNh7t8szQjT2DJoOt2i8HtNwlVkqQdCZQlVYElEVRTADYIMeGsbgL7vJ/HcDXV1I3J+4YO1gO/dXgOSrVa4frYaVzvTwwDFzsfLym6Jw3kjIJjdiOYbCsetz1xYGoLnC7gJTTiwM3ngLxf3ZybkudR8e1HvA7HfG48dTUsWOWladrCgzq2Dy2DEK6s81lou+OxbKO4jeDJNsUA4ROp1XuOMZreA7YVhIsFCRVJ1of/9bprgZ5F+TCb9xTrByY6/i7YKizQV5GYUbL4PuRuYUReP+6pvgN1G7xc7RUZcbAsyHngIfrzTk9n4uyORBzIOpAl/uZA3JlXlXrpNxn1y637ZgF6y7nYGyATg6ug/OrXrldzafjxvW56Rg7Tk4OthPtv5bT85zcDjGDpmNxLhDxJ/IpK81Nu52OWX4FbC3HGzduRhX53S2UqhsOAJyHx5HL5PwmF2BTSsueXcC65/b0nTuRR+vPrXu6y46qCNaqAUqg9vYqPa7OyCwu8rg1QadEHWxWooJVRDwFYI5ObAzMhlfSNz5XBUzc59x0NxIzG046xYljWpbbVJnVifW+au4TsykNVDkW4IA69z8SX0AdbfP4uH7xxwUbltHZn55jO1Mb0dskWe56U10FKh5XbkNZnGOqrCsGd64/6lUQ4aUh1i8zWtYHzxqYMaq+eNoMLL/zS6fmrFuWaSVBnXe4VV9VIAmsIqBkJajzaiRy09c4z0p3TpZjXQ4gAD/lyEVDoPYujVwE1z7xtI1ZqvaB63TnuT5mQxw8eDdY++Pkq5cce+T/VWBYjzG7MVM9qyzRLt9x4YKWuyND++KYqQN+7ZcCqYJjfGt+rsOxTpVZxzny6JsVFfQc8+Pj2jd3OZGOgwuWDORanvUQs8kgOmqj0Q73me9ay42fm5Gw7qr8SNOqAcpIVetVOghALZXWMi5/5NX/ytLUSR0g6bmUUs1FuiGfA2KNouqU+u2Sk7FKR1G/BpOQ0e0mcjk9roaW01G0yXp3RqrBIz58Z5LbnNJ1Q8fSmGFyf91UjMsx2ITeArCUdXJZBT8eR9aBTuWZ8bm+8HjprYj8ABUOqCqfk9ElbqPKDjk/583ZkbLZYKYBgJGUKccx157avbvEj2dunDel9A9/KMY/VdIpSw6ItJPA+8xEnVadn8twu1w/l+ekxp4DYTaEaJ+dUCMZg0iuPQWyOB8PPdAdxXAanc5w21GOncEZEhu9/mZZ45iL8NxPBWhl5wwYWqc6ux7jcXeOqpsGoRt35wsDCYMVT5fjHI+xXp6iYMXgFh9lUAqSLoCrzG5jiP0otwmlY+HGMco529d6InGfOX88/i90ryCuAcy1yYxUl9v4maUcpBVf1Nfr3e+9qoCSvyOxU/KaSOStos8MLM7htC0HEhz12eCcsUT5MAA35cita+ox7rd+az91Ss7rOtzfqNMtlvO0n5MCJMuhSwzqFMqo2FjjuOqHyzl5colZmepMmZked7rIAZX2T+sJeRnIHNi5tT0+z2WYQXL7Wpb1zADv2HH0xbHNKubobEH9JqbPOv5A7ROPory+z4h16PQNoObSKS4Tds+XnYV8/IAX7s9K0qoBSk4O/Bz15zwKUlzWOWwkdVquO/67ayq5PidfjsUyoKjz6PlcX1ku1ZcCQiR1jlwKHVad1zUpJ3NEfdU5y6H1xHllpblor8w2gFUZlSunAMKg76bgVXXoOOsxzccOHh/3Lmouz+PO488Mle3KgR/LEfrRjTC2syp5dDxZ3xqwuN9sL5xSev8hyRroOY+Oee7SMD0fdu1mQtzvXFoVQKmAF8kBjesYG446BxuPYy1chwMnHRRXF9epbDLXX9c3p4/cQDpdcZ3OwXP6VBlyoKoOG3kcU1SGpO063eVALZe0rLKacHZmVi6gON3FMb0O0k3JFQT1W2cMUYafs5kDLAZvBR21ZwfY8Z/BKTc+OV3mAjuDPLM2BsRcG/pbZ0O55HxUWSovhTBjjSUOfQfVSmYuqwIodZCVDXAkUcPVAYn8UReXVWdhVsNlXT63zqey60JzlXG46X+OnbI83K5G0XrGr310DqF53DWWOdDW8WJD1DHL7Ujrbn0O8KvsRfXD+ghn5rI6xQ754mVemgLk+K2EKpsD1Gib19d4hzoH3tyuWw7Qso6FM9iFPSu7z60LuxkeJ9040fHnJTMdF/U3tpfc+Efi2YOyUPbZkLnqwvKqdoBVApSAnwJHZ3V6l9vhdoCpYMmJlapTltwlSnosN01jg8vJnTMs1YW2ze06QHCyVk0tcmyE5dL1Vq5XZa5iEyEr36LGdbindqsjuTU51QUHQZbNbZzwGDDg8DuAeJ1QbSRkiHU5BTTtjwZoZoeR113iw/UGG9XbK7XfDuTcRfcuqCsI8bmcL/G4cN1ajyZ9iK5u6Ci4OzmY3Tq/VJbuMCOXVg1QsmKBWiPl7xx7BGqnOQqSynYiPzupYyacGPzcWqjWE8dUBh5EZarKZKMud10kt62/Wd7cf43srp7I465tY8NTORRMGSBdUFLA1XYc81TmwdNZDVRx0beCLy+n8HluR0GAA1/I6x53xsDG7eonjrv7vnPBUmdW0abaH8us46TndFz5280AFLRyx7gPrh4dW11LdGDGx3TZwvlZ1MtBSbGiKq0aoFzJ4/edAlk5Cjpcjy76VrGe+K/OEHXHMS2vgOxAT8FbwToHZG4Xm4+7QKA61KmdgrgaGh9XvefkVJ1zmSpj1Db0mkfOp3K53V9lgVo2/muZaDPnfI5ZRz2O3UX73D9uh+VUIGHwVN3yOa2Xf3OQcuugDoT1cjHHyNyYq53x1RXxyV1Vofdks76d3bvkfMhdFsVysgxVadUApTpYlUPnIkxu8LhcLlI55qmD6oA16uQ1GF7X4msd1ZBy8kTdvPajg+umWZy4Db5Ugh0r2nIgqf1UZ1Cw0vHSaz6dDnVNKeeIOgtQRudkUxthEFI9qT41HwOqHuf24+2VzN4UMBUAub0cOPPvsDPesGD5oj1ltKEvnea7gMfXlioQ8+28TkZun2VgQqDXnMa76jloKDiGPUX/Q9dxTsclpVTzahTtOx9bSVo1QOmYFSte0Z8/zMwc+9E2IrEB8SBp3Qqi7pgyx8inSYFFp/GOdbgBZcd3Gy7cti5JxFqO5nWgyM5fxcBX0qbW7aY/zkm4/tx4xNipLhiI2Ln4vAJaJHcPucoY59n5uF6Wz91uyP3QTQ8FDpaNp/mcGBxz+lS7jhTtM1PNXXPL/eA64zf3O5IG0FzwdsDK5dUunD9F+zy2uUuw2C5yadUAJQ9cbmobiYFM6+C6HGPLTVPUgTl6ssFVteucXaf8DqB4kNxaoJseqjE64+D8yjzdN3/Uiaoeme+WAgIctZ5IOrbKJJwcGkC5HWZrOpV148Ufdzwnq07hXJlwaL2LRwOJ++8+3GZOV8z0FEwUBFhfrDPuk+rdba4oCDIgh51HOWV9vAehty5GfTqLUt0q4OrSBwdOZ6NRx0pY5aoBSqB2vcIpSqOfRloGMhctFXzrrZvo+qReR8n5HNtw0wcAy4xTN65YvqKoncorw1bwr2J9VSCgfQpjZz0oqLigwbK7y2ZUXmVdmlfzOcatwOHqUKaVCxwMvlwXjyVfFsTTQScPB1iViYMo1+/6pP0LPbq1PbYNnZFoAGL9MLCwPbjx5zaifdUTgxX3R+25qamplFXr1fEP9seyRTl3PaXKq7KzD1elFd0fllL6L1NKb6WUDqaU/n1KqS2ltD6l9HRK6d2l7z7K/7WU0tGU0pGU0qdX2IYFrCrH5jx8e5KCmU4d2en4mG4+xHFVrK4b6hQqyrERcV53n25cl8dOEmVZZnUIBefIw6DE9TGwqxw6LVIwUZahulLdcllds+QPy+GYlgMO7Rs7f9Sp64AqM9sG51GncnXwMf7PNsNsV0FI9RUyq4yOEfF51j1fXsN1K+hzWR2vnC0riGlw1mCQC4whZywd5IK72p36msqvPqmECPCEhO0gF/SBFQBlSmkbgH8N4J6iKG4F0IjF19Z+FcAzRVHsAfDM0n+klPYtnd8P4DMA/jCltPyq3eXtZJ0/FOCcLX5zJ7XD7n/Up5em5Jisc0ZVMA8I18nGpIagwKPH2UAjKZhwn5QpxDkXRaM97gPX4xijK89603ZzrNABPMsTl9vk1mGd7sIB2QkZ6PSNgbngrMHKga86ooK+Amv0RWdHXEbZncqi4MV1cL0KmDmADkaqfsZ2rvWz7rRuF8x03JW0cN9cu7my7hPl9f32TnYX6NWnOK305WJNANpTSk0AOrD4VsXHAXxz6fw3Afz60u/HAXyrKIrpoiiOAzgK4N56DShIxjF1Cu6Yq4O/uY4oF+DIAOkiWuRX8NTEAM9JpyJxzMmZAxrXP46oeumDOowOvHM67a/qQkHZjQ+XZcDQ8XMsWvum7Slz4bKRNwdayka1LU08lsr81ZnDcWdnZ2uAWYOc6iHkzwEvEwKWRYlE1K/gz/1TPemHL2oH/JtP3fXCmj+OsW/Nzy++blfH0pXRcQz9xpQ8UtStfeUyoRf1O52GKxbUm3rXXaMsiuJsSum/x+ILxCYBPFUUxVMppc1FUQwu5RlMKW1aKrINwPNUxZmlYzUppfQVAF8BgIGBgRrjY8Tn9QbOk5F1WRTj76amphoA07xLctXUpfSe69TfkZhJVEUr7q9jYq4/nLgs/3byqaPr9ZeqPyeLRnPtB5dza8IsR0PD8ifMsM74P+tdAc4xELYdFyycPhVcVG4GXLVBnVpzcpfw8JRdwULldgDMgKl94QDh9M79ieNONwyQASQBytpP7R+v3SpgxRIZ64LP59YVc8ntWEf9s7OzNYAY53K6qEp1gTItrj0+DmA3gKsA/peU0r+sKmKOLfPyoii+AeAbAHDrrbcWS8esY/E5dQSRtcbI1eHV2XMK0st1mCXldpfjGBugY0Nclo/rlFD/a14HjA7wVXfMQtW5XdBQAONjGnDcNEgdIufAWlbHjL+rpmdsGywbg1MueLl64jhfjZCzmxz48FgFy2KgdDaggOcAzZ1nG2T9sJ6rLoBXXfDzI1k3jr0FKLlAw/UrsFZdURFJwZfZZ5znujUx6OeAtSqtZNf7UwCOF0VxcUmw/wjg4wCGUkoDxSKbHABwYSn/GQA7qPx2LE7VK5MCB/D+IPCumDrqkkxlHfFfHZbb0FQFOKx4ZW7shHE+B5A5Gbn/bHxVTuvq475wXt0McPU5ENW8DHzxrQ9JZR3rOOYcwMmsjEodTFmHGyPHxLhv7vpAlUXbd7LyMcfQdfqvIMkBh4+5ujg5x1YmqfVGngAzlVkvhYtjWqcCZ9gJA2S0oww+tzmp8gYL1LHVy4hcsAdqn1cZZV0KplwvrWSN8hSA+1NKHWmxtUcBvA3guwC+tJTnSwD+cun3dwF8MaXUmlLaDWAPgBfrNZIDNKXKburkHCL+626YGre2pwOj1D3XXiTtA0c/zVO1Qx+J+8sbAipz1Muyq+4c8Cm4s4waOBhUQ7c6DYt2derl2CjrgC/zqGIE0d/Z2VnMzs7W7FSHfpwduE0hBpGcXfF4OSBSGUMOx+p4o0fHTvOyrtQedHdeGRrL4vQBvB8ocmtz7De8W+zW3llutqvIG0ss2pZupKpNqe0A/j31oVOdMbBtsw1zG0WxePmde0oUp5WsUb6QUvo2gFcBzAH4ORanzF0AnkgpfRmLYPpbS/nfSik9AeDQUv7fLYqi8jnrykJyYFY1BVDHdwOqrIOP5ZzTMSWWhSOntuPk1DZZbi7D7eect8qxlZUpYDpZ4hzrh9t3+ndAzJdq8doW9ykXUFQeBgEto/rRKazqiMdInSp+6yaRMi4FS66PzymT4Smjlg9mroFS+67gzn3TQMJ90lkQM1pmcBq8FGS13hhrHf/ocwAkb8qoPlX/IRPfmBGyxAvIVDesD96o1fdXaT9c+7m0ogvOi6L4twD+rRyexiK7dPm/DuDrK6k7kjPAMMyGhoaax3IpsAK14KJ03rXFhhFtaX51NMcEnXNz0ojJ9ehxvd9VQZAd2AF3OEWOVXEAcaDJeeN3TKe4fgUON3UKHUc+Z9RVYxTHWT6eIqketA4Nko4BaZ/Dsd3MhutxIOUCVq5PutPLAKFTdQXX3HqaBgnVA/sPy1zPPpWJ85Q9R2wUUMN+qqa4zr7V39y6Lv/n9uMSvebm5hq9pZTKaXluNuXSqrozhyNfJH6+XA6YclHJKSCOayTUwciBGScFJmd0CrwcEPh6Sz5e5Xg5g+ZzyhLqTRer2K7Wq0xJ2aHTj7bv2quacgdgc/5cgFJbcAGYwY6ZEZ9nWVlGYPm7atwF86F73RGfnZ21AU+BMlh4XKnB7ase2G5z+lD7jGNKGFSmAJhgvXpebd4txYTsuu4Yx3L2wHVp8I9jwVbZf7lPXF/c663TbF5rz6VVBZQKkqxIdXoFCT2WMwCXX6eUOdn0PztI7pIHHTAFPb6MRqfLOtVTtsMy5YBEQcEFlWhL138UILhf7DRcnxo3n3PMxbE3zqfXW3LfQ2eqBw0YIS/rQqf6nEIXfCmL1snOzsc4EMexekwoF4hX8kSd6LPqkW1bL9PJgRrXHflc32NNzz2Yg+tV9susVmceHGB42q0A62Yq7O/apvPx2dnZmtfTOh1pWjVAqUYUqM/OkIt87ISOXTGAMtAq+HE7VVGZ62f5+Zvzc5u6s6h5dOBdfa7/Loo65wuWwnWzwfF3yOfeeewcQJ9O7hblVc4AQ86j+nHTzSjHY6/LDgy8LLcGgKq+uWWQSNyXXBkFBg0G8VuZntq7/mbd6Hl+So6rU4OxtqFlo05dMmDwYhCN5MiFBhe9d5uBzdlQjvkxqWJmyWV5tuiWB6rSqgTK6KA6biQ+r0CoQKbGrExUF4f5WM4QIz8PsDLYGDgXcXN16W6cAoormzP2SOqE0R6fV+dVlpxjRNqmnlcG5ACcHc1NMVVPaicMRG6K7y5Q14DI/ddNGW1bx8XdIaTyhA5VHl4a0Ck356ma6uvDcZkchHxRjwNNnbFEnRxwqjYueex0ms3ycV6uk+WN8jE9Zhm5XbXXmIKzTcQ6JOtRgZ7P8TNbXVo1QJlb24ikhsj53LSO61X2pkxKmQcfc+CqZRRgtD8azRcW3l9c52DABqdMhNtgo2LjrKdfBySsB15rZR2pDKx7d77KkdTYOUBx/zgxCKiBA8tvh+Qx4va4j1yez6suGRz1dkSdEitQar1xnDfIFOidDFy/CwgshzJP1z+2NZdf+6NjoX3kxDIqODFoB7i5OnQpjMGW63PyqOwagN2yQEy9PxJAGUlZxcLCQs3UzwGTMhUGT10Ezq3lxDF1Zl2Lyzmgk6EK3F19OSDOAbPKzf2L5Iye+65TeMfiNZqzsfG45e6f1ot/VVdufFxfuH2nbw1UVWOkOlOQya0Ps54C9GLsGEid/vlyFQY2x0qr1hY58Vpo/OfxibXE+B95mFnnZm65FLLHJkhsqKhvhN+Gbrh+XqbhGSL3STd6mpqaaq5MYFsM4NXdeZZZ9xHcbnourRqgVFYF1K5p6bcDCtdRBZmq/Ow4ygJcfbk82o7mjW+3iZRjAdqnquSAVtvW9vQ/Oz+w/FKgyBNrV3q1gDJfHl8FrFz/9JwGNmZhPKbMXnj5QKfPqof46M6sAoBL7hyXq5o2shwhK59zbbMec8ExzrkAzL7EoK9jEL+ZITKYK2uM4/Pzy1/3qxtjwPLbF7k+ndJHfr2Ma25ursbOdaYTx1TXWuYjwShVSI1yHAlccvdha71qpMoQlFlxXTrtURakO3XA8uUENyg5dpVzTMdUVX6Vm2Vi2RSoHCtQhq8PIWAn4nar+sNOqsDJ8rpAVhSFfZgG657L8ZipkyoY6O6xAofqNNaU+ZwCmAKJgrnKHfnjMiK2K31gsNMBb+S46aRbs9QNMAcaOpPhfLm8mphd696BLs1wGX2KkPpOBGseQ53tKB7w+dwrIjitGqAE8nc+8G5WDjQ5OjqFc9LBzhkrs1vHfLSuSApa3Dct41iTgiRPk1zb/FuBNCe76tm9oiLyuT6EbnRKFeXddFbb5/4FeAKwtwCybqrAOICU9eYu49JAyXKsZK1PyzhQjd/6GDbOw+PKjz2L4zw9zwVwBWrVlU6vdf2QbzF0AMfluM86ltoPoHYazSDNTN/NpqIMH9cX9TkGzIQp6ohybOORcjNRTasKKHkROJKChIvaURao3RGP/HHe7UKroUVyLEcvY9GUczqNvAqEKn/VupiCO8uqAK0BQwHCTW1yyU1ldF0tjvM51XdOHq1Pp2n6IInQqwMxzcf2kgPg+NbxdVPw3PWmKrPLq+VUhijPQJILqBoMXd+iTp0tBFPTseENFh6b0EluF12PMdFw5EdnElEv32CiPs2EievWdhko3Uwq/jNW1APLVQOUTumc+OJfPR8GEOfdO2Y0Aun0JepX2s71c7vcPt/bzPlz00g2WgVWpxOe3jH7yoGRMu8qHWhbjokrSHJ+basoimXrlqpDdlYH0rq4ruDjACYSy+GWGZzO3AzFBU4F6yijwK75VT4HGu4aPwdAcVzZsGP+um6sfdJHx7lxZjm0fZ7Gqv1wfbp8w0nvlgnbdrLkArr6kTvOdfAGWG4TUlNaCe38RaeU0iiAIx+2HP9IaSOASx+2EP8Iaa0fqyut9eMXn64riqLfnVgtjPJIURT3fNhC/GOklNLL/2voy1o/Vlda68eHm1b6zpy1tJbW0lr6/9u0BpRraS2tpbVUJ60WoPzGhy3AP2L6X0tf1vqxutJaPz7EtCo2c9bSWlpLa2k1p9XCKNfSWlpLa2nVpjWgXEtraS2tpTrpQwfKlNJnUkpHUkpHU0pf/bDlqUoppR0ppb9NKb2dUnorpfRvlo6vTyk9nVJ6d+m7j8p8balvR1JKn/7wpF+eUkqNKaWfp5S+t/T/I9ePlFJvSunbKaXDS+PywEe0H//lkk0dTCn9+5RS20elHymlP04pXUgpHaRjH1j2lNLdKaU3l87931PuCvMPI7k7Hv6pPgAaAbwH4HoALQBeB7Dvw5SpjrwDAO5a+t0N4B0A+wD8XwF8den4VwH8d0u/9y31qRXA7qW+Nn7Y/aD+/FcA/gzA95b+f+T6AeCbAP53S79bAPR+1PoBYBuA4wDal/4/AeA/+6j0A8BDAO4CcJCOfWDZsfha6wcAJAB/A+CzH/bYxOfDZpT3AjhaFMWxoihmAHwLwOMfskzZVBTFYFEUry79HsXi+823YVHmby5l+yaAX1/6/TiAbxVFMV0UxXEAR7HY5w89pZS2A/hVAH9Ehz9S/Ugp9WDRSf8dABRFMVMUxVV8xPqxlJoAtKeUmgB0ADiHj0g/iqJ4DsBlOfyBZE8pDQDoKYriZ8Uiav4JlfnQ04cNlNsAnKb/Z5aOrfqUUtoF4E4ALwDYXBTFILAIpgA2LWVbzf37HwD81wD4htqPWj+uB3ARwP+0tITwRymlTnzE+lEUxVkA/z2AUwAGAVwriuIpfMT6IemDyr5t6bceXxXpwwZKtwax6q9XSil1AfhzAP9FURQjVVnNsQ+9fymlzwG4UBTFKystYo596P3AIgu7C8D/qyiKOwGMY3Gal0ursh9L63ePY3EquhVAZ0rpX1YVMcc+9H6sMOVkX9V9+rCB8gyAHfR/OxanHKs2pZSasQiSf1oUxX9cOjy0NHXA0veFpeOrtX8PAvi1lNIJLC53PJJS+p/x0evHGQBniqJ4Yen/t7EInB+1fnwKwPGiKC4WRTEL4D8C+Dg+ev3g9EFlP7P0W4+vivRhA+VLAPaklHanlFoAfBHAdz9kmbJpaRfu3wF4uyiKP6BT3wXwpaXfXwLwl3T8iyml1pTSbgB7sLhg/aGmoii+VhTF9qIodmFR5/+pKIp/iY9eP84DOJ1Sumnp0KMADuEj1g8sTrnvTyl1LNnYo1hc//6o9YPTB5J9aXo+mlK6f0kH/1sq8+GnD3s3CcCvYHH3+D0A/+2HLU8dWX8Ji9OBNwC8tvT5FQAbADwD4N2l7/VU5r9d6tsRrKJdPJLvYby/6/2R6weAOwC8vDQmfwGg7yPaj/8LgMMADgL4/2JxV/gj0Q8A/x6La6uzWGSGX/77yA7gnqX+vwfg/4GlOwdXw2ftFsa1tJbW0lqqkz7sqfdaWktraS2t+rQGlGtpLa2ltVQnrQHlWlpLa2kt1UlrQLmW1tJaWkt10hpQrqW1tJbWUp20BpRraS2tpbVUJ60B5VpaS2tpLdVJ/z/m2HPeb+7EiQAAAABJRU5ErkJggg==",
+ "text/plain": [
+ "