code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from typing import NamedTuple, Tuple\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.patches import Circle, Arc\n",
    "plt.style.use('bmh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GaussianMixture(NamedTuple):\n",
    "    \"\"\"Tuple holding a gaussian mixture\"\"\"\n",
    "    mu: np.ndarray  # (K, d) array - each row corresponds to a gaussian component mean\n",
    "    var: np.ndarray  # (K, ) array - each row corresponds to the variance of a component\n",
    "    p: np.ndarray  # (K, ) array = each row corresponds to the weight of a component\n",
    "\n",
    "\n",
    "def init(X: np.ndarray, K: int,\n",
    "         seed: int = 0) -> Tuple[GaussianMixture, np.ndarray]:\n",
    "    \"\"\"Initializes the mixture model with random points as initial\n",
    "    means and uniform assingments\n",
    "\n",
    "    Args:\n",
    "        X: (n, d) array holding the data\n",
    "        K: number of components\n",
    "        seed: random seed\n",
    "\n",
    "    Returns:\n",
    "        mixture: the initialized gaussian mixture\n",
    "        post: (n, K) array holding the soft counts\n",
    "            for all components for all examples\n",
    "\n",
    "    \"\"\"\n",
    "    np.random.seed(seed)\n",
    "    n, _ = X.shape\n",
    "    p = np.ones(K) / K\n",
    "\n",
    "    # select K random points as initial means\n",
    "    mu = X[np.random.choice(n, K, replace=False)]\n",
    "    var = np.zeros(K)\n",
    "    # Compute variance\n",
    "    for j in range(K):\n",
    "        var[j] = ((X - mu[j])**2).mean()\n",
    "\n",
    "    mixture = GaussianMixture(mu, var, p)\n",
    "    post = np.ones((n, K)) / K\n",
    "\n",
    "    return mixture, post\n",
    "\n",
    "def logsumexp(ns):\n",
    "    mx = np.max(ns)\n",
    "    ds = ns - mx\n",
    "    sumOfExp = np.exp(ds).sum()\n",
    "    return mx + np.log(sumOfExp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estep(X: np.ndarray, mixture: GaussianMixture) -> Tuple[np.ndarray, float]:\n",
    "    \"\"\"E-step: Softly assigns each datapoint to each of the Gaussian distributions\n",
    "\n",
    "    Args:\n",
    "        X: (n, d) array holding the data, with incomplete entries (set to 0)\n",
    "        mixture: the current gaussian mixture\n",
    "\n",
    "    Returns:\n",
    "        np.ndarray: (n, K) array holding the soft counts\n",
    "            for all components for all examples\n",
    "        float: log-likelihood of the assignment\n",
    "\n",
    "    \"\"\"\n",
    "    n, _ = X.shape\n",
    "    K, _ = mixture.mu.shape\n",
    "    post = np.zeros((n, K))\n",
    "\n",
    "    ll = 0\n",
    "    for i in range(n):\n",
    "        mask = (X[i, :] != 0)\n",
    "        for j in range(K):\n",
    "            log_likelihood = log_gaussian(X[i, mask], mixture.mu[j, mask],\n",
    "                                          mixture.var[j])\n",
    "            post[i, j] = np.log(mixture.p[j] + 1e-16) + log_likelihood\n",
    "        total = logsumexp(post[i, :])\n",
    "        post[i, :] = post[i, :] - total\n",
    "        ll += total\n",
    "\n",
    "    return np.exp(post), ll\n",
    "\n",
    "\n",
    "def log_gaussian(x: np.ndarray, mean: np.ndarray, var: float) -> float:\n",
    "    \"\"\"Computes the log probablity of vector x under a normal distribution\n",
    "\n",
    "    Args:\n",
    "        x: (d, ) array holding the vector's coordinates\n",
    "        mean: (d, ) mean of the gaussian\n",
    "        var: variance of the gaussian\n",
    "\n",
    "    Returns:\n",
    "        float: the log probability\n",
    "    \"\"\"\n",
    "    d = len(x)\n",
    "    log_prob = -d / 2.0 * np.log(2 * np.pi * var)\n",
    "    log_prob -= 0.5 * ((x - mean)**2).sum() / var\n",
    "    return log_prob\n",
    "\n",
    "\n",
    "\n",
    "def mstep(X: np.ndarray, post: np.ndarray, mixture: GaussianMixture,\n",
    "          min_variance: float = .25) -> GaussianMixture:\n",
    "    \"\"\"M-step: Updates the Gaussian statistics (p,mu,var) for better fit\n",
    "    of the weighted dataset\n",
    "\n",
    "    Args:\n",
    "        X: (n, d) array holding the data, with incomplete entries (set to 0)\n",
    "        post: (n, K) array holding the soft counts\n",
    "            for all components for all examples\n",
    "        mixture: the current gaussian mixture\n",
    "        min_variance: the minimum variance for each gaussian\n",
    "\n",
    "    Returns:\n",
    "        GaussianMixture: the new gaussian mixture\n",
    "    \"\"\"\n",
    "    n, d = X.shape\n",
    "    _, K = post.shape\n",
    "\n",
    "    n_hat = post.sum(axis=0)\n",
    "    p = n_hat / n\n",
    "\n",
    "    mu = mixture.mu.copy()\n",
    "    var = np.zeros(K)\n",
    "\n",
    "    for j in range(K):\n",
    "        sse, weight = 0, 0\n",
    "        for l in range(d):\n",
    "            mask = (X[:, l] != 0)\n",
    "            n_sum = post[mask, j].sum()\n",
    "            if (n_sum >= 1):\n",
    "                # Updating mean\n",
    "                mu[j, l] = (X[mask, l] @ post[mask, j]) / n_sum\n",
    "            # Computing variance\n",
    "            sse += ((mu[j, l] - X[mask, l])**2) @ post[mask, j]\n",
    "            weight += n_sum\n",
    "        var[j] = sse / weight\n",
    "        if var[j] < min_variance:\n",
    "            var[j] = min_variance\n",
    "\n",
    "    return GaussianMixture(mu, var, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "''' The main algorithm to estimate the parameters/statistics of the Gaussians'''\n",
    "def run(X: np.ndarray, mixture: GaussianMixture,\n",
    "        post: np.ndarray) -> Tuple[GaussianMixture, np.ndarray, float]:\n",
    "    \"\"\"Runs the main algorithm. Iteratively perform E-step \n",
    "                                and M-step until convergence\n",
    "\n",
    "    Args:\n",
    "        X: (n, d) array holding the data\n",
    "        post: (n, K) array holding the soft counts\n",
    "            for all components for all examples\n",
    "\n",
    "    Returns:\n",
    "        GaussianMixture: the new gaussian mixture\n",
    "        np.ndarray: (n, K) array holding the soft counts\n",
    "            for all components for all examples\n",
    "        float: log-likelihood of the current assignment\n",
    "    \"\"\"\n",
    "\n",
    "    prev_ll = None\n",
    "    ll = None\n",
    "    \n",
    "    while (prev_ll is None or ll - prev_ll > 1e-6 * np.abs(ll)):\n",
    "        ''' This algorithm will iteratively run until convergence \n",
    "            The convergence is measured by relative improvement in the \n",
    "            log-likelihood measure of the Gaussians\n",
    "        '''\n",
    "        prev_ll = ll \n",
    "        post, ll = estep(X, mixture) # Perform E-step\n",
    "        mixture = mstep(X, post, mixture) # Perform M-step\n",
    "\n",
    "    return mixture, post, ll"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Predicting the incomplete entries from the calculated Statistics:\n",
    "\n",
    "Let $G_i$ be the $\\mathrm{i^{th}-Gaussian}$ distribution, for $1 \\leq i \\leq K$. Then prediction for the $n^{th}$ missing element is: \n",
    "\n",
    "$$ pred_n = \\sum_{i=1}^{K} \\mathbb{P}(X_n \\in G_i) \\times Mean(G_i)$$\n",
    "\n",
    "\n",
    "where $\\mathbb{P}(X_n \\in G_i)$ is probability that $X_n$ was sampled from $G_i$, and $Mean(G_i)$ is the mean of the $\\mathrm{i^{th}-Gaussian}$. \n",
    "\n",
    "Below code implements this formula to fill each of the missing matrix entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fill_matrix(X: np.ndarray, mixture: GaussianMixture) -> np.ndarray:\n",
    "    \"\"\"Fills an incomplete matrix according to \n",
    "    previously estimated Gaussian mixture statistics \n",
    "\n",
    "    Args:\n",
    "        X: (n, d) array of incomplete data (incomplete entries =0)\n",
    "        mixture: a mixture of gaussians\n",
    "\n",
    "    Returns\n",
    "        np.ndarray: a (n, d) array with completed data\n",
    "    \"\"\"\n",
    "    n, d = X.shape\n",
    "    X_pred = X.copy()\n",
    "    K, _ = mixture.mu.shape\n",
    "\n",
    "    for i in range(n):\n",
    "        mask = X[i, :] != 0\n",
    "        mask0 = X[i, :] == 0\n",
    "        post = np.zeros(K)\n",
    "        for j in range(K):\n",
    "            log_likelihood = log_gaussian(X[i, mask], mixture.mu[j, mask],\n",
    "                                          mixture.var[j])\n",
    "            post[j] = np.log(mixture.p[j]) + log_likelihood\n",
    "        post = np.exp(post - logsumexp(post))\n",
    "        X_pred[i, mask0] = np.dot(post, mixture.mu[:, mask0])\n",
    "    return X_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(X: np.ndarray, mixture: GaussianMixture, post: np.ndarray,\n",
    "         title: str):\n",
    "    \"\"\"Plots the Gaussians and the data points for 2D data\"\"\"\n",
    "    _, K = post.shape\n",
    "\n",
    "    percent = post / post.sum(axis=1).reshape(-1, 1)\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.title.set_text(title)\n",
    "    ax.set_xlim((-20, 20))\n",
    "    ax.set_ylim((-20, 20))\n",
    "    r = 0.75\n",
    "    color = [\"r\", \"b\", \"k\", \"y\", \"m\", \"c\"]\n",
    "    for i, point in enumerate(X):\n",
    "        theta = 0\n",
    "        for j in range(K):\n",
    "            offset = percent[i, j] * 360\n",
    "            arc = Arc(point,\n",
    "                      r,\n",
    "                      r,\n",
    "                      0,\n",
    "                      theta,\n",
    "                      theta + offset,\n",
    "                      edgecolor=color[j])\n",
    "            ax.add_patch(arc)\n",
    "            theta += offset\n",
    "    for j in range(K):\n",
    "        mu = mixture.mu[j]\n",
    "        sigma = np.sqrt(mixture.var[j])\n",
    "        circle = Circle(mu, sigma, color=color[j], fill=False)\n",
    "        ax.add_patch(circle)\n",
    "        legend = \"mu = ({:0.2f}, {:0.2f})\\n stdv = {:0.2f}\".format(\n",
    "            mu[0], mu[1], sigma)\n",
    "        ax.text(mu[0], mu[1], legend)\n",
    "    plt.axis('equal')\n",
    "    plt.xlabel('X axis')\n",
    "    plt.ylabel('Y axis')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmse(X, Y):\n",
    "    ''' Root mean squared error'''\n",
    "    return np.sqrt(np.mean((X - Y)**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def missing(data):\n",
    "    ''' Measure the percentage of elements \n",
    "        missing in the data matrix \n",
    "    '''\n",
    "    return 100*(data==0).mean()\n",
    "\n",
    "def run_test():\n",
    "    '''\n",
    "    Estimate the statistics for:\n",
    "        Experimental dataset 1\n",
    "    '''\n",
    "    K = 5 # Number of clusters / classes\n",
    "    seed = 0 # Setting numpy random seed=0\n",
    "    mixture, post = init(X, K, seed)\n",
    "    mixture, post, ll = run(X, mixture, post)\n",
    "    title = \"Plot for Experimental set 1\"\n",
    "    plot(X,mixture,post,title) #PLotting the Gaussians and the data\n",
    "\n",
    "def run_matrix_completion():\n",
    "    '''Evaluating the technique on the Netflix data'''\n",
    "    K = 12\n",
    "    seed = 1\n",
    "    mixture, post = init(X, K, seed)\n",
    "    mixture, post, ll = run(X, mixture, post)\n",
    "    X_pred = fill_matrix(X, mixture)\n",
    "    X_gold = np.loadtxt('netflix_complete.txt')\n",
    "    print(\"RMSE Error:\", rmse(X_gold, X_pred))\n",
    "\n",
    "def Test():\n",
    "    ''' Evaluating the technique on \n",
    "        Experimental dataset 2 '''\n",
    "    X = np.loadtxt(\"test_incomplete.txt\")\n",
    "    X_gold = np.loadtxt(\"test_complete.txt\")\n",
    "    print(\"X_incomplete:\\n\" + str(X)) # Printing the incomplete data matrix\n",
    "                                      # 0's indicate missing entries\n",
    "    K = 5 # Number of Gaussians \n",
    "    n, d = X.shape\n",
    "    seed = 0 # Numpy random seed \n",
    "    print('Missing %:',np.round(missing(X),6))\n",
    "    mixture, _ = init(X, K, seed)\n",
    "    print()\n",
    "    post, ll = estep(X, mixture)\n",
    "\n",
    "    mu, var, p = mstep(X, post, mixture)\n",
    "    (mu, var, p), post, ll = run(X, mixture, post)\n",
    "\n",
    "    X_pred = fill_matrix(X, GaussianMixture(mu, var, p))\n",
    "            # Filling the matrix uing the above computed statistics \n",
    "        \n",
    "    error = rmse(X_gold, X_pred) # Evaluating RMSE error\n",
    "    print(\"X_complete:\\n\" + str(X_gold)) # Printing original complete \n",
    "                                         #  data matrix \n",
    "    print(\"X_pred:\\n\" + str(X_pred)) \n",
    "    print(\"RMSE Error: \" + str(error)) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading and testing on Experimental Dataset 1\n",
    "\n",
    "Experimental dataset 1 can be found in \"exp_data5.txt\" file. This data set is a set of 2D points (1000 such points) on the cartesian plane, hence can be represented by a 1000 x 2 matrix. It is also to visualize on the 2D cartesian and hence we will use this to demonstrate how the algorithm converges tying to estimate the statistics for the Gaussians.  \n",
    "\n",
    "This dataset was originally contructed by randomly sampling points from 5 different Gaussian distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEXCAYAAACqIS9uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydeXxNZ/7430/2RRLZ7CSqQTFIrQliSRVVrdZPqzU680VRWqN00apSptrR0UUrStFNtdOiY2y1llChVbHTIUgQWUUi601yn98f9+ZO1pubK8m9T3ver9d95Z5znvOc91nyOec+51mElBINDQ0NjT8ODrYW0NDQ0NCoX7TAr6GhofEHQwv8GhoaGn8wtMCvoaGh8QdDC/waGhoafzC0wK+hoaHxB0ML/H8QhBBSCPHnOsrbWwjxvRAiy7id4LrYjj0jhBhg3PcWtna5U4QQV4QQr9naQ6Pu0AL/7wAhxGfGoCOFEEVCiHghxMdCCP87yLNvDYL4M0AY0AdoCly1drtVuMwvtX/lPwG1ua074BCGfU+0tQiAEOLPQgi7baQjhNgthPjMgnRNhRBfCSHOGK/t3fWg97vHydYCGrXGAeAxDOe0G7AKaAkMr4dthwBnpJSn7iQTIYSLlFJXxeIrGG4u5Um/k23WBqW8k2zt8jvEFbgJvAuMRotZtYL2xP/7QSelTJJSXpNSbgLeB4YKIdwrS2x8kvpGCHFLCJEnhNgnhOhuXBaM4UYCcNn4ZL2vinyuABOAQaXTCSG8hBArhBCpQoh8IcRRIcT9pdYLNqYfK4TYJoTIARaZ2b9i4/6V/0ghhKsQIlYI8e9S+bsLIU4LIf5VbnvjhBB7jPt8WQgxttz+NDb+gkoVQtwWQvwkhIgotbykSGe4EOKgECIfmFS+qKfU9ANCiBjj9n4VQnQ0fg4KIXKFED8LITqUc+gmhNgphMg2emwUQgSVWj5fCHFRCPGwEOK8ECJHCPGjEKJNybaBL43fS34ZfWacHmw81zeFEJlCiP1CiJ5mjnsFhBDOQoh3hRDXhBAFQogbQohvyqUZI4Q4bjz3V4zpPY3LPgMigb+U8htQ2baklFeklM9JKVej3VhrDS3w/37Jw3B+KzwhCSEE8G+gPfAg0BNIBnYZi06uAg8bk/fEUITxaBXb6QF8i+FGUTrdGmAI8GcgFPgJ2CKEaF9u/X8A64A/ActqupMAUsoC4HEgUgjxrHH2UsADmFTJ9tYAXYGvgC9L3fDcgR8BL2CY0XsbhuNyT7l8lgCLgXswHMuqeBOYg+FXmA74GlgOzCs179OSxMabwH4gBugODAKKjQ5upfJtiqGIbSwQDjQ07hcYip2eLZWuKfA343QDDMe5t3G9C8APombFgs9h+HX5Zwy/9h4CDpfah78a93EJ0AF4CrgP+NiY5G8YrpdvS/kdqsH2Ne4UKaX2UfwDfAbsLjXdAYgDDpeaJ4E/G79HGqc7lFruCtwAXjdO9zWmCbZi+3cb132gXLpjwBrj92BjmrkW5D8f0APZ5T4nyqX7C5APLMAQUHuWWlayvYXl1jkErDV+/ytwDXAql2Yv8L7x+wBjPuPKpSmZ36Lc9MhSaUYb540qNe8R47wGpY7lN+XydgVyS/IyHo8iILBUmjHGY+RmnP6z4d+72mPrAGQAY0vNuwK8ZmadD4zHRFSx/Aowpdy8CON++hqndwOf3cl1rn2s/2jlZb8fBgghsgFHDIFiDzC5irQdgXQp5dmSGVLKAiHEEeOyO6Wk6CK63PxoKpbT/2xhnlcx3LBKU+Z9gJTycyHEA8BcYLaUsrK8Y8pN/1Qq3x5AE+CW4UeRCVcMv6Cs8T5R6ntJUcXJSuY1wnAz6wHcbTyXpXHD8HRdQqKUMrXU9HVAGPNJqEpGCNEaw40xzJjWAcMvo6Cq1qmET4FdwEUhxC7j981SSp0QItCY17tCiH+W3rTx793ALzXYlkYdoAX+3w9HMDzxFgE3pKH4wxyV1fgQVcyvLSrLP8fCdQullBfNZi5EA+BeDEUjbWvgVIIDcA7DU3h5cstNW+xd6rs0M8+h1N8vgbcryav0i+zyL8HL51MVW4A0YBqGm6kOOAi4VLPe/zYk5XHjDWQwMBDDL4CFQojepbb/NwzFZuW5Zul2NOoOLfD/fsirLjCW4gwQIIToUPLUL4RwxVCeH2VMUxJYHK1wOWP8G4GhjLyEfkCsFflZynIMQX8QhjLxXVLKb8ql6V3OKQxDsAc4iqE8OktKmVKHnuY4CnQG4qSxfMNKdABCCEcpZbHxuz+GX2MPSCl3GOe1wPDkXyOklNnA98D3QohFGIoJ+0spNwshrgLtpJSfVONnzbWlUQtogf+PyV4MRRXrhBDTgEwMxSNuGIInQDyGMuMHjDVjCqSUmZZkLqWME0J8B0QJISYb83oG6AQ8aaWzoxCiSSXz06SURcLQOG000Nv4RPoqsEIIcURKeblU+glCiPMYAuyfMQT+GcZlXwHPA1uFEHOA/wKNMdxIzkkpzb3ErS0WYTg3a4UQHwCpGN5PjAQ+kFJesjCfkn1+SAhxEENRVYYxv6eFEHGAP4YX1OWLscwihHgRQ3uF4xh+CT2B4Yb7X2OSOcBqIcQtDC++CzG8BB8mpSwpfrwMDDTWRMoEMqWUpX8Jld5eV+NXP6BBybSU8nhNvDX+h1ar5w+I8UlyJHAe2IqhzLUJMFhKmWZMkwy8AszG8DS3qYabmQjsANZiKOfuAzwopTxvpXaw0aP8p6sQ4m4Mv1ReLBUM3sXw4vYbIYRzqXxmY6jpcxLD0/1fpJS/AEgp84H+GG4Kn2IIZBsx/BKKt9K7Rkgpz2GobdMAw/E7C3wCuAO3apDPLxiKYD7GUGPrIymlHsPNsQ2G/f8MQ7XfGzXUzAJmYnhfcgpD0dgoKeVvxm1/iaHWz3AMN7FfMLyQvl4qjyUYipxOYLgZ9TGzvVjjZwTQq9S0hpWIO/s1qaGhBsLQNuEy0E9KedC2NhoatkV74tfQ0ND4g6EFfg0NDY0/GFpRj4aGhsYfDO2JX+MPgRBipBDideP3KUKIU8a+ZA6W7itHCFFsnH9cCPGfKvJqZewbJ1YIcdLYaKz0cm8hxHUhxEcWeFWZl3HZTiHEOSHEWVFJT6lCiCBh6HvopLEPnhbllldwEYY+mkLK56Xxx0F74tf4QyCEOAQ8JKVME0J4SymzjPMfAqZKKYcap7OllA2qyWslECulXG68aWyTUgaXWv4BEAjclFI+W0U21eYlDB3evSml3GVsnKaXUuaWW/87YIux1fIg4P+klOPMuQgh+mPovuNp80dN4/eK3Qf+ffv2SVdXV1traNwhN27c4Pnnn6dLly6cPn2akJAQhg8fzqpVq8jIyGD+/Pl06NCBVatW4eHhwZNPGqr7jx07ln/+8580bdrU6m0nJCSwePFiPvqo4gP4zp07+eGHH3j33XcBiIyMZM+ePWbz+8c//kGzZs0YN24cp06d4sMPP2TlypUAnD9/nq+++orevXtz/vx5Zs2aZVVely9f5h//+Acff/yx2fXHjh3Le++9R6NGjZBSMnjwYHbv3m3WRa/XM3r0aP71r3/h5KQ15fm9kpubmxYZGRlY2TK7P+uurq60b1++Q0fbEx8fT1BQTbo3sR324Orh4cH169dZt24d7du3JzIykiNHjrBv3z62b9/OunXr6NatGwEBAXh6eprOuaurK23atKFVq1Zl8hs/fjwXL1ZsqDx16lTGjBlTZt6vv/5K7969y1xHq1atIioqCp1Ox6ZNm2jTpg0AOp2OqVOn4ujoyIwZMxg+vOJwBm+99RYPPfQQ//73v8nNzeX777+nffv26PV6XnzxRZYvX050dDTJycnVXrtvvfUWo0aNqpBXXFwcTZo0YdGiRcTHx9O/f3/mzZuHo2PZxq733nsvZ8+eJSIigs2bN5Obm0ujRo1o2LAhL774Iq+99hpxcXEVXNq2bUtRURGdOnUy61ef2MN1agmqeB47dqzKtid2H/jtFWdn5+oT2Qn24hoUFESHDobi9Pbt29O/f3+EEHTo0IGEhIQaea5Zs6b6REaSk5MJCCg7UNfEiROZOHEi69evZ8mSJURFGXqqOHnyJE2bNuXKlSs8/PDDdOjQgdatW5dZd8OGDYwcOZJXX32Vn3/+mSlTpnDo0CFWr17N4MGDadHC8tEXN2zYwBNPPMGzzz5bJq+ioiJiYmLYv38/LVq0YPz48axbt45x48aVWX/BggW8/PLLfP3114SFhdG0aVOcnJxMLi1btiQuLq7CdgMDA0lKsq/u7e3lOq0OVTzNoQV+K/Hx8bG1gsXYi6uLy//6AXNwcDBNOzg4UFRUhI+PD05OTpQufiwoqLyvuZo88bu5uZGVlVVpPo8++miZ4piSIqXg4GD69u3LyZMnKwT+tWvX8sUXXwDQs2dPCgoKSE9P55dffiEmJobVq1eTk5ODTqfD09OTefPmVXlM1q5dy3fffVchr2bNmtG5c2eCg4MBGD58OEePHq2wftOmTU0u2dnZbN68GW9vb5PLqlWryM3NreCSn5+Pm5tbhfxsib1cp9Whiqc5tMBvJWlpaXh6etpawyJUcU1LS6NVq1bs2LEDgBMnThAfX/mv1Zo88bdt29YUXAHi4uJMRTs7d+40fb916xbu7u64urqSnp7OkSNHeO655yrk16JFC3744QdCQkL47bffKCgoICAgwFTOD7Bu3TqOHz9uCrQLFizg3nvv5cEHH6yQV3R0NE8++WSZvPz8/Lh16xZpaWkEBAQQHR1NaGhoBZf09HR8fX1xcHDg/fffZ+xYw4BiJS7x8fH89NNPZVxKjoG9FaGqdJ2q4GkOLfBbiUp3fVVcfXx8GDFiBN988w0RERGEhoaagvKdEB4ezty5c5FSIoTgk08+Yf/+/Tg7O9OwYUOWLTMM/PXbb78xc+ZMHBwc0Ov1/O1vfzMFx0WLFhEaGsqwYcNYuHAhzz33HN9++y1CCD766CPK9d9fgbNnzzJ06NAK8xcuXMiMGTNYvnx5mbwcHR1ZsGABI0eOREpJ165deeqppyq4HDx4kIULFyKEICwsjHfeeafCMS1PSkoK7u7uNGlSWZ93tkOl61R17L5WT0xMjLS3JxMw/PM0alTj3mxtgiqudek5e/Zshg4dyoABA2olv5q6jho1ig0bNtTKtmtCZZ5RUVF4eXlVeF9ga7TrtHY5duzYr5GRkd0rW6Y14LKSvLwa9WRrU1RxrUvPmTNn1mr+Nc3LFkEfKvf08fHhiSeesIGNebTrtP7QnvitpKCgAFXaF6jiqoonqOOqiieo46qKp/bEXwfYW1U4c9Snq16fR0rKh+Tn17zb/ZUrV3LjRuVdwx88eLBCbZ3a5O9//zudOnWiZcuWVab59ddfiYiIMH22bNkCQG5uLo8//ji9evUiLCyMN954o848a4p2ndY+qniaQwv8VlK6aqK9U5+uWVk7cHfvTlraqhqvu2XLFpv9Uw0ZMsTU4rUq7rnnHvbu3Ut0dDSffPIJM2fOpKioCIBnn32WI0eOsH//fo4cOcKuXbvqQ7tatOu09lHF0xxarR4r8fLysrWCxdSnq5fXIK5c+SsuLi3IytpFYeF18vLO4ef3GG5uHcjP/y8uLh3529/+RmxsLEIIxo4dS/PmzTlz5gyTJ0/Gzc2NHTt28NNPPzFnzhz8/Pzo0qULYOhuIDQ0lOjoaFPtim7durF9+/Y7euHWo0ePatN4eHiYvjs5OZlq8nh4eNCvXz/AEBQ6d+5MYmKi1S61iXad1j6qeJpDC/xWkp6eToMGZvvyshvq09XR0ZuGDR8hP/8Mjo4+5OWdwdGxAcXF2SQnv4uzc2OOHfuJxMREDh06BEBGRgq+vo346KOPWLx4MaGhoeTn5zNjxgw2bdrEXXfdxfjx4wFDY69hw4axZcsWxo4dy9GjR2nVqlWFoH/gwAHmzJlTwc/d3d3UTsAajh49ynPPPUdCQgIrVqyo0NdNZmYmO3bsYMqUKVZvozbRrtPaRxVPc2iB30p8fX1trWAx9e3q7/+/aoKurn/i5s3PyMk5hF6fS37+Re6+ezTx8at4+eWX6dkzl549PXBxGV0miF64cIGgoCBTPf7HHnuMzz//HIBHHnmEd955h7Fjx7Jx40ZGjhxZwaFfv35ER0fX+r51796dmJgYjh07xksvvcR9991nagFbVFTExIkTmTRpkqnFra3RrtPaRxVPc2iB30ry8vLw9va2tYZF2Mq1uDiLixeHASCED3r9TVxcWtK8eVeio6PZtet7vvxyGVu3wkcfPY5ery+zflWNonr27Mnly5dJS0tj27ZtvPDCCxXS1NUTfwnNmzfHw8ODc+fOmVrUzpgxgzZt2vDMM8/ccf61hXad1j6qeJpDC/xWkp+fb2sFi7GV65UrEykuzkSvz0evj8PFpTlCOBIbOwoPj9b06tUQH58ezJ9/HA+Pe3F3dyc7OxuAkJAQ4uPjuXz5Mq1bty5TD14IwfDhw3nttddo27Ytfn5+FbZdF0/88fHxNG/eHCcnJy5dusTFixdNvYa++eabZGVlsXTp0lrd5p2iXae1jyqe5tACv5XYW3N3c9jCtbAwlcLCawA4OLji6OiFm1sHGjd+kXPnpvLmm6eREoRw4Y03PgBg3LhxzJo1y/Ry97333mPMmDH4+fnRu3dvzp07Z8r/kUceITIy0tTdwp0yb9481q9fT25uLh07dmTcuHHMnj2b7du3Exsby6uvvsrhw4d5//33Tb0zvvPOO/j7+3P9+nWWLFlCSEiIqWXwxIkTTV0s2BLtOq19VPE0h9aAy0pU6ZMbbOOamDiP3NzjFBRcoqgoFWfnxri4BOHrOwoh3CkqSsHFJYiGDUfY1NNaVHFVxRPUcVXFU2vAVQfYW5e25rCFq6dnL1xcWqHX5+Ls3JTi4lsUFd0kOXkxHh6hCOFETs5hm3taiyquqniCOq6qeJpDC/xW4u7ubmsFi6lPVykl6emf4ejoR6tWH9Ku3UE8PcPx9h4OFOHr+2fAGUdHL7y977eZ552iiqsqnqCOqyqe5tACv5VkZGTYWsFi6tM1P/8s4ERm5mYAXFyaEhS0DAcHF9zc2uPu3h43t2D8/Mbi5dXfZp53iiquqniCOq6qeJpDC/xW4u/vb2sFi6lPVze39hQX3ywT1IuLc5BSh7NzU+OTf+Vox7T2UcUT1HFVxdMcWuC3ktu3b9tawWLq01UIRxo1ml6mGCcl5X0aNOiHp2dPHByq7udEO6a1jyqeoI6rKp7m0AK/leh0OlsrWIytXT08QsnLO427+5/MprO1Z01QxVUVT1DHVRVPc2j1+K1Epbq8tnb18XkAH58Hqk1na8+aoIqrKp6gjqsqnubQnvitRKU+uVVxVcUT1HFVxRPUcVXF0xxa4LcSlap0qeKqiieo46qKJ6jjqoqnObTAbyUqDcagiqsqnqCOqyqeoI6rKp7mqLfAL4R4XghxRghxWgjxtRDCTQjhJ4TYJYS4YPyrTH+nmZmZtlawGFVcVfEEdVxV8QR1XFXxNEe9BH4hRHNgOtBdStkJcATGALOBPVLKEGCPcVoJAgICbK1gMaq4quIJ6riq4gnquKriaY76LOpxAtyFEE6AB5AIPAx8blz+OVBxRA07RaW7viquqniCOq6qeII6rqp4mqNeqnNKKa8LIf4JJAB5wE4p5U4hRGMp5Q1jmhtCiAqDpqakpDBhwgScnJwoLi7m0UcfZdq0aSQlJeHp6YmjoyNZWVkEBgZy8+ZNpJQEBgaSnJxsGh4tOzubxo0bk5qaihACPz8/UlNT8fb2pri4mJycHJo0aUJSUhLOzs74+PiQlpaGj48POp2OvLw803IXFxe8vLxIS0ujQYMG5OXlkZ+fb1ru5uaGu7s7GRkZ+Pv7c/v2bXQ6nWm5u7s7Li4uZGZmEhAQQGZmJoWFhabldbFPt2/ftmif0tPT8fX1tdk+paWlmfajts5TXe1Tbm4u7u7uNrn2arJPaWlp+Pv72+zaq8k+FRQUkJycrP0/1dI+maNeumU2lt1vAB4HbgHfAeuBj6SUDUuly5BSlinnt9dumQsKCnB1dbW1hkWo4qqKJ6jjqoonqOOqiqc9dMt8H3BZSpkqpSwENgLhQLIQoimA8W9KPfncMSrV5VXFVRVPUMdVFU9Qx1UVT3PUV+BPAHoLITyEYSDVSOAc8B/gL8Y0fwE21ZPPHePp6WlrBYtRxVUVT1DHVRVPUMdVFU9z1FcZ/xEhxHrgGFAExAIrgQbAt0KICRhuDqPrw6c2cHR0tLWCxajiqoonqOOqiieo46qKpznqrVaPlHKelLK9lLKTlHKclLJASpkupYyUUoYY/96sL587JSsry9YKFqOKqyqeoI6rKp6gjqsqnubQWu5aSWBgoK0VLEYVV1U8weBaVFTEP/7xD7744gtb61SJasdUBVTxNIcW+K3k5k1lfpwo46qKJxhck5OTEUIQExPDnDlz2LJli621KqDaMVUBVTzNoQV+K6mParC1hSquqniCwbV58+b07NmTTp060bFjR06cOGFrrQqodkxVQBVPc2j98VuJSj/3VHFVxRP+5zpgwAD69+/P4cOHGTFihI2tKqLiMbV3VPE0h/bEbyXJycm2VrAYVVzt3VNKaRp2r7SrEIKwsDC8vLxspVYl33zzDYsXLzZNf//99/Tu3ZuwsDCefvrpMmmzsrLo2LEjL730UqV5vfrqq0RERBAREUGPHj0IDg42LZs/fz7h4eGEh4ezceNGi9zKu5Qc02vXrvHoo4/Sq1cvevfuTUJCQpV5bNq0CT8/P2JjY03zvv76a7p370737t35+uuvTfMnTJhAXFycRW7msPfr1BK0J34rsaRZtL2giqu9e77//vvcvHmT4cOHExISUmv53rx5E3d39zrp5/3LL7/k22+/BSAuLo7333+fH374gYYNG5Kamlom7aJFiwgPD68yr0WLFpm+r1y5kpMnTwKwc+dOTpw4QXR0NAUFBYwYMYL77rsPb2/vKvOqzMXBwfAc+swzzzBz5kwGDhxIdna2aX55bt++zcqVK+nWrZtpXkZGBosXL2bv3r0IIRg4cCDDhg2jYcOGjB8/nqVLl/LBBx9Uc9TMY+/XqSVoT/waGhbi4eFBfn4+bm5uFq+TkJBAr169mD59OuHh4UyaNIl9+/YxdOhQunfvzq5du3jrrbcYMWIEH374oWm98PBws0+6lnDx4kVcXFzw9/cH4IsvvmDChAk0bGjoJaV0kcXx48dJTU1l4MCBFuW9YcMGRo0aBcD58+fp06cPTk5OeHp60rFjR/bs2WN2/apczp8/T1FRkcmjQYMGeHh4VJrHokWLmD59epnzsXfvXgYMGICvry8NGzZkwIABJpewsDD2799PUVGRRfv4e0YL/FaSnZ1tawWLUcXV3j0nT57MO++8Q9euXWvkeunSJaZMmcLBgwe5cOEC69evZ/v27SxYsIBVq1ah1+st7vtl/PjxpuKW0p9vvvmmQtojR46U+WUSFxdHXFwcQ4cOZfDgwezevRsAvV7P3LlzeeONNyxyuHr1KgkJCURERADQqVMndu/eTW5uLunp6Rw8eJDr16+bzaMyl+zsbOLi4vDx8eGpp56if//+vP766xQXF1dY/+TJk1y/fp0hQ4aUmZ+YmEjz5s1N082aNSMxMREABwcHWrduzenTpy3az6qw9+vUErSiHitp3LixrRUsRhVXVTyhZq5BQUF06NABgPbt29O/f3+EEHTo0IEbN27wr3/9i7ffftuivNasWWPxdpOTk2nZsqVpuqioiEuXLrF582YSExN54IEHOHToEN9++y2DBw+mRYsWFuW7ceNGHnroIVML1kGDBhEbG8vQoUPx9/enR48eODmZDy2Vuezdu5eioiJiYmLYv38/LVq0YPz48axbt45x48aZ1tXr9cyZM4dly5ZVyLeyGjeGXmIMBAYG3nFfOypdp1WhBX4rSU1NLfNPZc+o4qqKJ9TMtfRQfQ4ODqZpBwcHU6CLjo4mLCzMlK6goKDSvMaPH8/FixcrzJ86dSpjxowpM8/NzY1r166Zpps1a0b37t1xdnYmKCiIkJAQ4uLi+OWXX4iJiWH16tXk5OSg0+nw9PRk3rx5lTps3LixzAtjgFmzZjFr1iwAnn76adq0aWPukFTqcuzYMZo1a0bnzp1NL46HDx/O0aNHy6ybnZ3NuXPnTLWoUlJSGDt2LF999RXNmzfn4MGDprSJiYn07dvXNF3TorrKUOk6rQot8FtJ6acIe0cVV1U8oXZdd+/eTa9evdi/fz9z587lxIkTxMfHU1RUxNKlSwkNDaVfv35AzZ7427ZtS0xMjGn6gQceYOPGjTz55JOkp6dz8eJFgoODWblypSnNunXrOH78eJVB/8KFC9y6dYuePXua5hUXF5OZmYmfnx9nzpzhzJkzpjL6BQsWcO+99/Lggw+Wyacyl5YtW3LPPfdw69Yt0tLSCAgIIDo6mtDQ0DLrent7l7n5jRgxggULFhAaGkpwcDALFy7k1q1bAPz444+8/vrrprRxcXHcaTfvKl2nVaGV8VuJn5+frRUsRhVXVTyhdl2nTJnC3XffjZeXFxEREaxZs4Y2bdpw5swZWrZsyY4dO6zKNzw8nAsXLpiKPyIjI/H19aV379489NBDvPHGG9Xux6JFi9i+fbtpesOGDTz66KNlgl9hYSEPPPAAvXv3ZsaMGaxYscJU1HP27FkaNaowvlKlLnfddReOjo4sWLCAkSNH0qdPHwCeeuqpSl0qw9fXlxdeeIHIyEgiIyN58cUX8fU1DPGRkpKCu7s7TZo0qe7QmUWl67Qq6mUgljvBXgdiiY+PJygoyNYaFqGKqyqeUNF12bJl3Lp1i9mzZ99x741ZWVkcPHiQgQMH8tlnnxEaGkrv3r2tymvatGmMHj2aAQMG3JGTtYwaNYoNGzZYlLauz39UVBReXl5l3hdYgyrXqbmBWLSiHisxV0fZ3lDFVRVPqOh6+fJlnJycWL16tal+fEn9c71eDxha+YaHh1dbVLBs2TJ69erFunXreOaZZ+7Ic8aMGZW+E6gvLA36UPfn38fHh8cff/yO81HpOq0KLfBbSWVVzFRvAvkAACAASURBVOwVVVxV8YSyrnq9niZNmhAfH0///v1p165dpel3797NK6+8wpgxY+jatWuVeQ8YMIBt27aZijjuBF9fX4YNG3bH+dQHdX3+x44dWyv5qHSdVoVWxm8lOTk5tlawGFVcVfGE/7lmZ2fz6quvMnDgQD788MNKgz4YBu8YMmQIixYt4vjx42Ve0q5bt44bN26YpsPCwli4cCEhISEcPHiwQm0dazyr48knnzTbavfMmTPcf//9hIWF0adPH/Lz88nNzeXxxx+nV69ehIWFWdwO4E5dbY0qnubQAr+V3OkLovpEFVdVPMHgWlhYyIIFCwgICGDz5s3k5+eblut0OtasWcOFCxfKrOfg4MBf//pXWrZsyWeffQYY+papq3FcLTmmmzdvNjucYFFREZMnT+bdd98lJiaGzZs34+zsDMCzzz7LkSNH2L9/P0eOHGHXrl116moPqOJpDi3wW4lKAy6r4mrPnnq9nuXLl7Nt2zbA4Prxxx/z9NNPk5KSgpOTE2lpaab0u3fvJigoqMwgLcXFxUybNo3w8HDmz5/Pnj17+Oijjzh+/DiTJ08mIiKCvLw8U/XOYcOGmfr41+v1dOnShczMTFN+3bp1IyUlxax3dcc0OzubqKgoUx38yvjxxx/p2LEjnTp1Agy1WhwdHfHw8DBVM3VxcaFz586mVrLWYM/nvzSqeJpDC/xWUvLEowKquNqzZ3x8PE5OTvz000+AIWAKIQgJCWH69OmMGDGiTMvXvn378vPPP5cpXz916hSJiYkcOnSIn376iQ8++IDr16/TtWtXVqxYQXR0NEIIZsyYwbp169i2bZupJ0gHB4cyN4KjR4/SqlWrClUlDxw4UKYrhzFjxhAREVGha4MSFi1axLRp06rsDwcMff4IIRg1ahQDBgxg6dKlFdJkZmayY8cO+vfvb+ERrYg9n//SqOJpDu3lrpX4+PjYWsFiVHG1Z8/g4GBcXV1NAXTfvn2ml69nzpwhOjqaW7dumQKft7c3r7zyCmDoyTI1NZXJkycTHx/Pyy+/zF133cX169cJCgri8OHDpu1cuHCBoKAgU8vXxx57jM8//xyARx55hHfeeYexY8eyceNGRo4cWcGzX79+REdHm6ZzcnKqLMY5deoUly9fZtGiRWY7hCsqKuLw4cPs2bMHd3d3Ro4cSZcuXUz7WlRUxMSJE5k0aVKZrpprij2f/9Ko4mkOLfBbSVpamtlyUXtCFVd79hRClKll4+Pjg5OTE3PmzOHKlStERkby1Vdf4ezsXOEl6YULF/Dy8kKv1xMdHc3evXtZsGABISEhNG7cuMJQflVV9+zZsyeXL18mLS2Nbdu28cILL1RIc+DAAebMmWOa1ul0uLi44O7uXqEh2C+//MKJEyfo0qULRUVFpKWlMWLECDZv3lwmXbNmzejTp4+pl8/Bgwdz4sQJU+CfMWMGbdq0ueOqp/Z8/kujiqc5tKIeK1Hprq+KqyqeN2/epLi4mJs3b+Ls7EyzZs1o3rw5/fr1Y+vWrRXST5o0iWHDhlFQUMC7776Lt7c37733HqdPn+b//u//cHZ2NvX4GBISQnx8PJcvXwbK1oMXQjB8+HBee+012rZtW2kL0pIn/pLPDz/8QHR0dKWtf8ePH8/Zs2c5ceIE27dvp02bNhWCPhha2Z45c4bc3FyKioo4dOiQqduDN998k6ysrDJ99VuLKudfFU9zaE/8VqLT6WytYDGquKriefXqVZo2bUqLFi144okn8PPzw8/PjzfeeAOdTkdxcbGp9W5MTAz/+c9/6NGjB6dPn+b777/nyy+/xM3Nja5du7J//346derErFmzcHNz44cffqBnz55ERkbSrl07evfuzblz50zbfuSRR4iMjKy0Z8rKsPaYbt++ndjYWF599VUaNmzI1KlTiYyMRAjB4MGDuf/++7l+/TpLliwhJCTE1DJ44sSJVrc/UOX8q+JpDi3wW0leXp6tFSxGFVdVPKWUpsBeut7+vn378PPzY9u2baaeI0v6h7927RpTp05FCEGvXr34+OOPuXr1Kjqdjnbt2vHJJ58QGxuLk5MTDRo04OGHH+btt9+u0E9/aGhohaIhc1h6TFu1asWhQ4dM08OGDSvzYvqxxx7jscceK7NO8+bNa+RSHaqcf1U8zaEFfitRqS6vKq4qecbGxpKVlVWm+X5wcDBXrlyhR48epnmjRo1iz549bNu2jaioKKZPnw6Ap6cnYWFhhIWFcfr0aZYvX06rVq349ddfeeWVV2o0OEt1rqqgiqsqnubQyvitRKW6vKq4quLZpEkTbty4waJFi/jtt98AQ8dqQUFBTJw4sUJgCAgIwMXFpUwRwf33309wcDD+/v4IIQgNDeXo0aPce++9NGnShGbNmtWKqyrHFNRxVcXTHNoTv5WUHlzD3lHFVRVPMPzcz83NNT2VJyQk4O/vX6Y8voQuXboQEhJSZjD1kjLxXbt20bdvX/r162dqDFWbqHRMVXFVxdMcWuC3Ei8vL1srWIwqrqp4AowePZoTJ07QrFkzlixZwp/+9Cf69OlDq1atKk1fWQMpKSU//vgjb775Zp15qnRMVXFVxdMcWlGPlaSnp9tawWJUcVXFEwxV+hISEkxjw+7fv58ePXrUaDzWjRs3MnTo0Dod0UmlY6qKqyqe5tACv5WUjOqjAqq4quIJBtfp06fzww8/kJ2dzejRo2u0/qlTp7hy5QoRERF1ZGhAtWOqAqp4mkML/FaiUpUuVVxV8QSDq5eXFy+++CLnz58nMDDQ4nX37t3Lf/7zH55//vkab1en05Gbm1sjT1VQxVUVT3Nogd9KSnfBa++o4movnlJKNmzYwLFjx6pMU+LapEkT3njjDdauXcuKFSvIysqqcp0LFy4wf/58bt++zZw5c0wjdNWEV155hZdffpnbt29blN5ejqklqOKqiqc5tJe7VqJSXV5VXO3F88KFC6Snp5uqV1ZGaVcPDw9efvll4uPjWblyJTqdDjc3N/z9/SkuLubYsWNcvHiRhx56iBdffLFCPy9RUVGkpqbyyiuvVFtjxNXVlaKiIovH9bWXY2oJqriq4mmOegv8QoiGwCqgEyCB8cBvwL+AYOAK8JiUMqO+nO6EpKQkJQZcBnVc7cXzrrvuYuvWrfTs2bPKNJW5BgUFmTpOy8vL4+bNmzg6OhIfH8/jjz9OYmJipZ17Xb58GQ8PDzIyMqp9OVzTPnHs5Zhagiquqniaoz6f+D8AfpBS/j8hhAvgAbwK7JFSvi2EmA3MBl6uRyercXNzs7WCxajiai+eTk5O1Za/V+b63nvvcfPmTebOnYu7uzvNmzcHDJ207d27t0KXByVMnz6dzMzMGtUIshR7OaaWoIqrKp7mqJcyfiGENxABrAaQUuqklLeAh4HPjck+Byp2MG6nlG6MY++o4mprT51Ox9q1a7l69Wq1aStzvXbtGkVFRRXK35s3b864ceOq7IKhZcuWptGtahtbH9OaoIqrKp7mqK+Xu3cBqcCnQohYIcQqIYQn0FhKeQPA+LeRuUzsiYwMJUqkAHVcbe25detWmjVrxqpVq6pNW95VSslTTz3FhAkTTP3WlycrK4uLFy/Wiqul2PqY1gRVXFXxNEd9FfU4AfcCz0kpjwghPsBQrFMtKSkpTJgwAScnJ4qLi3n00UeZNm0aSUlJeHp64ujoSFZWFoGBgdy8eRMpJYGBgSQnJ9OgQQPAMExe48aNSU1NRQiBn58fqampeHt7U1xcTE5ODk2aNCEpKQlnZ2d8fHxIS0vDx8cHnU5HXl6eabmLiwteXl7odDqysrLIy8sjPz/ftNzNzQ13d3cyMjLw9/fn9u3b6HQ603J3d3dcXFzIzMwkICCAzMxMCgsLTcvrYp/0ej05OTnV7lN6ejq+vr422yedTkd+fn6tnqea7NM999zDpk2bGDx4MCkpKWb3ydHRkYyMDNM+7dy5E51OR0FBAc7OzpWepx9++IGsrCwGDhxI69at62WfSpxsde3VZJ/c3NxITk7W/p9qaZ/MIaSUNQngViGEaAIcllIGG6f7YQj8dwMDpJQ3hBBNgX1Synal142JiZElgz7YEzdu3KBp06a21rAIVVxV8YSKrjt27ODAgQP06dOnTHfGpfn444+5cuUKL730UqWDqNSHpz2jiqsqnseOHfs1MjKye2XL6uWJX0qZJIS4KoRoJ6X8DYgEzho/fwHeNv7dVB8+tYFKgzGo4qqKJ1R0HTJkCP/973+Ji4urcp0pU6bUtVYFVD6m9ooqnuaoz1o9zwFfGWv0XAL+D8M7hm+FEBOABKBm7d5tiEp1eVVxtUfPgoICPv74Y3r37k2vXr1M88u76nQ6rl+/jpSSoqIinJzso4mMPR7TqlDFVRVPc9Rby10p5XEpZXcpZWcp5UgpZYaUMl1KGSmlDDH+rb3hfOoYlfrkVsXVHj0PHz7MPffcU2Es2vKuLi4ujBo1itGjR9tN0Af7PKZVoYqrKp7m0LpssBKVqnSp4mqPnmFhYfz222+MHFm2pnFlrj169KB790qLVG2GPR7TqlDFVRVPc9jPo4liqDQYgyqu9ujp4uLCc889V+l8FVDFE9RxVcXTHNoTv5VkZmbaWsFiVHFVxRPUcVXFE9RxVcXTHFrgt5KAgABbK1iMKq6qeG7dupW1a9eWmbdp0yb8/PyIjY0tMz8rK4uOHTvy0ksvmc2z/PoHDhwgIiLC9GnatClbt261yK90XgEBARbn9emnn9KnTx8iIiIYNmwY58+fNy2bP38+4eHhhIeHs3HjRtP8CRMmmK3JVBNUOf+qeJpDC/xWotJdXxVXVTyXLl1apsz/9u3brFy5km7dulVIu2jRIsLDw83mV9n6/fr1Izo6mujoaDZt2oS7uzsDBw6s1q18XpmZmRbnNWrUKH766Seio6OZPn06r732GgA7d+7kxIkTREdHs2vXLj766CNT99Pjx49n6dKl1XpZgirnXxVPc2iB30oKCwttrWAxqrha6pmQkECvXr2YPn064eHhTJo0iX379jF06FC6d+/Or7/+CsDbb7/Nhx9+aFovPDychISEO3K8ePEirq6uZVpHLlq0iOnTp1fovOv48eOkpqZWG7CrWr+ETZs2cd9991U6bm91eZU/puby8vb2Nn3Pzc01DQl5/vx5+vTpg5OTE56ennTs2JE9e/YAhpff+/fvp6ioqFq36vi9Xaf2jBb4rUSluryquNbE89KlS0yZMoWDBw9y4cIF1q9fz/bt21mwYAHvvfdejbY7fvz4MkUhJZ9vvvmmQtojR47QuXNnk+vJkye5fv06Q4YMKZNOr9czd+5c3njjDbPbrmr90nz//feMGjWq2v2oLK/yx7S6vFatWsW9997LvHnzePvttwHo1KkTu3fvJjc3l/T0dA4ePMj169cBcHBwoHXr1pw+fbpav+r4PV6n9opWq8dKVOqTWxXXmngGBQXRoUMHAFq2bE/Pnv0RQtChQ4caP9WvWbPG4rTJyckEBASQlJREy5YtmTNnDsuWLauQbvXq1QwePJgWLVpUmZder69y/RKSkpI4e/YsgwYNMutVVV6lj6kleU2cOJGJEyeyfv16lixZQlRUFIMGDSI2NpahQ4fi7+9Pjx49yrRVCAwMrJW67b/H69Re0QK/lVQ2oIa9ooprTTxLqtRduya4cMGZvLwGPPWU4Qm0pNjBycmJ0n1RFRQUVJrX+PHjK+01c+rUqYwZM6bMPDc3N7KysvD09CQ7O5tz584xYsQIwNCh4NixY/nqq6/45ZdfiImJYfXq1eTk5KDT6fD09GTevHmmvMytHxoaCsC///1vhg8fjrOzs9njUVVey5YtMwUpS/MCePTRR5k1a5ZpetasWabpp59+mjZt2piW5efn10of9b/H69Re0QK/lVg69J09oIqrNZ7e3hK9Hho2rLisVatW7NixA4ATJ04QHx9faR41eeJv27Yt3333HY6Ojnh7e5e5YYwYMYIFCxYQGhrKypUrTfPXrVvH8ePHywR9g3vV65ewYcMGXn/99TLrLViwgHvvvZcHH3yw2ryCg4PN5lWauLg4U0DfuXOn6XtxcTGZmZn4+flx5swZzpw5U+a9RVxcHLXRkeLv+Tq1N7QyfisxN6i2vaGKa2nPzz93YcuW6p9Mvb2he/diHnywYsdZI0aMICMjg4iICNasWVPmKdVawsPDOXny5B3V7Fi0aBHbt2+vNl1CQgKJiYn06dOnzPyzZ8/SqJFlQ1eUHNOq8irt8sknnxAWFkZERARRUVGmYqPCwkIeeOABevfuzYwZM1ixYoWpqCclJQV3d/daKfdW8TpVlXrplvlOsNdumXNzcy2qZWEPqOJa4pmaKvjqKxdSUx148808MjIEXl6S6rrAycmBoiKBj0/dXtOzZ89m4MCBZl/I1iWjRo1iw4YNFqWt63MfFRWFl5cX48aNu+O8VLtO7R1z3TJrT/xWcvOmMv3JKeNa4hkYKPHzk/TpU8R//+vA4sVu/P3v1fePMmeOB3PnuvPzz47s3193pZgzZ84kJSWlzvKvDkuDPtT9uffx8eGJJ56olbxUu05VRivjtxJ7/6VUGlVcS3s+9ZSh6ObaNUFBgaBRI32163t5SXJzBd9+64K/v6R9+2IaN679fW/UqBEDBgyo9Xzrgro+92PHjq21vFS8TlXFoid+IcRMIURX4/feQogEIcQlIURY3erZL4GBgbZWsBhVXCvzbNFC8u67uTz/fH616y9cmMeSJbl06lTM2bNfUlSUWGm6gwcPVqitUxuuYCgGePzxx+nVqxdhYWFV1uNPSEigWbNmpjYDM2fOrJDmySefrLbVr7We9ogqrqp4msPSop7ngcvG728B7wJvAu/XhZQKJCcn21rBYlRxTUpK5vvvnUlIqHhZXrsmmDnTnagoF06ccKSyxpOffebCnDnuPPRQIbdufUFKSt31m27umD777LMcOXKE/fv3c+TIEXbt2lVpuuDgYFNXCu+++26ZZZs3b66VaoOqnHtQx1UVT3NYGvh9pJSZQggvoAvwoZRyNdCumvV+t1gyoLG9oIrr5cs+7NzpzLhxnqSmCtP8NWtcmDjRk+Jiwa+/OnHypCNjx3qyZIkrOTnFTJs2jfDwcF57rRf79n3IV1/9h+PHjzN58mQiIiLIy8tj9+7d9OrVi2HDhrFlyxbA0OipS5cuZWrodOvWzaLy+6qOqYeHB/369QMMbQ06d+5MYmLlvzyqIjs7m6ioqDL16K1FlXMP6riq4mkOSwP/VSFEODAGiJZSFgshvIHiulPT+D1y7Zpg6VJXMjNFhWV3363nv/91pHXrYtLTBVLCU0958MEHbnTuXMx99xXyz3/mceWKA1On5uPgAM8991+uXUtk795DDBhwDG/vvzJ16gi6du3KihUriI6ORgjBjBkzWLduHdu2bTM9sTk4OJS5ERw9epRWrVpVqCpZvnfLiIgIRowYUW2tnszMTHbs2EH//v0rXZ6QkED//v158MEHiYmJMc1ftGgR06ZNU6LmiIaaWBr4XwTWA3OAhcZ5DwI/14WUCmRnZ9tawWLsyfXTT10ZMqSQr7+uOJiFXn+bTZtuM39+Pu3b67l2TXDtmiN+fpJLlxxwc5N4ekpSUgQvvOCBp6fkpZeacepUAq+//jJJSTto3LgB+eVeB1y4cIGgoCDatGmDEILHHnvMtOyRRx7h+++/B2Djxo0VRtqCsj1llny+/vprU+OwyigqKmLixIlMmjSpTCOqEho3bszJkyfZv38/f//733n66afJysri1KlTXL58uUzjrDvBns59dajiqoqnOSyq1SOl3AY0Kzf7O+PnD0njxo1trWAx9uT62GM61q93YeLEit0nNG7cGDc3uHRJMGCAF0VF8P77uWzd6szx447MmuXBmDE6kpIcadFCUlwsaN/eh//7vyO0a7eNmJjlJCZ+h6fnBxXyLulpsjw9e/bk8uXLpKWlsW3bNl544YUKaQ4cOMCcOXPKzNPr9Xh6elYZ/GfMmEGbNm145plnKl3u6uqKq6srAF27dqV169bExcURGxvLiRMn6NKlC0VFRaSlpTFixIgKY/5aij2d++pQxVUVT3NUGfiFEMFSyivG73eZyeNSbUupQGpqKi1btrS1hkXYk2u7dnrmzKm8hk6J586dzjRtquf4cSc6diymZUs9Dz/shZubpE2bYtq0KSYqyhUHB0l6ejoDBrhx+/YooqJaM23aNMBQDlvyZBYSEkJ8fDyXL1+mdevWZerBCyEYPnw4r732Gm3btsXPz6+CV8kTf2muXr1a5TF98803ycrKMttPfVpaGr6+vjg6OnLlyhUuXbpEcHAwoaGhjB8/HjAUBY0ZM8bqoA/2de6rQxVXVTzNYe6J/xTgZfx+EZBA+ccmCajfcYUVVPUEaY+o4lriOWVKPlu3OnP33cW4u4O7u2TXriwKCgRLl7qxfr0zjRsXs3ixGzdupLFjx2RyciQNG0rmzp0LwBNPPMGsWbNwc3Njx44dvPfee4wZMwY/Pz969+7NuXPnTNt95JFHiIyMNNtLZlWu5bl+/TpLliwhJCTEVNd/4sSJPPXUU2zfvp3Y2FheffVVDh06xFtvvYWTkxOOjo4sWbIEX19fK49czT3tEVVcVfE0h9Zlg5Wo0mwb1HEt77lwoRtz55b9ddCvnxfp6YLsbIGHh8Td3VAcJCUkJTmQkODAzJn51XbvUNuu9ooqnqCOqyqed9xlgxDCv4r5d97rlaKkpqbaWsFiVHEt79mypZ49e5xYvNiNmzcFhYXQq1cRwcF6pk7No3FjSdeuRbzzjhs//uhMbKwjrq6S69frvicSVY+pPaOKqyqe5rD0P+S0EGJY6RlCiGeAI7WvpAalh6mzd1RxLe8ZEKBn6lQP7r23kI0bXbh2zYG8PMHUqfnMnq3j66+zeeaZAvR6KCiAfv0KadtWT1BQ9d071LarvaKKJ6jjqoqnOSwN/OOBVUKIKCHE3UKI7cAUwPywQL9jiovVacKgimtxcTHFxfDhh67GmjxODBlSyKRJDWjXroh161x44YV8Ll1yRK+H48ed+PFHZ1asyOHuu4u5cMGJ8+fr55WTSsdUFVRxVcXTHBYFfinlduBPQF/gNyAd6CGlPFmHbnZNTk6OrRUsRhXXnJwcYmMd2bXLmRUrXJk+PZ/cXMGaNTksX+7GrVuC55/3ICXFgQUL3PHwkLz0Uj4tW0pCQ4tJThY0bVr3T/slriqgiieo46qKpzksegUmhGgA/BPwAd7D8Avgr8BKM6v9rlFpwGVVXJs0aUJOjjS+tJV4e8O8eXksXerG++/n0KgR6PUghOEjJezY4UTDhpLevYvp1i2vXl1VQBVPUMdVFU9zWFrUcwJwBjpLKV/AUMTznBBia52Z2Tm1Mbh0faGKa1JSEh066Fm4MI8PPsgFYPNmF4YOLeSLL9yIi3Pg/HkHSmrTHT7sxCefuPL66x4cO+ZIRobg+PH6KepR6Ziqgiquqniaw9LA/4qUcpyUMhNASnkc6IGh2OcPiSUDVtsLqriWeIaE6Cl5f9avXxE//uhM376FfPSRG1984UpJi/l77ilGSoGXl55GjfQsWeLGvn1OxMbWffBX7ZiqgCquqniaw9IuG76tZF4+ULET8T8IPj4+tlawGFVcK/P805+K+dOfDEU4MTF6cnOhpAp1w4aSjRv/129Ku3bFnD/vSKtWdV/Or/IxtVdUcVXF0xwWN3MRQjwE9AcCKNWCV0r5VB142T1paWm10l96faCKa3We1Q3GMm5cxQHX64rfyzG1J1RxVcXTHJY24JoHrDCmH42hVs8Q4Fbdqdk3Kt31VXFVxRPUcVXFE9RxVcXTHDWpxz9YSvk8oDP+HQEE12RjQghHIUSsEGKLcdpPCLFLCHHB+Lf2OyupI3S6+nu6vFNUcVXFE9RxVcUT1HFVxdMclgb+hlLK08bvOiGEs5TyZwxFPzXhb8C5UtOzgT1SyhBgj3FaCfLy6q/q4J2iiqsqnqCOqyqeoI6rKp7msDTwxwkhOhq/nwaeEUKMAzIs3ZAQogUwHFhVavbDwOfG758DFUfBsFNUqsurimtpzytXHHj9dXeOHLHPzl9VPKb2jiquqniaw9LA/xpQ0lHbbGA68A41q9XzPvASULrKRWMp5Q0A499Gla1oj6hUl1cV19Kee/c6cf/9hezebZ9V51Q8pvaOKq6qeJqjJiNwlXz/Gbi7JhsRQjwIpEgpfxVCDKjJuikpKUyYMAEnJyeKi4t59NFHmTZtGklJSXh6euLo6EhWVhaBgYHcvHkTKSWBgYEkJyebBkXOzs6mcePGpKamIoTAz8+P1NRUvL29KS4uJicnhyZNmpCUlISzszM+Pj6kpaXh4+ODTqcjLy/PtNzFxQUvLy+ys7PJysoiLy+P/Px803I3Nzfc3d3JyMjA39+f27dvo9PpTMvd3d1xcXEhMzOTgIAAMjMzKSwsNC2vi33Ky8sjJyen2n1KT0/H19fXpvuUn59PamoqkZGC2NjG/L//d4OMDC+rz1Nd7ZNOpyMjI8Mm115N9ik7O5uCggKbXXs12Se9Xk9ycrL2/1RL+2Q2JtdHf/xCiLeAcUAR4AZ4AxsxNAIbIKW8IYRoCuyTUrYrva699sefnZ1t0QG2B1RxVcUT1HFVxRPUcVXF8477479TpJSvSClbSCmDgTHAXinln4H/AH8xJvsLsKk+fGqD9PR0WytYjCquqniCOq6qeII6rqp4mqNeAr8Z3gYGCyEuAION00pQF8Pk1RWquKriCeq4quIJ6riq4mkOs4FfCFHrw8lLKfdJKR80fk+XUkZKKUOMf2/W9vbqCpWqdKniqoonqOOqiieo46qKpzmqe+I/Y6y2qVGO/Hzz3QfYE6q4quIJ6riq4gnquKriaY7qAv8o4DUhxFYhRPP6EFIFleryquKqiieo46qKJ6jjqoqnOcwGncCpLQAAIABJREFUfinlfqAzcBI4LoR4VggxqPSnXiztEJXq8qriqoonqOOqiieo46qKpzmqrccvpSwQQiwEOgBvAWmlFwN31ZGbXePm5mZrBYtRxVUVT1DHVRVPUMdVFU9zVBv4hRCRGIZY/BVoI6VMqXMrBXB3d7e1gsWo4qqKJ6jjqoonqOOqiqc5qqvVsxr4EnhZSvmYFvT/R0aGxd0U2RxVXFXxBHVcVfEEdVxV8TRHdU/8rkAnlapZ1hf+/v7VJ7ITVHFVxRPUcVXFE9RxVcXTHNW93P2zFvQr5/bt27ZWsBhVXFXxBHVcVfEEdVxV8TSHrVvuKotKgzGo4qqKJ6jjqoonqOOqiqc5tMBvJSrV5VXFVRVPUMdVFU9Qx1UVT3Nogd9KVKrLq4qrKp6gjqsqnqCOqyqe5tACv5WoVKVLFVdVPEEdV1U8QR1XVTzNoQV+K3FxcbG1gsWo4qqKJ6jjqoonqOOqiqc5tMBvJZmZmbZWsBhVXFXxBHVcVfEEdVxV8TSHFvitJCAgwNYKFqOKqyqeoI6rKp6gjqsqnubQAr+VqHTXV8VVFU9Qx1UVT1DHVRVPc2iB30oKCwttrWAxqriq4gnquKriCeq4quJpDi3wW4lKdXlVcVXFE9RxrYlnYWIhNxbeIO+MbUaY+j0eU3tFC/xWolJdXlVcVfEEdVxr4nlr8y0aPtKQW9/fqkMj2Lp1K4sXLwbg2rVrPPTQQ/Tv35+IiAh27doFwKlTp7j//vsJCwujb9++bNy4sdK8Pv30U/r06UNERATDhg3j/PnzpmXz588nPDyc8PDwKtcvzbJly+jduzd9+/Zl5MiRXL161bTs66+/pnv37nTv3p2VK1fW2AUgKyuLjh078tJLL5nmTZgwgbi4uGrdahshpaz3jdaEmJgY2b59e1trVCAtLU2ZlzyquKriCfbvWvJ/nZ6ebrFn8e1iMr7JwOdBH5ybOteZ25AhQ1i3bh3+/v7MmDGDzp07M378eA4fPszkyZM5ceIEFy9eRAhBmzZtuHHjBoMGDeLw4cP4+PiUySsrKwtvb28Atm/fzurVq1m/fj07d+5k+fLlfPfddxQUFDBixAj+/e9/m9JWxoEDB+jWrRseHh6sWbOGgwcPsmbNGjIyMhg0aBB79+5FCEFERATR0dE0bNjQIpcSZs+eTXp6Or6+vqYb308//cS3337LBx98UCvHtjTHjh37NTIysntly7QnfitxdHS0tYLFqOKqiifYt2txVjHX/naNxNcSK/VMSEigV69eTJ8+nfDwcCZNmsS+ffsYPno4Q1cM5WTiSQDefvttPvzwQ9N64eHhJCQk3JHbxYsXcXV1NfVwKYQwdXqWk5NjKka5++67adOmDQBNmzYlICCAtLS0CvmVDuS5ubkIIQA4f/48ffr0wcnJCU9PTzp27MiePXvMuvXr1w8PDw8AunfvTmJiIgB79+5lwIAB+Pr60rBhQ/r27VtpXlW5ABw/fpzU1FQGDhxYZp2wsDD2799PUVGRWbfaptqBWDQqJysrC19fX1trWIQqrqp4gn27Sr0EAbJQVul56dIlPv30U9q3b09kZCTr169n+/btbN++nffee4+1a9davL3x48dz8eLFCvOnTp3KmDFjysw7cuQInTt3Nk2//PLLjBo1ipUrV5Kdnc2mTZsq5PPrr7+i0+lo3bp1pdtftWoVUVFR6HQ60/qdOnVi8eLFTJ06lby8PA4ePEi7du0s3qe1a9dy3333AZCYmEjz5v8bctzX19d0U7DERa/XM3fuXJYvX050dHSZ9A4ODrRu3ZrTp0/TtWtXi/3uFC3wW0lgYKCtFSxGFVdVPMG+XZ0aOtH8H80RToI8XeUvaoOCgujQoQMA7du3p3///ggh6NChQ42f6tesWWNx2uTk5DJFTxs2bOCJJ57g2Wef5cCBA0yZMoVDhw7h4GAojEhKSuKZZ55h2bJlpnnlmThxIhMnTmT9+vUsWbKEqKgoBg0aRGxsLEOHDsXf358ePXrg5GRZuPv222+JjY1ly5YtwP+KzUrw9PQs8zRfncvq1asZPHgwLVq0qHSdwMDAen9npAV+K7l586bpZ6G9o4qrKp5g/64OboYgeTOpck8XRxdkoUQ4Cxz+f3tnGh5FlS7g9/SW7gTIzhIQCAgDjKIIAmEJKKhEBNlUELdhrssFH0RQAXVklHFf0Dug11FxHBUdkc2Ni4KyiePGvsiwhCCQhGwkJPTe5/7oTpOQTqfThnSfSb3P00/SVaeq36qu/rr61Knv0+n8aQh0Op2/28FgMFQLena7PeBr1eeM32w2U1ZW5n/+3nvvsXTpUgDS09Ox2+0UFRWRmppKWVkZEydO5OGHH+byyy+vc5vHjRvHrFmz/M9nzZrlf37nnXf6u46CsX79el588UU+++wzYmJiAGjbti2bN2/2tzly5AjXXHNNyC4//vgj3333HW+99RYVFRU4HA7i4uKYN28eADabrdHr+GqBP0yi/aJ4VVRxVcUT1HEN5Ok44sBV6CJ3fi5pT6TVumz79u1Zs2YNADt27CAnJydgu/qc8Xft2tUf6AHatWvHxo0bufnmmzl06BB2u52UlBQcDge33XYbN910E2PGjKl1fYcOHfIH9C+//NL/v9vtprS0lKSkJPbs2cOePXv8/etPPPEEl112Gdddd121de3cuZOZM2eydOnSar/orrzySubPn8+pU97RTlu2bOHZZ58N2aXqKKAlS5awfft2f9CvXK6xB7BogT9Movmn/rmo4qqKJ6jjGshTWARSSvQJwS9Qjxo1ig8//JDMzEx69eoV0hlzXQwYMIA//elPSCkRQjB//nxmzJjBa6+9hpSShQsXIoRg5cqVbNmyheLiYj744APAO9zy4osv5qmnnqJXr15kZWXxxhtvsGHDBoxGIwkJCSxatAjw3mR17bXXAtC8eXNef/11f1fP3r17GTFiRA23efPmUVFRwR/+8AfA+6W0ZMkSEhMTeeCBBxg2bBgADzzwgP+6SSguwTh58iQWi6XR7w3QhnOGSU5ODh06dIi0Rkio4qqKJ6jjWtXTnm2n6B9FJN2UhLmbGXeZG+mWVGyp4MzPZ0idnooh4fyfC86ZM4cRI0YwdOjQWl3PJ+PHj2fZsmVhL9+Qnq+++irNmzfn1ltvbZD1VUUbznkeaNasWaQVQkYVV1U8QR3Xqp6nVp0i8YZETq30dlmc+NMJcv+cS/mmcmK6xGDd2Th37M6cOROrteZrNdY+/S1BHxrWMz4+nkmTJjXY+kJFC/waGk2E5DuSOf3NaRLGJyCdEnTgLHKS/F/J6FvoiRsU1ygeLVu2JCsrq1FeK9qZPHlyyKONGhIt8IdJeXl5pBVCRhVXVTyhcV3t2XYqfqwIa9ny8nLsh+zkPZNHwV8LEDpB+fpyjk47yul1p3EXufl1+q/Y9tk4Pus4+S/lA96LkLm5uQHXuXnz5hqjdRqCyn06YcIEBg8eTEZGBjNnzsTtdtdoe/ToUdLS0sjMzCQzM5OZM2f65y1fvpxBgwaRkZFR7SJqQ3uqjHZxN0xatWoVaYWQUcVVFU9oXNeC1wqIb1+IJ609urZJdbaXLknZ2jIsF1twLXBx+qLTpE5LRd/87MXc2N6xeKwebL/YEEaB9RcrxhQjHqsH8Oam6d69O61btyb3T7lIlyTtyTSEPvD49Yagcp8uXryYFi1aIKXk9ttvZ+XKlYwfP75G+44dO9a4Iaq4uJh58+bxzTffkJKSwtSpU9mwYQNDhgxpcE+V0c74w6SgoCDSCiGjiqsqntC4ruauZspLWkFqQt2NgZKlJUi35MCIA7iGu/BUeKoFfYDYy2JJfzedrv/qyoIOC7jh6xu4Yd0NrEhcwapVq9i+fTt33303Q4YMoaK4gk0HN9Gvfz+ysrL8NzZ5PB4uueSSavnpe/fuzcmTJ8Pazsp9Wpn6wOVy4XQ6a71ZKhBHjhzhwgsv9N8kNmTIED799NOwfOryVBntjD9M6nMwRhpVXFXxhMZ1Tfmv+iWDi8uI49h9x7wXcvNO4S5zU/hWISl/rLmeXbt2cdJ2ku93f0/uE7k0m9yM+Ph4/rbob8x/cj6XXX4ZFScrGDt0LKs+XUWnTp2YMmUK4L3Zq/KLYPLkyfz000+0b9+eli1bVnuNTZs28cgjj9R4bYvF4r9PAKrv0/Hjx7N161aGDx/O9ddfH3A7jx49ypAhQ2jevDmPPPIIGRkZdOrUiX//+9/+rqDPP/+8wfPnq3Sc1oYW+MMkKanun9zRgiquqnhCdLsaUg1YLrGQ9lga8eXxnHrqFPZsOx6rG52l+pl/y/KW5GTnMPfxuWTEZTD8wHDKHeW48l0Uv1sMl8Ph/MN0SO/gH8d/44038s477wAwduxYnn/+eSZPnszy5csD3mw1ePDgGl0ygai6T5ctW4bNZuOuu+5i48aNNZKbtWrVip07d5KUlMT27du55ZZb2LJlCwkJCbz44otMmTIFnU5H3759a73xLFyi+b0PlUbp6hFCXCCE+EYIsU8IsUcIcZ9vepIQ4ishxAHf3+jMehUAlX7uqeKqiidEr6v0SI5MPoJtvw3p9FBYVIgxzYh1h5Vjs46fbeeU5L+cj2GbgY+u+4iBAwfy4Z4Pue/++zCkGJAuiT7l7JdEbWe5ffv2JTs7m8LCQr744gtGjRpVo82mTZv8F2GrPs5Ne3DuPjWbzWRlZbF69eoa64yJifEH4EsvvZT09HR/XvsRI0awdu1avvzySy688EI6deoU4t4LjWh97+tDY53xu4BZUsqtQojmwM9CiK+AO4B1UspnhBBzgDnA7EZy+k0Ey+sdbajiqoonRLeru9SNLkFH7l/y0PXT4cx14i514453YztgI6ZzDIV/K8Rd5ib7/WySrkpi9OjRpKenc9eEuzB3NZPcMxnDFd7w0KVLF3JycsjOziY9Pb3aOHghBCNHjuTRRx+la9euAc+GQz3jb9GiBeXl5ZSXl9O6dWtcLhdfffUVGRkZNdoWFhaSmJiIXq/nyJEjHD58mI4dOwLewJyamsqpU6dYvHhxvVJKhEI0v/eh0iiBX0qZC+T6/j8thNgHtAWuB4b6mr0DrEeRwB9oiFm0ooqrKp4Qva5CJ2hxdQs8FR5iusSQ/30+MbtNGFL0WHpaKFtThuOog9hesZR/W461n5U7v74TXab3x//9Q+5HSsmkSZOYNWsWZrOZNWvWsGDBAiZOnEhSUhL9+/dn3759/tccO3Ysw4YNCylFQTDcbjdOp5PJkydjt9txu91kZmb6UyisXr2abdu28fDDD7NlyxaefvppDAYDer2eF1980Z9GYe7cuezevRuABx98kAsvvPA3eQXyVJ1GT9kghOgIbAQuAo5KKROqzCuRUlbr7lm1apWcO3cuBoMBt9vNuHHjmDZtGnl5ecTFxaHX6ykrKyM1NZXi4mKklKSmppKfn++/w668vJxWrVpRUFCAEIKkpCQKCgpo0aIFbrfbXwAiLy8Po9FIfHw8hYWFxMfH43A4sFqt/vkmk4nmzZv7zzCsVis2m80/32w2Y7FYKCkpITk5mdOnT+NwOPzzLRYLJpOJ0tJSUlJSKC0txel0+uefj206ffo06enpdW5TZXWgSG3T8ePH6dq1a4O+T+drm86cOUPr1q0b5diL2RiD0/M1pm6difv97+vcpl8/+5W0fmmUZJdwSneKuDficFzpoNX1rTjz2RnKRTmtu7fG0d6BNddKy8SWFOuLiTXHUraiDEOWodE/T0VFRdjtduLj47XPUwNtU05OTq0pGxo18AshmgEbgCellMuFEKfqCvzRmqvHbrf707ZGO6q4quIJjeua/3I+ssJNs4Fmmg2t+8Ji/sv5JN+RjCHBgN1uRxQIdHE6DIneNMvSKtHFes/wj806hjB6k7ZJm0TECNo9Fzhv/PlGlfdfFc+oyNUjhDACy4D3pZSVlY/zhRBtfPPbAOENAI4AqhTbBnVcG9NTFBVhWrwYquSGrw+N6dryvpak3NOyRtB3nXJR8a+KGqmXE8cnUvJRid/T1M7kDfpub7vKoA94a+sKEEaB+7QbU1vTed6a2tGO08ajsUb1COAtYJ+U8qUqsz4Bbvf9fztQs+5alGI0nr9i1A2NKq6N6Rnzt7/hysggJswLf43pKoTAkFzzctzJBSex7yvF8e2v1aabLjDhPO5EuqXf037IzvEHj5P/XD7SJSn6exG2fTZaPdCKtKfSMKQYcJe4SZ6S3CjbFAjtOG08GuuMfyBwK3ClEGK773Et8AxwlRDiAHCV77kSxMfHR1ohZFRxbUxPx9ixxHz0Ec4gRT6CUW9XKTF+8gm6X34J6/UCYelpwXpQouvcpsa8pElJ5D+f7x+BIj3epGy4JeWbyzGkGCh+v9g7zyY58/MZpFuii43czUnacdp4NErgl1JullIKKWVPKeWlvscXUsoiKeUwKWUX39/ixvBpCAoLCyOtEDKquDamp6dbN6zz5uHxDQGsL/V11e/a5e1eeucfYb1eIBLHJ9L2ybbe7ppzMHcz03xoc3JW5iCdEnMXM9sHbecfhneJvTwW6x4rGywb6N+3P/0u6cejvzxK2p/TEHodx44dY9y4cfTr14/+/fsHrMG7aNEi+vfvz6BBgxgzZgy//ur91XHumP02bdrw+eefB92OY8eOMXr0aIYPH86gQYP46quvArZbtmwZAwcOZNCgQUyYMIGioqKgLoWFhUyYMKFe+zQUVPk8BUMrxBImp06dIiEhtNwpkUYVV1U8IQxXmw3zggW4u3XDOXbs+RM7h4I9BbiWudC10HHrqlv54OMPSE5O5sCeA9wx6Q7+J/l/aNW3JfZeDjpP9N6ZO2rUKGbOnMkVV1xBeXk5Op2uRt3eTZs20bt3b2JjY1m8eDGbN2+uMV6+pKSE3r17s3v37qD1iWfMmEHPnj0ZN24ceXl53HTTTezYsaNaG5fLRY8ePfjuu+9ITk5m3rx5WCwW5syZE9Rl2rRp3HrrrfTv378hdiegznEaFRd3/9NwOByRVggZVVxV8YQwXM1m/j15Mpc98wzTp09nwIAB3HXXXaxfv54RI0bQp08ffv75ZwCeeeYZ/vrXv/oXHTBgQMCz7lCQqZI2j7Xh1KBT6E7rcL7hJO/pPF5/8HWmTJ5C5xmdEQYdFwy4AIBffvkFl8vlT5HQrFmzgEF78ODB/ul9+vThxIkTNdqsWrWK4cOH11mUXgjhH9JYVlYWsAyhlBIpJWfOnEFKyenTp/3tgrmMHDmSjz/+OJRdFTIqHae1oeXqCZNAFYSiFVVcVfGE0FxNS5cicnOx33sv6LznWIcPH+btt9+mW7duDBs2jI8//pjVq1ezevVqFixYwHvvvReyw5QpUzh48GCN6VOnTvXny6/0/Gn/T/QZ0YfWc7zBMm9PHnFn4pj42kTcbjezr5jN8PbDOXToEPHx8dxy0y3kHMnhimuuYN68eej1tdfnfe+99xg+fHiN6StWrGDq1Kl1bsfs2bMZP348r732Gna7nRUrVtRoYzQaeeGFFxg4cCBxcXF06tSJ559/vk6XSy+9lCeffLJOh/qg0nFaG1rgD5PGLo78W1DFVRVPOMe1rAwC3Mav37oVT2oqIjcX2bYtAB06dKBHjx4AdOvWjSFDhiCEoEePHvU+qw8lFUGr+Fbkzs8l+2A2Kb3OZue0F9vZn7+fT7/4lBMnTnDttdeyZcsWXC4X3333Hf8c+U/aXdyOOdvmsGTJklprwn700Uds27bNn6q5kry8PPbu3cuVV15Zp+OyZcuYNGkSd955Jzt27OCee+5hy5Yt6HRnOyScTidvv/02GzZsoGPHjsyePZsFCxbwwAMPBHVJTU1t8OGXKh2ntaF19YSJSmN5VXFVxRPOupoWL8a8cCHGlStrtHF364bhX//CVOUM1mQ6O05ep9P5n+t0OlwuFwAGg6Ha2Hy73R7QYcqUKQGTn3344Yf+Nsf2HQMjGIp12Gw2//Sk0iSGthuK57iHDh060KVLFw4dOkRaWho9e/bkwl4XorPpGDl6JDt37gz4+uvXr+fFF19kyZIlNW5oWrlyJSNHjgxp6ON7773HmDFjyMvLo2/fvtjtdv+F20p27doFQHp6OkIIxowZww8//FCni81mw2w21+lQH1Q6TmtDO+MPk6of4GhHFVdVPKGKa1wcorwcGVezXq27Rw/0+/ZBPQZQiNJSOp45w2pfLpwdO3ZUTyvsciHy85Ft24Z0xu/8yYl9TwU9R1/Eqp/Ojq7JGp7FqvWruLvd3RQVFXHw4EE6duxIfHw8p06dgrGQlpLG5mmb6dWrV4317ty5k5kzZ7J06VJSU1NrzF+2bBmPPfZYtWlPPPEEl112Gdddd1216e3atfOnXt6/fz92u91fSKWSNm3asH//fgoLC0lJSeGbb76ha9eudbocOnSI7t2717mf6oNKx2ltaIE/TJo3bx5phZBRxVUVTzjr6hgzBtOZM3h8Qagq7ssvx3rJJVCPQBHzyiuMad+RD9auIzMzk169evnz4AOYn3sOj86Ec+hwZP9La11P+aZyXIUuLB4LhjYGMsdfyfzFTyOlRAjB9X++nh8f/ZGMQRno9Xoef/xxf2bNJ554gjFjxiClpEdKD4Y3G451t5UFnyygV69eZGVlMW/ePCoqKvwJ1Nq1a8eSJUsAb4GUEydOMHDgwGpOe/fuZcSIETVc58+fz4wZM1i0aBF6vZ6FCxf600BnZmayceNG2rRpw0MPPeT/FXHBBRf4k8IFc9m0aRNXX311yPs/FFQ6TmtDG84ZJjk5OXTo0CHSGiGhiqsqnnDW1fjxx8g2bTB++SXWxx+v1ka/axeGDRsQpaXei7tCgNuNTEjAecUVeHx9/VUxfvYZxjVrkElJ2GbORPpuFhLHjmF56SVEfj5He8whbkIXzL+rvQvj+NzjNGt/mpIhOjr28H4pzZkzhxEjRjB06NCQtzP/+XwMbQzoTDoSb/xt5TLGjx9fLaXzuZyP93/kyJG8//77DTr8UpXjNNhwTu2MP0wqU8CqgCquqnjCWVfdyZOYXn0V2513+ucZ1q7FuGED7p49cdx6qz94VyJKSjCuWYP+gw9wXnUVrsxM/zzndddh+PlnbBMmoF/xCa47vBdVhdWKNJkovPwukmdcXKdf7GWxOApMJLY14TnjIXd+Lrc1v4Uc6691LluV1Omp2PbasFxqqddygQgW9KHh3//CwkKmTp3a4GPuVTpOa0O7uBsmKg3pUsVVFU846yqKi5FJSWA2Q0UF5vnzETYb1vnzcdxwQ42gDyATE3FMnIj1iScQBQWYn3wSqlx4tU+cSMxnn+EePdI/zdOlC7Y5c4i9d2CN9QUi8YZEUqemYrPZcJe6QQdJNgv9SvuR90weHocnpPXoYnTE9optlDqzDf3+p6SkMHLkyLob1hOVjtPa0AJ/mFQdIRHtqOKqiiecdZUJCVjvuw/9tm1Y/vxnnNdcg6iogFCKdQiBc/x47HfcgWXePPAFFM/vfodt9mzvF0oVZEICpnffxfz00+AbARSKp7GNkZS7Uki5vwPW7VYMzd24fi2v3wY3Aqq8/6p4BkML/GGi0lheVVxV8YSzro6xYzGuX4+usBDbQw9h+uQTpMmEIYRSg5XItm2x3XcflueeY8mSJeTm5gZst3nzZm54+WVslva4Cyrq5RnTIQZDkoGk25IwdW2OqXPN+w5GjRpF3759/cNCa6stu2DBAnr37k3fvn1Zt25dvZevyzXaUcUzGFrgDxOVxvKq4qqKJ5x1lW3b4u7RA1dGBrpjx3AOGYJ+715cvXvXa30yLQ1nRgYfLlwYdD94LrgAQ9909G1CyxB57rosPSy0uKr2mrGvv/46GzduZOPGjQGHaf7yyy8sX76cLVu2sHTpUh588MFqpQjrWr4+rtGKKp7B0AJ/mDT0TSHnE1VclfB0uxHFxWddXS4MO3ei37ED00cf4Ro4ENsjjwS8k7f6atxMmzaNAQMGMHDgQF599VWWWa1sO3iQu2+/nczMTKxWK2vXrqVfv35kZWXx2WefIZs1w9G/P5dccgmlpaX+9fXu3ZuTJ2vWMWrofbp69WrGjRtHTEwMHTp0ID093Z9j6LeixPuPOp7B0Eb1hInF8ttHOTQWqrhGzNPhwPzyy7gvvBDnuHFBm5qfeQZhsxE3ahSGffsQpaU4srIwfPstxo0bMXz7La6rrgJAl51NzDvv4JgwAfdFF1Vbz65duzhx4gRbtmwBoLS0lPj4eN6eOZOnu3Sh+yuvYBOCGTNmsGrVKjp16sSUKVO869Xp/F8EkydP5qeffqJ9+/a0bNmy2mts2rSJOXPm1MizY7FYWLNmTcDtu/fee9Hr9YwaNYoHHnigxkXd3Nxc+vQ5O0IwLS2tWtdUXcsHQztOGw/tjD9MSkpKIq0QMqq4RspTd+gQMjERw48/1tlWJiSA1UrZ3r3gcBDz5psYduzA07Ur9ttvx7h+vb+t8dNPcYwbh/GTT2qsp2PHjuTk5DB79mzWrl3rvylIxsejy8+H1FQOHDhAhw4d6Ny5M0IIbrzxRv/yY8eO9SczW758OWMCFJQZPHgwS5Ys8Xe9VD5qC/qvv/463377LZ9//rk3X88//1lz+wPc91MZ3ENZPhjacdp4aIE/TJKTI1eirr6o4hopT0+3bnhatsQxeXKdbe3TpmF94QUSfvc7jBs3IqxWnBkZ3nw9n36Kvco67FOmoN+xA/s999RYT0JCAhs3bmTgwIG89fLLzLjlFq9L27a4evUCX76Z2s6Y+/btS3Z2NoWFhXzxxReMGjWqRptNmzYxefLkGrl8rrnmmoDrTEtLA7x3pk6YMIGtW7cGbHP8+HH/8xMnTvgvdoayfDC047Tx0AJ/mJw+fTrSCiGjimvEPIXAef31NbpjglHasSPWv/wFV9++xHzwAe7OnSEpqfrduHFx4PGgDxAAi4qK8Hh3OImCAAAQSElEQVQ8XN+rF385eZKd330H5eU0a9aMct9QzS5dupCTk0N2djZQ/QYoIQQjR47k0UcfpWvXrv50C1UZPHgwS5cuDemM3+Vy+ROjOZ1O1qxZEzDHzYgRI1i+fDl2u52cnBwOHz5M7969Q14+GNpx2nhoffxholIxBlVcVfGEs66yeXNsc+ei27sXzsnhojt2DGG1YvzmG1zn5KvPzc3l3nvvRRYWIsrLebJFC4iLY9KkSUy//35iMjNZs2YNCxYsYOLEiSQlJdG/f3/2+ZK3gbe7Z9iwYf6cNcE868JutzNhwgScTidut5shQ4Zw2223Ad4Lutu2bePhhx+me/fujBkzhoyMDAwGA8899xx6vZ6Kiopalw8VVd5/VTyDoeXqCRO73V4jFW20ooqrKp643XjWrUPfty8xixZ5R/EEIOaZZzCuW4dz+HDss2cjiovBZkP6ukQAsNmIu/lm3N27Y581C5mUhPmpp7A9/HCDqCqzT1HHVRVPrfTieUClsbyquKriafr4Y07odMQsXFhrG93Roxi++w59djbG778HlwvL/PmYX3rJW7ilErMZ59ixyFatkJU5YOo6GZMS85NPYp4/v862quxTUMdVFc9gaF09YaLSkC5VXFXxdPXvT7PPP8c1YAD6/fsReXnIc+7m9KSk4LnoItxnzuBJTwch8LRsiSgr8+b1qYKjSnUrXU4Onqq/CM7B/PzzYLWiKy5G6vRIhwMR5OxTlX0K6riq4hkMLfCHiUrFGFRxVcXT06ED3HwzroQE3D17YlqyBPv06dUbxcbi9OWedw4ZAno9trlz61y3aelSbHfdFXimlOiOH8djNGGdOhVhNgcN+qDOPgV1XFXxDIbW1RMmVe+ajHZUcVXFE866ypQUREUFIkBeGt3Bg7i7dEEfoCC6fts279l7lREi4sQJb87+2u76FQLbjBk47roT2bUrnvbtQ/ZUAVVcVfEMhhb4w+Tc0nDRjCquqnhCdVfbtGmYX3mlRn+74447cLdvjy43F8PmzdXmmZYvxzl4MMavv/ZO8Hgw//Wv2P77v4O+rqdjRzxduoTlGe2o4qqKZzC0wB8mKn3rq+KqiidUdzUtXw5lZZhfeKFa8De9/z6WhQvxtGiB8f/+z5tz32bD9Pe/4xg2DMMPP+C8+mrvxdpnn8Vx880QG1vzxZxOYl55pV4ZPwN5RjuquKriGQwt8IeJ0+mMtELIqOKqiidUd9Xv24ds0wbX73+PZd48b6lFAJcL18UXIwoKECUlWB57DNOKFbi7dcO0fj326dMRViuWRx/FmZWF++LAlbX0W7fi6dABYy2pFkL1jHZUcVXFMxjaxd0wUSkntyquqnhCdVf7HXegKy7GNXAg7n79ML/0Eu7OnXHccIO31q7LRey8ecgWLXAOH07M//4vjiuuwPT22+iPHsX24IPeHEC14O7TB8Obb+IIkI+nPp7RjiquqngGQ7uBK0xUKbgM6riq4gl1u+r+/W9MK1Z4K3FZLHh8mTN1+fneLh+jEcf48Xg6dfIuYLej37kTd+/e3i+LRvKMJlRxVcVTK7Z+HoiLi4u0Qsio4qqKJ9Tuqt+9G9M//4lj3Dhss2d7J545g66w0DuWPyUFAowDNy9ciKd9e/R79+K4/faA6za9+y76PXu8aZ77BPw8h+wZjajiqopnMLQ+/jA5N8d5NKOKqyqeULur4euvcY4YgbFKSUJiY/G0b4/nggsCBn0AV8+e6Ldvx33ppbW+piguRsbHI06d+s2e0Ygqrqp4BkM74w+TsrIyEitvsY9yVHFVxRNqd7XfeSfGdevqHJZ5Lq6rrvIXcKkN+/TpiMJCZD1KGv4n7NNoQxXPYGiBP0zqW080kqjiqoonBHG1WHBed935eVEh6hX04T9kn0YZqngGI+JdPUKIEUKI/UKIg0KIOZH2CZXi4uJIK4SMKq7R7Gn48kvvnbZ2OxDdrlVRxRPUcVXFMxgRDfxCCD2wCMgCegCThBA9gi8VHUT7aKiqqOIazZ7GDRu8F1937QJCcI2SbYnmfXouqriq4hmMSJ/x9wUOSikPSykdwIfA9RF2CgmVfu6p4hrNno6JE8HpxH3ZZUBwV93Ro8Tefz/mZ59tLL1aieZ9ei6quKriGYxIB/62wK9Vnh/zTYt68vPzI60QMqq4RrOn++KLcdxyi3+MfVBXlwtpNEIUVGqK5n16Lqq4quIZjEhf3A1USbra76iTJ0/yxz/+EYPBgNvtZty4cUybNo28vDzi4uLQ6/WUlZWRmppKcXExUkpSU1PJz8+nWbNmAJSXl9OqVSsKCgoQQpCUlERBQQEtWrTA7XZTUVFB69atycvLw2g0Eh8fT2FhIfHx8TgcDqxWq3++yWSiefPmWK1WysrKsFqt2Gw2/3yz2YzFYqGkpITk5GROnz6Nw+Hwz7dYLJhMJkpLS0lJSaG0tBSn0+mffz62yeFwUFFRUec2FRUVkZiYGLFtqnzdhnyfztc2eTweSkpKAm+TXk+zmTPBaKQ8J6fBj736bJPVasVut0fs2KvPNul0OvLz87XPUwNtU9DAG8n+KiFEBvBnKeU1vudzAaSUT1e2idY7d4uKikhOTo60Rkio4qqKJ6jjqoonqOOqimc0l178EegihEgXQpiAicAnEXYKifLy8kgrhIwqrqp4QiO62mzot28P+2Kxtk8bHlU8gxHRrh4ppUsIcS+wBtADi6WUeyLpFCqtWrWKtELIqOKqiic0nqt5wQI8rVujO3oU5+jR9V5e26cNjyqewYj0GT9Syi+klF2llJ2llE9G2idUCgJUXIpWVHFVxRPOg6uUGD/9FP3u3dUmu7t3R3/gAO6uXcNabZPep+cJVTyDEfHArypCBLouHZ2o4qqKJzS8q27/fnQnTmB6//1q051jxmB96ik8YV7nasr79HyhimcwtMAfJklJSZFWCBlVXFXxhIZ39XTujCgtxRVi1s1Qacr79HyhimcwtMAfJir93FPFVRVPOA+uRiO2hx7COX58g662Se/T84QqnsHQAn+YtGjRItIKIaOKqyqeoI6rKp6gjqsqnsHQAn+YuN3uSCuEjCquqniCOq6qeII6rqp4BkML/GFSUVERaYWQUcVVFU9Qx1UVT1DHVRXPYGiBP0xUKrisiqsqnqCOqyqeoI6rKp7B0AJ/mOTl5UVaIWRUcVXFE9RxVcUT1HFVxTMYWuAPk5UrV0ZaIWRUcVXFE9RxVcUT1HFVxTMYWuAPk+XLl0daIWRUcVXFE9RxVcUT1HFVxTMYWuAPE5fLFWmFkFHFVRVPUMdVFU9Qx1UVz2BENC1zKKxbt64AyIm0x7kUFxenJCUlFUbaIxRUcVXFE9RxVcUT1HFVxRPoMGzYsIDlwqI+8GtoaGhoNCxaV4+GhoZGE0ML/BoaGhpNDC3wh4EQYoQQYr8Q4qAQYk6kfSoRQlwghPhGCLFPCLFHCHGfb3qSEOIrIcQB39/ESLsCCCH0QohtQojPfM+j1TNBCPGxEOIX377NiEZXIcT9vvd9txDiAyGEOVo8hRCLhRAnhRC7q0yr1U0IMdf3+dovhLgmwp7P+977nUKIFUKIhEh7/la0wF9PhBB6YBGQBfQAJgkhekTWyo8LmCWl7A70B6b53OYA66SUXYB1vufRwH3AvirPo9XzFeD/pJTdgEvwOkeVqxCiLTAd6COlvAhvRbuJRI/n34ER50wL6OY7ZicCv/ct86rvcxcpz6+Ai6SUPYF/A3OjwPM3oQX++tMXOCilPCyldAAfAtdH2AkAKWWulHKr7//TeANUW7x+7/iavQOMiYzhWYQQ7YCRwJtVJkejZwsgE3gLQErpkFKeIgpd8ZZStQghDEAscIIo8ZRSbgSKz5lcm9v1wIdSSruUMhs4iPdzFxFPKeWXUsrKMZz/AtpF2vO3ogX++tMW+LXK82O+aVGFEKIj0Av4HmglpcwF75cD0DJyZn5eBh4CPFWmRaNnJ6AAeNvXLfWmECKOKHOVUh4HXgCOArlAqZTyS6LM8xxqc4vmz9gUYLXv/2j2DIoW+OtPoLprUTUmVgjRDFgGzJBSlkXa51yEENcBJ6WUP0faJQQMwGXAa1LKXkAF0dMF5cfXP349kA6kAXFCiFsiaxU2UfkZE0I8grc7tbI+ZlR6hoIW+OvPMeCCKs/b4f1JHRUIIYx4g/77UsrKe8vzhRBtfPPbACcj5edjIDBaCHEEb1fZlUKI94g+T/C+38eklN/7nn+M94sg2lyHA9lSygIppRNYDgwg+jyrUptb1H3GhBC3A9cBk+XZm5+izjNUtMBff34Euggh0oUQJrwXdz6JsBMAwlsF+i1gn5TypSqzPgFu9/1/O7Cqsd2qIqWcK6VsJ6XsiHf/fS2lvIUo8wSQUuYBvwohfuebNAzYS/S5HgX6CyFifcfBMLzXeKLNsyq1uX0CTBRCxAgh0oEuwA8R8AO8o/iA2cBoKeWZKrOiyrNeSCm1Rz0fwLV4r+4fAh6JtE8Vr0F4f2ruBLb7HtcCyXhHTRzw/U2KtGsV56HAZ77/o9ITuBT4ybdfVwKJ0egKPA78AuwG3gViosUT+ADvtQcn3jPlPwZzAx7xfb72A1kR9jyIty+/8jP1v5H2/K0PLWWDhoaGRhND6+rR0NDQaGJogV9DQ0OjiaEFfg0NDY0mhhb4NTQ0NJoYWuDX0NDQaGJogV9D4zwjhGgvhChXJYGXxn8+WuDXaJIIIZoJIY4IIW6uMq25EOKoEGJCQ76WlPKolLKZlNLdkOvV0AgXbRy/RpNFCHE13rwrPaSUBUKI1/AmDhsXYTUNjfOKdsav0WSR3uyVnwP/I4QYCtwITKutvRDiD75CLKeFEIeFEHdXmTdbCPEvX0pkhBD/7SuKYhZCdBRCyCrz7vAtf1oIkS2EmHxeN1RD4xy0M36NJo0vq+VewAg8KKV8O0jbkXhTIhzGm6N/NTBISrlVCKED1gNf4/0V8QNwpZRymy9FdrbvNWLwpgS4XEq535ecLElKuef8bKGGRk20wK/R5BFCrMWbybKNlLK0HsutBL6RUr7ie94R2ArkA/+QUj5dZXrVwH8cbw6YL6SU1obaDg2NUNG6ejSaNL6c9R2BtcCzdbTN8nXnFAshTuFNgJdSOV9KeQT4xre+RYHWIaWsAG4C7gFyhRCfCyG6/fYt0dAIHS3wazRZhBAtgQXAncDdwI1CiMxa2sbgrXPwAt4LwAnAF1QpxiGEuBbIwJtp8vnaXldKuUZKeRXQBm/X0RsNskEaGiGiBX6NpsxCYKWU8hvpLf33EPCGL8ifiwlvN00B4BJCZAFXV84UQqTgrYXwX3hzy4/yfRFUQwjRSggx2le+0Q6UA9owT41GRQv8Gk0SIcQYvPULHqycJqV8E28O9sfObS+9xeunAx8BJcDNVC/A8zdglZTyCyllEd4+/DeFEMnnrEoHzMJbqakYGAJMbaDN0tAICe3iroaGhkYTQzvj19DQ0GhiaIFfQ0NDo4mhBX4NDQ2NJoYW+DU0NDSaGFrg19DQ0GhiaIFfQ0NDo4mhBX4NDQ2NJoYW+DU0NDSaGFrg19DQ0Ghi/D9rQBPnaOsodQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.loadtxt(\"exp_data5.txt\")\n",
    "run_test()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circles are the different Gaussian distributions. The centre of the circle is the Mean (mu) and the radius of the circle is equal to the Standard Deviation (stdv) of the corresponding Gaussian distribution. \n",
    "\n",
    "All the Gaussians are initialized randomly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"result.gif\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "   Displaying a GIF to demonstrate how the Gaussian statistics \n",
    "   iteratively converge to the data points. \n",
    "   Obtained by saving the above plot for iterations from 1 to 30\n",
    "   of the algorithm. \n",
    "'''\n",
    "from IPython.display import HTML\n",
    "HTML('<img src=\"result.gif\">')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that after each iteration of the algorithm, the Gaussian statistics (the circles) adjust more and more to the data points. After some iterations (30 in the above case) the algorithm converges. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Missing(%) vs Mean Absolute error for Experimental dataset 1\n",
    "\n",
    "In this experiment we try to obtain a plot between the percentage of missing data, and the error obtained after recovery (in terms of Mean Absolute Error) for this method. We use the Experimental dataset 1 for this purpose. \n",
    "\n",
    "$$ \\mathrm{Mean Absolute Error} = \\frac{\\sum_{i}^{N} \\mid X_i - Y_i \\mid }{N} $$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEXCAYAAACqIS9uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dd3hcxbm4309lrS6rrCVXuVBsTIIhBEIoBhQIXEhInOQSSCVO901vhPRGcm/uTfsluQkJhDQgXELgxqFDqNdAsIGAsQ0YW64rrYp3VVZaafX9/jgrIYTKarW7c0aa93n2sU6bfWd2/e2cmTkzoqo4HA6HY/aQZ1rA4XA4HLnFBX6Hw+GYZbjA73A4HLMMF/gdDodjluECv8PhcMwyXOB3OByOWYYL/A7rEBEVkXdmKK2rReSuTKTlcNiCC/yOtBCRt4vIbhE5JCK/EZHCEcfyReQREbkwhXTemwzkoZFpJI8FRaQvefyUEYfmAzdkKCufAN6WobRyTvIz0DFeG027OfyLC/yOKSMiNcBvgC8BpwKvBT444pRPAwdU9U8pJpkABoA3jNp/CXBw9MmqGlLV3ql6j4WqRlS1IxNpGeTf8X4MR77GvSMSkcBU9k+GiOSJSH461zrM4AK/Ix2WAxFV/aOqPgXcBBwFICKH49WiPzLFNK8CPjC0ISICvB+4cvSJo5t6ROT9IrJNRHpFpE1E7heRRcljFck7klDy7mGviPxgxLUvaeoZ2haRD4pIk4hEReRmEQmOcvikiOwTkR4RuV1E3pX0WjRW5kTkAyISEZHiUfu/ICL7k8GzUER+kEy3T0QOish1KZRdV/LHcOTr0Kjy+riIXCMiEeCPIrI0uf8dInKLiHQDlyfPf02yDGMi0pG8bt6I9L4uIs+LyIUish2IA6tS8HT4BBf4HenwPFAiIseLSCmwFng8GayvBC5T1dAU0/w18DoRaUhunwEEmaRJR0ReBfwC+C5wJHA68LsRp3wbOA64ADgcuBDYNonLq5Pvfx5wDrAG+M8R77kuuf194BjgWrxa90RcDwSAN43a/y7gD6o6CHwM+Fe82vrhwBuBhydJN1W+BmzCK4svjdj/78A1wCuAn4lIPXAHsA84Ae8u7Gjgz6PSWwB8FHgv3o9+U4Y8HblAVd3Lvab8wgsITwIvAD8E8vEC19+Aery7gBfwgnDZBOm8FxhI/n0L8I3k39cBPwGWAgqcMuIaBd6Z/PvNQASoGCf9m4GrJ3j/q4G7Rm2HgTkj9l0KHByx/RDw+1HpfC/ptWiC97oOuHXE9nHJa1Ynt38M3APIFD6H3UAf0DXq9flR5XXlqOuGyvUro/Z/Cy/oB0bsOyZ57mnJ7a8Dg8AS099D90rv5Wr8jrRQ1b+q6jGqulxVPwUsxguQH8IL2NuAI/BquV9JMdkrgPeJSB1eQP9VCtfcifcDs0tErks20dSOOP5z4K0i8rSI/FhEzhWRyb7321S1b8T2fqBuxPZRvLwmvikF198BZyVr1eDV9jer6tbk9m/wat7Pi8gvROQtKba7/wzvrmTka3TZPTrOtaP3rwYeVtX40A5VfRLvx3X1iPOaVXVPCm4OH+ICvyNT/Ar4uqruA16H13wxAPwhuZ0KG/G+k38EtqjXfzAhqtoFHI/3Q/Es8GG8wPmq5PHbgSXAd4CipM89k3RGxkdtKyBj7Jsqt+PdTbxDRAqAixjRLKWqTwDLgM8mHX4MPCEiFZOk266qz496je6w7h7n2rH2j5e3kfvHS89hAS7wO6aNiHwAr3liqJaZBwwNzQyQ4vcs+UNxFdBIarX9oesSqnq/qn4VeBXeSKCLRxxvV9VrVfVDeO32a0l2RqfJM8BJo/a9JhVPvPb0dwNnA9V4/QMjz+lS1b+o6sfxftBWJX1zxVbgpJF3GiJyDFCZPOaYARSYFnDYjYgsxOs4HDnO/n7g0yLyHWADcN8UkvwmXk03pSGWInIB3iij+/Fq06/Ca3Z6Jnn8O8BmvKA1CLwDrw18Os0U/wX8SUQeBW7FG8767uSxye4Efos33PU7eO394RF5+RxwAHgC6MG7I0jg3clMRNmI5qMh+lW1LYW8jOaneKOyrhaRy4G5eM1lD6rqA2mk5/AhrsbvmC6/BL6nqrtH7Ps4XufhY3hNAt9INTFV7VfV1mTtOBU68Dqab8MLkP+BN5LnquTxXrwfk81Jn1cC56pqJFWnMRxvBD6P16fxFN6PyVAeJ3y+QFX/iRfY1/DS0UcAUbwfhU3JdN8MvEVVd0yi9AW8u5yRr6n82I70a8a7G1kE/AOv+e1p4C3ppOfwJ6LqVuByOKaLiHwV+ISq1ph2cTgmwzX1OBxTJDm1xGfwhp924435/xze6BqHw/e4Gr/DMUWSI3I24vUnlAO78Jptvp/soHY4fI0L/A6HwzHLcJ27DofDMctwgd/hcDhmGb7v3L333nt1zpw5pjUcDofDKnp6elobGxuDYx3zfeCfM2cOK1euTPv6pqYmGhoaJj/RZ9jqDc7dFM7dDH5137Jly7gzps74pp7CwsLJT/IhtnqDczeFczeDje4zPvBXVlaaVkgLW73BuZvCuZvBRvecBH4RuUpEWkTk6TGOfTa5ElDtWNdOl9bW1mwkm3Vs9QbnbgrnbgYb3XNV478abyWjlyAii4GzmN6EWRNi468x2OsNzt0Uzt0MNrrnJPCr6v1A+xiHfog32VXWniKLx0dPrW4HtnqDczeFczeDje7GRvWIyBuB/ar6pLdU69i0tLSwfv16CgoKSCQSrFu3jg0bNhAKhSgtLSU/P59oNEowGKS9vR1VJRgM0tzcTFlZGe3t7cRiMerq6giHw4gI1dXVhMNhKioqSCQSdHd3U19fTygUorCwkMrKSlpbW6msrCQejxOLxYaPBwIBysvLaWtro6qqilgsRm9v7/DxoqIiiouL6ejooKamhs7OTuLx+PDx4uJiAoEAkUiE2tpaIpEI/f39w8eH8hQKhYb9R+cJoKury7d5amlpeUn6qXxOfslTIpEgFoul/Dn5KU/d3d0UFRVN+7tnIk9Dn0m2/j9lM0+tra1GYsRkeZow/uZqygYRWQpsVNWjRaQE+DtwtqpGRGQ3cLyqvqyxbNOmTTqd4Zx9fX3Y+ByArd7g3E3h3M2QLfcd4W56+wc5ZkF5Wtdv2bJlc2Nj4/FjHTM1qmcF3hJzTyaD/iJgyxiLSUybUCiU6SRzgq3e4NxN4dzNkA33HeFuLr11J1++4wV2tvVkPH0jTT3JtVTnDW1PVOOfLoFAKmtV+w9bvcG5m8K5myHT7kNBvzue4NRlc2moKs5o+pC74ZzX4q0qdKSI7BOR9bl4X4Dy8vRuk0xjqzc4d1M4dzNk0n100P/iGUspyBu/DzRdcjWq5yJVna+qhaq6SFWvHHV8aTZq+wBtbeksO2oeW73BuZvCuZshU+65CvowC57craqqMq2QFrZ6g3M3hXM3Qybccxn0YRYE/lgsZlohLWz1BuduCuduhum65zrowywI/L29vaYV0sJWb3DupnDuZpiOu4mgD7Mg8NfXZ3yEaE6w1RucuymcuxnSdTcV9GEWBH5bxwfb6g3O3RTO3QzpuJsM+jALAn9RUZFphbSw1Rucuymcuxmm6m466MMsCPzFxZl/+CEX2OoNzt0Uzt0MU3H3Q9CHWRD4Ozo6TCukha3e4NxN4dzNkKq7X4I+zILAX1NTY1ohLWz1BuduCuduhlTc/RT0YRYE/s7OTtMKaWGrNzh3Uzh3M0zm7regD7Mg8Nu4SALY6w3O3RTO3QwTufsx6MMsCPy2jg+21RucuymcuxnGc/dr0IdZEPhtHR9sqzc4d1M4dzOM5e7noA+zIPDbOkzMVm9w7qZw7mYY7e73oA+zIPDbusCDrd7g3E3h3M0w0t2GoA+zIPBHIhHTCmlhqzc4d1M4dzMMudsS9GEWBP7a2lrTCmlhqzc4d1M4dzPU1tZaFfRhFgR+W2sStnqDczeFczfD402tVgV9MLTYei7p7+83rZAWtnqDczeFc889O8LdfP/RDmIDak3Qh1lQ47d1fLCt3uDcTeHcc8uu9hiX3rrTuqAPsyDw2zo+2FZvcO6mcO65ozue4Jt37aI7nuC4eYVWBX2YBU09paWlphXSwlZvcO6mcO65QVX54QN72B/tY3l1ER8/scaqoA85qvGLyFUi0iIiT4/Y930R2S4i/xSRv4jI3Gy8d35+fjaSzTq2eoNzN4Vzzw3/+0wr9+86RElhHl9uXEZxoNC00pTJVVPP1cA5o/bdCRytqq8EngW+mI03jkaj2Ug269jqDc7dFM49++wId/PLR/YD8OlTl7Cossga95HkJPCr6v1A+6h9d6jqQHLzYWBRNt47GAxmI9msY6s3OHdTOPfs0tk3wLfv3s3AoHLBUUFOW14F2OE+Gr907r4PuDUbCbe3t09+kg+x1Rucuymce/YYVOX79zXR3BXnyGAJHzhxwfAxv7uPhfHOXRH5EjAA/HGs4y0tLaxfv56CggISiQTr1q1jw4YNhEIhSktLyc/PJxqNEgwGaW9vR1UJBoM0NzdTVlZGd3c3TU1N1NXVEQ6HERGqq6sJh8NUVFSQSCTo7u6mvr6eUChEYWEhlZWVtLa2UllZSTweJxaLDR8PBAKUl5fT1tZGVVUVsViM3t7e4eNFRUUUFxfT0dFBTU0NnZ2dxOPx4ePFxcUEAgEikQi1tbVEIhH6+/uHjw/lqaOjg+rq6jHzBNDV1eXrPCUSiZflaaLPyS95GhwcpKmpKeXPyU956unpIRqNTvu7ZyJPHR0dVFRUZO3/03TzdOPWVh7eE6OkQHj/6hL6e2McTOapq6uLpqamnMeIyfI0YdxV1enG7pQQkaXARlU9esS+9wAfBhpVtWes6zZt2qQrV65M+317e3spKipK+3pT2OoNzt0Uzj07PBXq4nN/e45BhW+evZzXLKl8yXG/um/ZsmVzY2Pj8WMdM9bUIyLnAF8A3jhe0M8Ezc3N2Uo6q9jqDc7dFM4983TE+rn8nt0MKvzrK+e9LOiDf90nIlfDOa8FNgFHisg+EVkP/BQoB+4UkSdE5BfZeO9Ubnv8iK3e4NxN4dwzS2JQ+d7fm2jr6efoulLee/yCMc/zo/tk5KSNX1UvGmP3lbl4b4fD4UiHa54I8fiBTiqLCrjsTLuezJ0Mv4zqyRpdXV2mFdLCVm9w7qZw7pljy/4of9gSQoAvntFAben4C8X4zT0VZnzgr6urM62QFrZ6g3M3hXPPDK3dcb779yYUeOdx9Ry3sGLC8/3kniozPvCHw2HTCmlhqzc4d1M49+mTGFQuv2c3kd4Bjl1QzsVrJp811C/uU2HGB34RO9vlbPUG524K5z59rn7sAE83d1NTUsilZzSQn0K7vl/cp8KkgV9E8kVkp4jMyYVQpqmurjatkBa2eoNzN4Vznx4P74nwp3+2kCdw2ZlLqSpObfI1P7hPlUkDv6omgATgvycUUsDG2zCw1xucuymce/qEOvv4/n1NALzv+AW8oj71IZqm3dMh1eGcPwKuF5HLgX3A8OO+qvpCNsQyRUXFxB0zfsVWb3DupnDu6dGfGOQ79+ymsy/BiYsreOsr503pehvLPdXA/9Pkv2eN2q+AryfSTiQSphXSwlZvcO6mcO7pccUjB9gR7qGuLMDn1jaQN8U2exvLPaXOXVXNG+fl66AP0N3dbVohLWz1BuduCuc+de5/oYObnwlTkCd86cylVBRN/ZlWG8t9SrkUkSXAQmCfqu7NjlJmsXERZ7DXG5y7KZz71NgX6eUHD+wB4IMnLmTlvPSWf7Sx3FOq8YvIfBG5D3geuBHYKSL3i8jYk1f4CNsWcR7CVm9w7qZw7qnTNzDIt+/eRU//IKctm8sFR9WmnZaN5Z7qOP7/Bp4EqlR1PlAFPA5kZWK1TFJYaN96mGCvNzh3Uzj31Pn5pn280N7Lgoo5fOrUJdMai29juafa1HMKMF9V+wFUtVtEPg/sz5pZhqisfPk0qjZgqzc4d1M499S487k2bt3RRiBf+ErjUkoD0+uqtLHcU63xdwBHjdp3JHAoszqZp7W11bRCWtjqDc7dFM59cna1x/jJg1735IbXLmZFTcm007Sx3FOt8f8HcJeIXAk0AQ3AJcBXsiWWKWz8NQZ7vcG5m8K5T0ysP8G3795FX0J53eHVnHNEZp64tbHcUwr8qvorEdkJXAy8EjgAXKSq92RTLhPE43HTCmlhqzc4d1M49/FRVX704F72RvpoqCriY69dlLE5dmws90kDv4jkA1cBH7Qh0I8mFouZVkgLW73BuZvCuY/P37a38fedHRQV5PGVM5dRXJi5R5BsLPdU5+o5GxjMvk7msXGMLdjrDc7dFM59bJ5r7eG/N+0D4JOnLGZJVWanHbOx3FPt3P0h8A0RsW7cko1jbMFeb3DupnDuL6erb4Bv372L/kHlvJU1nHlY5mfStLHcU+3c/RhQD3xaRMK8dJK2JdkQyxSBwPhLpvkZW73BuZvCub9IfGCQe1/o4IanWjjYGeewmmI+8ppFGX2PIWws91QD/zuzapFFysvLTSukha3e4NxN4dwh3B1n4zOt3LKjjUjvAAC1JYV8uXEZgYLsrDtlY7mn2rn7PrzO3b7sK2WWtrY2yspSn1vbL9jqDc7dFLPVXVV5urmbm7eGeXD3IQaT7RGH1RTzptVBTl9elbWgD3aW+6SBX1UTIjKtzl0RuQo4H2hR1aOT+6qBPwFLgd3Av6pqR7rvMR5VVVWZTjIn2OoNzt0Us829b2CQv+/0Ztfc2eaNrMkTWLt8Lm86KshRdaU5WRbRxnLPVefu1cA5o/ZdCtytqocDdye3M46NQ63AXm9w7qaYLe4tXXGu/McB3nHt0/zggT3sbItRWVTAxWvq+MPbV/OlM5exur4sZ2vh2ljuOencVdX7RWTpqN0XAKcn//4tcC/whRR9Uqa3tzfTSeYEW73BuZtiJrurKk+Furhpayv/1/Ric87htV5zztpl2W3OmQgby91k526dqh4EUNWDIjK19c5SxMYxtmCvNzh3U8xE996BQf7+fDs3PxPmhXYvwOYLnLGiiguOCrJqXknOavbjYWO5pzplw33ZFhmPlpYW1q9fT0FBAYlEgnXr1rFhwwZCoRClpaXk5+cTjUYJBoO0t7ejqgSDQZqbmykrKyMcDlNcXExdXR3hcBgRobq6mnA4TEVFBYlEgu7uburr6wmFQhQWFlJZWUlrayuVlZXE43Fisdjw8UAgQHl5OW1tbVRVVRGLxejt7R0+XlRURHFxMR0dHdTU1NDZ2Uk8Hh8+XlxcTCAQIBKJUFtbSyQSob+/f/j4UJ727t3LYYcdNmaeALq6unybp127dlFeXv6yPE30OfklT4lEgvz8/JQ/Jz/lqbu7m0WLFk37u2ciTy+88AJLly4dzlNeWTXXP76Ph/b30dXvVe8r5+Rz+pIiTqkvZOXS+YRCIcJh83k6cOAApaWlOY8Rk+VpIkRVxz8ocrOqXjBi+xuq+rUR2/9Q1VenEsCTTT0bR3Tu7gBOT9b25wP3quqRo6/btGmTrly5MpW3GJPm5mbq6urSvt4UtnqDczeF7e7z5s3jyYNd3Lw1zKY9keHmnCODJVxwVJDTls8lkG+mOWci/FruW7Zs2dzY2Hj8WMcmq/GfMWr7Y8DXRmynH5Hhf4H3AN9L/nvzNNIal+Li4mwkm3Vs9Qbnbgpb3WP9CR480M/tD25nd4fXnFOQJ5y+fC4XrA6yKs0lEXOFjeU+1ZWFRzemjX+7MPIikWvxOnJrRWQf3o/H94DrRWQ9sAd42xRdUqKjo4OKiopsJJ1VbPUG524KG913tce49Nbn6Yh5D1tVFxdw3qpazltZS3WJHTPE2FjuUw38KQX6l12ketE4hxrTSW8q1NTUZPstsoKt3uDcTWGbe0tXnC/dtpOO2ADLquZw4TH1nLpsLoU+bM6ZCNvKHSYP/IUicgkv1vTniMj7pnC9cTo7O617qg7s9Qbnbgqb3Dv7BvjS7Ttp7enn6PpSPrGmjIZFmZ9ALRfYVO5DTBa4HwHePWL7UeBdo477GhsXSQB7vcG5m8IW9/jAIF+/cxdNHb00zC3iG2ctpz3k++W7x8WWch/JhIFfVU/PkUfWsHGMLdjrDc7dFDa4D6ry7/c18VSoi5qSQr5zzgrK5xQQsMB9PGwo99HY1ZiWBjbOlQ32eoNzN4Xf3VWVXzy8nwd2HaKkMI/vvH4F88q8KY397j4RNrrP+MBv41ArsNcbnLsp/O7+P0+1cNPWMIV5wtfPWs7ymhd9/e4+ETa6z/jAb+MiCWCvNzh3U/jZ/e7n2/n1owcA+NzaBtYseOkc9n52nwwb3Wd84I9EIqYV0sJWb3DupvCr+5b9Uf7r/j0AfOjEhZy+4uXTGPvVPRVsdE95OKaIrALeCtSr6gYRWQkEVPWfWbPLALW1taYV0sJWb3DupvCj+862Hr551y4GBpW3HB3kLa8Yey5GP7qnio3uKdX4ReRtwH3AQl4czlkG/CBLXhnDxl9jsNcbnLsp/OYe6uzjS7ftpKd/kLXL5/KBExeOe67f3KeCje6pNvV8EzhbVT8MJJL7ngSOyYpVBunv7zetkBa2eoNzN4Wf3KO9A1x2207aYwMcM7+Mz61tIG+C6ZP95D5VbHRPNfDPwwv08OK0DUqaUzjkEhvH2IK93uDcTeEX976BQb56xwvsi/SxrKqIr5+1fNJZNf3ing42uqca+Dfz0id2Ad6O9ySvr7FxjC3Y6w3O3RR+cE8MKpf/fTfPtHQTLPUe0CoN5E96nR/c08VG91Q7dz8O3JGcSbNURG4HjgDOzppZhigt9feUruNhqzc4d1OYdldVfrZpH5uaIpQF8rn8nBXUlqY21NG0+3Sw0T3VFbi2J0fxnA9sBPbiLarSlU25TJCfP3ltw4/Y6g3O3RSm3a97spmN21opzBe+cfZyGqpSf7DJtPt0sNE91VE9P1HVHlW9XlW/r6rXqWqXiPwo24LTJRqNmlZIC1u9wbmbwqT7Hc+28ZvHDiLApacv5RX1U5ut0pV7bkm1jf+94+wf3e7vO4LBoGmFtLDVG5y7KUy5/2NvlB884D2g9dGTFnHqsrlTTsOVe26ZsKlnxNz7BaPm4QdYDrRmxSqDtLe3U1JSYlpjytjqDc7dFCbcn23t4Vt372JQ4cJj6rhgdXpB0JV7bpmsjX+oRh/gpbV7BZrx1sr1NRMtJu9nbPUG526KXLsfiPbx5dt20jswyOsOq+J9x89POy1X7rllsvn4zwAQkW+r6pdzo5RZbLwNA3u9wbmbIpfuh2L9XHbbTg71DnDcwnI+deoSZIIHtCbDlXtuSbWN/6sikjfWK6t2GaC5udm0QlrY6g3O3RS5co/1J/jKHS9wINrHYTXFfLVx2bTXyXXlnltSHcc/wPhP6fp6LJNta2EOYas3OHdT5MI9Mahcfs9udoR7qCsL8O3Xr6AkhQe0JsOVe25JNfAvG7U9H7gU+GtmdRwOh19RVX7y0F4e2RulYo73gFZ1SaFpLUcapHR/pqpNo14P43XsfiG7etOnq8v3z5iNia3e4NxNkW33Pzwe4tYdbczJF771+hUsnluUsbRdueeW6TTMVQDT7tUQkU+JyFYReVpErhWRzH2bgLq6ukwmlzNs9Qbnbopsut+yvZXfbwmRJ3DZmctYNS+z0xS4cs8tqT65+3sR+d2I1w14E7f9YTpvLiIL8eYBOl5Vj8brL3j7dNIcTTgczmRyOcNWb3DupsiGu6py53Nt/OShvQD822sXc1JDZcbfx5V7bkm1jf/5UdvdwC9U9a4MORSLSD9QAhzIQJrDTGeImUls9QbnbopMux+K9fOTh/bx4O5DALzj2HrOX5Wd1aZcueeWVCdp+0Y23lxV94vIfwJ7gBhwh6reMfKclpYW1q9fT0FBAYlEgnXr1rFhwwZCoRClpaXk5+cTjUYJBoO0t7ejqgSDQZqbmykrKyMvL4+mpibq6uoIh8OICNXV1YTDYSoqKkgkEnR3d1NfX08oFKKwsJDKykpaW1uprKwkHo8Ti8WGjwcCAcrLy2lra6OqqopYLEZvb+/w8aKiIoqLi+no6KCmpobOzk7i8fjw8eLiYgKBAJFIhNraWiKRCP39/cPHh/LU29tLT0/PmHkCr13Rr3kaGBigqanpZXma6HPyS54qKytpampK+XPyU57y8vKIRqPT/u5Fo1Fe6CviZw8fJBpXigvyeNvhRZzTUEhbW1tW8tTb20s0Gs3a/6dsfk7AS77vuYoRk+VpImS8p87GmKJhvOB9VSrnjfMeVcCfgQuBQ8D/ADeo6nAT0qZNm3TlypXpvgVNTU00NDSkfb0pbPUG526KTLh39Q3w84f3c9dz7QAcM7+Mz5y2hPryOZlQHJfZXu7ZYMuWLZsbGxuPH+vYRDX+VCZgUyDtwA+8DtilqmEAEbkReC3T7DsYSUVFRaaSyim2eoNzN8V03R/b50221trdTyBfWP/qBVywOjjhkomZYjaXuwnGDfxD0zVkmT3Aa0SkBK+ppxF4LJNvkEgkJj/Jh9jqDc7dFOm6x/oT/OrRA2zc5s25uDJYwufWNmR0uOZkzMZyN0nKwzlFpEpE3i0iX0z+WzXdN1fVR4AbgC3AU0mfK6ab7ki6u7szmVzOsNUbnLsp0nF/KtTFh2/czsZtrRTkCZccP58fvuGInAZ9mH3lbpqUOndF5CTgb8B2oAlvJa4fich5qrppOgKq+jXga9NJYyJsXAgZ7PUG526KqbjHBwa5evNB/vxUCwosry7m82sbWF6T+qpZmWS2lLtfSLXG/yPgo6r6WlW9SFVPBj4C/CR7apnBxoWQwV5vcO6mSNX92XAPH71pBzc81YIIXLSmjv93wRHGgj7MjnL3E6mO4z8CuH7UvhuAX2RWJ/MUFto5l4it3uDcTTGZ+8Cgcs3jIa55IsSgwqLKOXxubUPGn8JNh5lc7n4k1cD/HN4TtdeM2Pc2YGfGjTLM0Dhb27DVG5y7KSZy39Ue4/v3NfF8WwyANx8d5H3HL2BOgT9mVp+p5e5XUg38nwQ2isjH8dr4lwKH47X1+5rW1lZKS83XaKaKrd7g3E0xlujryoMAACAASURBVHtiUPnzUy38dvNB+geVurIAnz1tCccsKDdkOTYzrdz9TqpP7v6fiKwAzgMW4E3HfIuqtmdTLhPY+GsM9nqDczfFaPf9kV6+f98enmnxRp2ce2QNHzpxYUbmz880M6ncbSDVGj+q2kHywSoRWQ6UA74P/PF43LRCWtjqDc7dFEPug6r89ZlWfv3ofvoSSnVJAZ8+dQknLPZvgJoJ5W4Tqc7Oea2IvDb59yXAVuAZEVmfTblMEIvFTCukha3e4NxNEYvFaOmK88Vbn+dnm/bRl1DOXFHFFetW+Trog/3lbhup1vgb8RZeAfg03lQLh4CbgCuz4JUxbBxjC/Z6g3M3garydFcRv7p3Gz39g1QWFfCxkxdx2rJpP2eZE2wtd7DTPdUu/YCqxpPz51er6kOquhXw/QoENo6xBXu9wbnnmkOxfr5+5y5+vOkAPf2DnNRQyRXrVloT9MHOch/CRvdUa/xPiMgXgQa8J3iHFlGJZkssUwQCAdMKaWGrNzj3XPLYvij/eV8T7bEBSgqEDScv5nWHVVs3R7xt5T4SG91TDfzrgW8B/cDnk/tOAv6YDalMUl7ur2FrqWKrNzj3XBBPDPKbfxzgz097qz8dXV/Kx14dZFmdPbX8kdhS7mNho3uqwzl3AheP2ncD3tO7vqatrS2lhQn8hq3e4NyzzZ6OXr577252tsXIE3j3cfO58Jg69u3dA9gZ+G0o9/Gw0T3l4ZzJhVkuwhvHfwC4DrhKx1vJxSdUVdn5H8FWb3Du2UJV+dv2Nn75sDdiZ355gEvPWDo85YKf3SfDueeWVGfn/A/gArzJ2pqAJcBngSN5senHl8RiMSsXSrDVG5x7Noj0DvCDB/awqSkCwOsOr2bDSYsoHfEwll/dU8G555ZUa/zvBY5T1X1DO0Tkb3jz6Ps68Pf29ppWSAtbvcG5Z5rH93fyH/c10dbTT0lhHp84ZTFnrKh+2Xl+dE8V555bUg38ncnX6H2+H9Vj4xhbsNcbnHum6E8McvVjB7khOWf+6rpSvnB6w7jr3/rJfao499wy7jh+EVk+9MJr4rlRRM4SkVUicjbewug/zJVoutg4xhbs9Qbnngn2RXr55F+f5X+Sc+a/67h6/vO8wydc9Nwv7ung3HPLRDX+5/EWUx85IHj0OrxnAj/NtFQmKSrK7RJymcJWb3Du00FVuW1HGz9/eD99A4PUlQW49IwGVtdNPmrEtPt0cO65ZaLF1v0xUfc0KS42t6rQdLDVG5x7ukR7B/jRg3t5cPchAM5cUcXHTl78kg7ciXDlbgYb3dMO7iKyOjnax9d0dHSYVkgLW73BuafDkwc6+fBftvPg7kOUFObx+bUNXHrG0pSDPrhyN4WN7imP4wcQkVq8B7neDRwL3JoNqUxSU1NjWiEtbPUG5z4VBgaV320+yJ+ebEaBVfNKuPT0pcyvGL8tfzxcuZvBRvdJA7+IFAJvwAv25wJ78R7ierWqbsmu3vTp7Oy07qk6sNcbnHuq7I/08b17d7Mj3EOewMVr6nnnsfXk56U3z44rdzPY6D5h4BeRn+KttduPN4pnrao+LCIHgX0TXZsqIjIX+DVwNF5n8vtUdVMm0gY7F0kAe73BuU+GqnLnc+38bNM+Yv2DzCsr5AunL+UV9dMLHq7czWCj+2Q1/o/grbL1deA6VY1kweHHwG2q+lYRCQAlmUzcxjG2YK83OPeJ6Oob4McP7uW+XV4H7trlc/nEyYspmzOlVtcxceVuBhvdJ+vcXQH8P+BzQEhE/iwib0nhupQQkQrgNJKLuahqXFUPZSLtIWwcYwv2eoNzH4+nQl18+C/buW/XIYoL8/jsaUu47IylGQn64MrdFDa6T/iNU9XdwDeBb4rIqXjt/FcCFcB3ROSHqvrMNN5/ORAGfiMixwCbgU+oavfQCS0tLaxfv56CggISiQTr1q1jw4YNhEIhSktLyc/PJxqNEgwGaW9vR1UJBoM0NzdTVlZGIpGgqamJuro6wuEwIkJ1dTXhcJiKigoSiQTd3d3U19cTCoUoLCyksrKS1tZWKisricfjxGKx4eOBQIDy8nLa2tqoqqoiFovR29s7fLyoqIji4mI6Ojqoqamhs7OTeDw+fLy4uJhAIEAkEqG2tpZIJEJ/f//w8aE8dXd309PTM2aeALq6unybp97eXpqaml6Wp4k+J7/kqbCwkKamppQ/p1TyFKio4apH9nL/vl4UWFqRz6dOqidYrOzZsydjeYrH40Sj0Wl/90x8Tt3d3USj0az9f8pmnvr7+1/yfc9VjJgsTxMhU51cU0SKgDfjLcV4pqqmvQqBiBwPPAycrKqPiMiPgaiqfmXonE2bNunKlSvTfQsOHTrE3Llz077eFLZ6g3MfItaf4M9PtXD9P1voHRgkT+DCV9bxrlfNpyDNDtyJcOVuBr+6b9myZXNjY+PxYx2bcpONqvaq6rWqeg6wdJpu+4B9qvpIcvsG4LhppvkSIpFsdEtkH1u9wbknBpVbtrdyyfXP8LstIXoHBjlpSSVXvGUVl7x6QVaCPrhyN4WN7tNqXFTVA9O8PiQie0XkSFXdgbeo+3Sajl5GbW1tJpPLGbZ6w+x1V1Ue3Rvl148eoOmQN2PjkcESPnDCQl45P/vD/WZruZvGRvfM9CpNj48Bf0yO6HkBuCSTiUciEUpLSzOZZE6w1Rtmp/uzrT386pH9PHmwC4D68gCXHL+Atcvnkpej9W9nY7n7ARvdjQd+VX0CGLMdKhP09/dnK+msYqs3zC73UGcfVz92kHt2eo/tl8/J5+I19bzhqFoC+bmd7mo2lbufsNHdeODPNjaOsQV7vWF2uHf2DXDtE83cvDVM/6BSmCdcsDrIRWvqKM/Q8MypMhvK3Y/Y6J5SlUREqkXkchG5RUTuH/nKtuB0sXGMLdjrDTPbPZ4Y5M9PtfDe65/hhqda6B9UzlhRxZVvW8UHT1xoLOjDzC53P2Oje6rf0muAOcD1QE/2dDKPbW1vQ9jqDTPTXVW574VDXPXYAUKd3iP6x8wv4wMnLOSIYEYfNk+bmVjuNmCje6qB/7VAUFX7simTDfLzU5/W1k/Y6g0zz/2pUBdXPLKfHWGvztMwt4j3n7CAExZXIDnquE2FmVbutmCje6qB/5/AImBnFl2yQjQapaqqyrTGlLHVG2aO+95Dvfz6HwfY1OSN064uLuDdr5rP64+oSXsGzWwyU8rdNmx0TzXw3wPcJiK/AV7SoKWqV2XcKoMEg0HTCmlhqzfY797R08/vHw9xy/ZWBhWKCvJ42yvn8dZXzKO40L+1O9vL3VZsdE818J+K95TtWaP2K+DrwN/e3k5JiT/aYKeCrd5gr3vvwCC//cdebtvdS6zfm2LhX1bW8K7j5lNTUmhab1JsLXdw7rkmpcCvqqMXWbeGqc5F5Bds9Qb73AdVueu5dn7z2EHaerwx2ScuruD9Jyygocqe9VRtK/eROPfcMuWxZ+L1Zg03cKrqYEaNMoyNt2FgrzfY5f7EgU6ueGQ/z7fFADisuogPvWYRxywoN2w2dWwq99E499yS6jj+hSLyFxFpAwbwVuQaevma5uZm0wppYas32OG+91AvX71jJ5+/5Xmeb4tRW1rI59c28Lnjiq0M+mBHuY+Hc88tqdb4f4E3fr8RuA9v8ZSvA7dkRytz2LYW5hC2eoO/3SO9A/xhy0E2bmsloVBcmMeFr6xj3SvmUVSQR1ubfbftQ/i53CfDueeWqYzjX6Kq3SKiqvqkiKwH/g/4Vfb0HI7MEE8McvPWMNc80Ux3PEGewLlH1vCeV82n2oKOW4cjk6Qa+BN4TTwAh0QkCESBhVmxyiBdXV3U1NSY1pgytnqDv9xVlQd2HeLX/3jxidtXLSzngycuZFn1yztu/eQ+VZy7GWx0TzXwPwL8C/AX4HbgT0AMeCxLXhmjrq7OtEJa2OoN/nHf1tLNLx/ezzMt3kqeDVVFfPCEhbx6ccW41/jFPR2cuxlsdE913th34bXtA3wS74Gup4GLsyGVScLhsGmFtLDVG8y7hzr7uPyeXXzif5/lmZZu5hYV8PGTF/OLN6+cMOiDeffp4NzNYKN7quP4D434OwZ8O2tGGcZPc6lMBVu9wZx7dzzBdU+EuHFrmP6EUpgvvOXoeVx4TB2lgdSeuHXlbgbnnltSCvwiMgf4KnARUKOqlSJyNnCEqv40m4LTpbq62rRCWtjqDbl3Twwqf9veyu+3hIj0el1RZ66o4n2vXsC8ssCU0nLlbgbnnltSber5IXA08A68aRoAtgIfyYZUJrHxNgzs9Ybcuasqj+yJ8KEbt/PT/9tHpHeAo+tK+ckbj+DSM5ZOOeiDK3dTOPfckmrn7puBw5LDOQcBVHW/iPh+VE9FxcRtun7FVm/IjfvOth6ueOQAjx/oBGBBRYD1r17IKUsrp3Xr7crdDM49t6Qa+OOjz00O6WzLuFGGSSQSphXSwlZvyK57W08/Vz92gDuebUeBskA+Fx9bzxsztMatK3czOPfckur/lP8BfisiywBEZD7wU+C6bIlliu7ubtMKaWGrN2THvakjxo8f3MN7/7SV259tJ0/gzauDXP2vR/HWV8zL2MLmrtzN4NxzS6o1/suA/wCeAkqA5/Ce2P1Glrwyho0LIYO93pA590FVHt0b5S9Ph4ebdABe21DJ+09YwKLKooy8z0hcuZvBueeWVIdzxvHG738y2cTTqhmci1RE8vEeBtuvqudnKl3wFkJuaGjIZJI5wVZvmL57dzzBHc+2cfMzYQ5Evadt5xTkcdbh1VxwVG1Wp0qezeVuEueeWyYM/CKyZJxDi4c60FR1TwY8PgFsAzLeS1JYaOc8LLZ6Q/ru+yO93LS1lTueayPW7832XVcW4I1H1XLOkTWUz5nyLOJTZjaWux9w7rllsv9Ju3lx+OZYQyUUmNZadCKyCDgP+A7w6emkNRaVlZWZTjIn2OoNU3NXVTbv7+SmrWEe3Rsd3n/M/DIuWB3kpCWVOV3fdraUu99w7rllssD/T6AI+C3wB+BAFhx+BHweyMok6K2trZSWlmYj6axiqzek5h7rT3DXc+3ctDXM3kgfAIX5QuOKai5YXcuKGjNL2c30cvcrzj23TBj4VXWNiBwNvAd4ENgO/A64MTl1w7QQkfOBFlXdLCKnj3VOS0sL69evp6CggEQiwbp169iwYQOhUIjS0lLy8/OJRqMEg0Ha29tRVYLBIM3NzZSVlSEiNDU1UVdXRzgcRkSorq4mHA5TUVFBIpGgu7ub+vp6QqEQhYWFVFZW0traSmVlJfF4nFgsNnw8EAhQXl5OW1sbVVVVxGIxent7h48XFRVRXFxMR0cHNTU1dHZ2Eo/Hh48XFxcTCASIRCLU1tYSiUTo7+8fPj6Up97eXnp6esbME3gzAvo1TwMDAzQ1Nb0sT9FolMHiSm78Z4j79/USG/BuJquL8jlreRknBoXDFs8jHA6zL9ZuJE8lJSU0NTWl/DlN9N3L9ec0ODhINBqd9nfPRJ56e3uJRqNZ+/+UzTyp6ku+77mKEZPlacLYm2ofrYjk4S22/l7gXOBMVd2S0sXjp/ldvAngBvDuLCrwflTeOXTOpk2bdOXKlWm/R0tLC/PmzZuOphFs9YaXu6sqTx7s4qatYR7eE2Ew+ZU7al4pb1od5JRlcynIYXPORMykcrcJ5555tmzZsrmxsfH4sY5NpbfscGAtcBLwONAxXTFV/SLwRYBkjf+zI4N+JojFpn1jYgRbveFF976BQe553mvO2dXRC0BBnnDmirm8afU8jgiaac6ZiJlQ7jbi3HPLZKN6qvEmZnsPXhv874HTMjSSJyfYOMYW7PUGyC+v4cp/HOCW7a109nlPNVYVF3D+qlrOW1nr6xWvbC53524GG90nq/EfAHbhBfyHk/sOE5HDhk5Q1XsyIaKq9wL3ZiKtkdg4xhbs8x5qzvnrtlYe2n1ouDnniNoS3rQ6yGnL52bs6dpsYlu5j8S5m8FG98kCfwiv7f0DyddoFFieaalMEghMfYZGP2CLd1ffAHc+187Gba3Do3PyBNYun8ubV89j1bwSq+Yrt6Xcx8K5m8FG98lG9SzNkUfWKC/PyijRrON37+dbe/jrtlbu2dlB34D3sFVNSSH/srKG0xYV0zBvrmHD9PB7uU+EczeDje7ZfxTSMG1tbSkNb/IbfvSODwxy/65D/HVbmG0tPcP7j11QxvmrgpzUUElBnjd8FuwM/H4s91Rx7maw0X3GB/6qqirTCmnhJ++D0T42bmvl9mfbiCY7a0sD+Zx9eDXnraplydyXTpbmJ/ep4tzN4Nxzy4wP/LFYzMqFEkx7Jwa9mTE3bmvlsX3R4Xk7Dqsp5g1HBTl9+VyKC8eercO0+3Rw7mZw7rllxgf+3t5e0wppYcq7I9bPbTvauGV7G81d3syYhfnC6curOH9VLSuDk3fW2lrm4NxN4dxzy4wP/DaOsYXceqsqW5u7+eu2Vh7YdYiB5FjMBRUBzltZy+uPqKGiKPWviq1lDs7dFM49t8z4wG/jGFvIjXdPPME9Ozv46zMvPlmbJ3DSkkrecFQtxy0sJy+NoZi2ljk4d1M499wy4wN/UVHmV2nKBdn03t0R46/PtHL38+30JOe9n1tUwLkrazhvZS3zyqY3LtnWMgfnbgrnnltmfOAvLs7eak3ZJNPe8cQgD+46xMbtrTwdenGN0FfUl3H+qlpOWVpJYYaerLW1zMG5m8K555YZH/g7Ojqs63GHzHkf7Ozjlm2t3PZsO5HeAQBKCvNoPKya81fVsqw6819aW8scnLspnHtumfGBv6amxrRCWkzHe7yhmMurizl/VS1nrqiiJDCthdMmxNYyB+duCueeW2Z84O/s7LTuqTpIz7utZ2goZivh7n7AG4q5dnkVb0hxKGYmsLXMwbmbwrnnlhkf+OPxuGmFtEjVW1V54mAXG7e18n+7D5FIVu8XVMzh/JU1nD3FoZiZwNYyB+duCueeW2Z84LdxjC1M7t05YlbMfSNmxTy5oZLzV9VybJpDMTOBrWUOzt0Uzj23zPjAb+MYWxjbW1XZEe5h47ZW7n2hg3iyel9bUsi5K2s498gaakvNTxFra5mDczeFc88tMz7w2zjUCl7qHetP8PedHWzc1srzbS8u8/aqheWct6qWk5ZUku+TNWvB3jIH524K555bZnzgt3GRBPC8d3fE+Nu2Vu587sUHrcrn5PP6I7wHrRZWzjFsOTa2ljk4d1M499wy4wN/JBJh7lx75oZv6+lnU1OE27Ye5NlDA8P7j5pXyvmrajlt2VwCBf5ewtC2Mh+JczeDc88tMz7w19bWmlaYlAPRPh7afYiHdkfY1tI9PO6+uDCPxhXVnLeqhhU1JUYdp4INZT4ezt0Mzj23zPjAH4lEKC0tNa3xElSVF9pjPLQ7wkO7Dw1PkAbeuPtXLSxndSW84dilWX3QKlv4scxTxbmbwbnnlhkf+Pv7+00rADCoyrbmbh5q8oL9wc4Xx/6WFOZx4pJKTm6o5PhFFZQE8mlqarIy6IN/yjwdnLsZnHtuMRr4RWQx8DugHhgErlDVH2fyPUyOse1PDPLkwS4e3H2ITU0ROmIvttnPLSrgpIZKTl5ayZoF5QRGTZBm49jgIZy7GZy7GWx0N13jHwA+o6pbRKQc2Cwid6rqM5l6g1yPsY31J3hsXycP7j7Eo3ujdMcTw8fqygKcvLSSU5bOZdW80gmHYNo4NngI524G524GG92NBn5VPQgcTP7dKSLbgIVAxgJ/Ltreor0DPLwnwkO7I2zeHx1+sApgaVURpyydy8lLK1leXZzyXDm2tRmOxLmbwbmbwUZ30zX+YURkKXAs8Egm083Pz047eUdPP/fvOsRDTYf458EuBl+M9ayaV8LJS+dycsPctMfaZ8s7Fzh3Mzh3M9jo7ovALyJlwJ+BT6pqdOSxlpYW1q9fT0FBAYlEgnXr1rFhwwZCoRClpaXk5+cTjUYJBoO0t7ejqgSDQZqbmykrKyMcDhONRqmrqyMcDiMiVFdXEw6HqaioIJFI0N3dTX19PaFQiMLCQiorK2ltbaWyspJ4PE4sFhs+3j2Yz517+7n9+UMMeM9UkS+wZn4pqyuVExeXs6CqjI6ODirz+zl4sJ14PD58fXFxMYFAgEgkQm1tLZFIhP7+/uHjQ3nau3cvc+bMGTNPAF1dXRnLUyAQoLy8nLa2NqqqqojFYvT29g4fLyoqori4mI6ODmpqaujs7JwwT/v37ycajb4sTxN9Tn7JUyKRIBqNpvw5+SlP3d3d5Ofnp/w5+SlPe/fuJT8/f9rfPRN5Onjw4Eu+75n+/5RuniaMuao66UnZREQKgY3A7ar6g9HHN23apCtXrkw7/Z6eHkpKpj8GvrU7zp+ebOaWHW30J5tyTlxcwdrlVZy4pILyOZn9Dc2Utwmcuxmcuxn86r5ly5bNjY2Nx491zPSoHgGuBLaNFfQzQXt7+7Q+lHAy4N+6vY3+ZHvOqcvm8o419Syvyd4cHdP1NolzN4NzN4ON7qabek4G3gU8JSJPJPddpqq3ZOoN0r2jaenyAv5tO7yAL8Bpy+byjmPrs7Jc4WhM34lNB+duBuduBhvdTY/qeRDI6rSSwWBwSue3dMW57slmbh8R8Ncum8vFOQr4Q0zV2084dzM4dzPY6O7v2b4yQHNzc2rndcb58YN7eO/1z7BxWysDg8rpy+fyy7es5EuNy3Ia9CF1bz/i3M3g3M1go7vppp6sM1kPd3NnnGufDHHHs+0MJGv4Z6yo4uI1dTRUmZtn27Y1PEfi3M3g3M1go/uMD/zjEers49onmrnj2TYS6i1beMaKKt6xpp4lVUWm9RwOhyNrzPjA39XVRU1NzfD2wc4+rhsV8M9cUcXFx9azZK5/Av5ob5tw7mZw7maw0X3GB/66ujoADkb7uOaJEHc91z4c8BsPq+LiNfUs9lHAH2LI20acuxmcuxlsdJ/xgX9rU4h7m4U7n2tnMBnwX3d4NRevqWNRpf8C/hDhcJjFixeb1kgL524G524GG91nbOA/EO3jmsdD3PV8x3DAPysZ8Bf6OOAPkepkbn7EuZvBuZvBRvcZG/g374tyx3Pt5AmcfXg1F62p9+3i5GNRXV1tWiFtnLsZnLsZbHSfseP4X39kDW8+Osi3XlvJZ9c2WBX0wbt9tBXnbgbnbgYb3WdsjT+Qn8dHXrOIjo4O0yppUVFRYVohbZy7GZy7GWx0n7E1/iESicTkJ/kQW73BuZvCuZvBRvcZH/i7u7tNK6SFrd7g3E3h3M1go/uMD/w2LoQM9nqDczeFczeDje4zPvCHQiHTCmlhqzc4d1M4dzPY6D7jA/9NN91kWiEtbPUG524K524GG91nfOC/8cYbTSukha3e4NxN4dzNYKP7jA/8AwMDphXSwlZvcO6mcO5msNHd+GLrk3H33XeHgaZ0r29vb6+trq5uzaBSTrDVG5y7KZy7GXzs3tDY2Djm8mC+D/wOh8PhyCwzvqnH4XA4HC/FBX6Hw+GYZcyYwC8iV4lIi4g8PWJftYjcKSLPJf+tMuk4HiKyWET+LiLbRGSriHwiud/3/iJSJCKPisiTSfdvJPf73h1ARPJF5HER2ZjctsV7t4g8JSJPiMhjyX22uM8VkRtEZHvyO3+SDe4icmSyvIdeURH5pA3uo5kxgR+4Gjhn1L5LgbtV9XDg7uS2HxkAPqOqq4DXABtE5Cjs8O8DzlTVY4A1wDki8hrscAf4BLBtxLYt3gBnqOoaVT0+uW2L+4+B21R1JXAMXvn73l1VdyTLew3wKqAH+AsWuL8MVZ0xL2Ap8PSI7R3A/OTf84Edph1TzMfNwFm2+QMlwBbgRBvcgUV4/1HPBDba9J0BdgO1o/b53h2oAHaRHFhik/so37OBh2x0V9UZVeMfizpVPQiQ/HeeYZ9JEZGlwLHAI1jin2wueQJoAe5UVVvcfwR8Hhgcsc8GbwAF7hCRzSLyweQ+G9yXA2HgN8kmtl+LSCl2uI/k7cC1yb9tc5/xgd8qRKQM+DPwSVWNmvZJFVVNqHf7uwg4QUSONu00GSJyPtCiqptNu6TJyap6HHAuXtPgaaaFUqQAOA74b1U9FujGhqaREYhIAHgj8D+mXdJlpgf+ZhGZD5D8t8Wwz7iISCFe0P+jqg49A26NP4CqHgLuxetr8bv7ycAbRWQ3cB1wpoj8Af97A6CqB5L/tuC1M5+AHe77gH3Ju0KAG/B+CGxwH+JcYIuqNie3bXIHZn7g/1/gPcm/34PXdu47xFut+Upgm6r+YMQh3/uLSFBE5ib/LgZeB2zH5+6q+kVVXaSqS/Fu2+9R1Xfic28AESkVkfKhv/Ham5/GAndVDQF7ReTI5K5G4BkscB/BRbzYzAN2uQMz6MldEbkWOB2oBZqBrwE3AdcDS4A9wNtUtd2U43iIyCnAA8BTvNjefBleO7+v/UXklcBvgXy8isT1qvpNEanB5+5DiMjpwGdV9XwbvEVkOV4tH7ymk2tU9Ts2uAOIyBrg10AAeAG4hOR3B/+7lwB7geWqGknus6LcRzJjAr/D4XA4UmOmN/U4HA6HYxQu8DscDscswwV+h8PhmGW4wO9wOByzDBf4HQ6HY5bhAr/D14jIL0TkK9O4/jIR+XUmnfyGiDwkIsemcV1dcnbMOdnwcvgXF/gdRkhOKxwXkdpR+58QEU3OWYSqflhVv5Xu+6jq5ar6/unZjk3Ss1tEukSkTUTuFpELp3D96SKyb5oObwA6VfXx5HajiOwSkYMjXZJTIW8ZevALIPnk6d+BD74sYceMxgV+h0l24T0FCYCIvAIoNqeTFseoahlwJN7U4D8Vka/l8P0/DPx+xPaPgDfgTZvx3yKSn9z/XeB7qto56vo/Ah/KuqXDV7jA7zDJ74F3j9h+D/C7kSeIyNUi8u3k37UislFEDolIu4g8ICJ5yWNfEJH9ItIpIjtEpDG5r270rQAABFNJREFU/+vJOXgQkaXJWvp7RGSPiLSKyJdGvFexiPxWRDqSTSCfT7VGrqqtqvp74CPAF5NPcyIilyTT6hSRF0TkQ8n9pcCtwILkHUOXiCwQkRNEZFMyjwdF5KfJScFeRnL/mcB9I3aXqurTqvokEAdqROQEYJmqXj9GMo8Ay0WkIZV8OmYGLvA7TPIwUCEiq5I10wuBP0xw/mfwJvkKAnV401poct6XfwNerarlwOvx5qsfj1PwauiNwFdFZFVy/9fw1nRYjrcewjvTyNPNeNMonJDcbgHOx5uH/hLghyJynKp24032dUBVy5KvA0AC+BTe1CMnJR0/Os57HQ4MqurIH6cWETlGRI7Bm/6jA+8u4ONjJaCqA8DzeAuiOGYJLvA7TDNU6z8Lb3K3/ROc24+30EWDqvar6gPqzTmSAOYAR4lIoaruVtWdE6TzDVWNJWvFT/Ji0PtX4HJV7UgG059MNTOq2g+0AtXJ7b+p6k71uA+4Azh1gus3q+rDqjqgqruBXwJrxzl9LjC66ebDeCtcXQG8C+8O5G6gSERuF2+Jz9HpdSbTcswSCkwLOGY9vwfuB5YxqplnDL4PfB1vARKAK1T1e6r6vIh8MnlstYjcDnx6aOriMQiN+LsHKEv+vQBvAq4hRv6dEuJNrx0E2pPb5+LdSRyBV9EqwZuMb7zrjwB+AByfPLcAGG/NgA6gfOQOVX0Cb7LCoSmC/wvvzuE+4JPAAeB+EWnQFyfqKgcOTS2nDptxNX6HUVS1Ca+T91+AGyc5t1NVP6Oqy/E6MD891Javqteo6ilAA97qVP+ehs5BvMVkhlicRhoX4K2h/GhymOSfgf/EW6VpLnALIENZGuP6/8a78zlcVSvwmrNkjPMAnsOb1XvhOMd/CHxZVWPAK4DHkncRQz9OiEgBcBjenY9jluACv8MPrMdbsL17opNE5HwROSy5fkEUr4knISJHisiZyUDbC8SSx6bK9Xgds1XJYPpvqV4oItUi8g7gZ8C/q2ob3rTDc/CWGhxI1v7PHnFZM17na+WIfeXJvHWJyEq8ppoxSTYr3cUYTUEichZQpKobk7t24S02szrp1JbcfwKwO/kD7JgluMDvME6yDfyxFE49HC/QdQGbgJ+r6r14gex7eG3rIbw1Ty9LQ+WbeJ3Hu5LvcwPQN8k1T4pIF14H6fuBT6nqV8G7Q8HrVL0er1nmYrxFO0ge3463oMcLyVE8C4DPJs/rBH4F/GmS9/8lXlv+MMkfwO8Dnxix+2PAL5L5+qiqDv0wviO53zGLcPPxOxzjICIfAd6uquN1rvoCEXkQ+NjQQ1xTuG4eXtv/saramxU5hy9xgd/hSJLsDF2OdzdxOPA34Keq+iOjYg5HhnGjehyOFwngNZ0swxvlch3wc6NGDkcWcDV+h8PhmGW4zl2Hw+GYZbjA73A4HLMMF/gdDodjluECv8PhcMwyXOB3OByOWYYL/A6HwzHL+P8NI1hej59FVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(0) # Setting random seed \n",
    "sp = []\n",
    "mae = []\n",
    "for spar in [0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75]:\n",
    "    data = np.loadtxt(\"exp_data5.txt\") # Loading the data\n",
    "    X_gold = np.array(data) # The original complete data\n",
    "\n",
    "    num = 2*len(data)*spar # Number of missing elements \n",
    "\n",
    "    while num>0: # Randomly selecting 'num' elements and setting them 0\n",
    "        i = np.random.randint(0,data.shape[0])\n",
    "        j = np.random.randint(0,data.shape[1])\n",
    "        if data[i][j] > 0:\n",
    "            data[i][j] = 0\n",
    "            num-=1\n",
    "\n",
    "    X= np.array(data) # The incomplete data\n",
    "\n",
    "    K = 5 # Number of Gaussians \n",
    "    seed = 0\n",
    "    mixture, post = init(X, K, seed) #Initializing the statistics\n",
    "    mixture, post, ll = run(X, mixture, post) # Estimating the statistics\n",
    "    X_pred = fill_matrix(X, mixture) # Filling the incomplete matrix \n",
    "                                     # using the above calculated stats\n",
    "    sp.append(100*spar)\n",
    "    mae.append(np.mean(np.abs(X_pred-X_gold)))\n",
    "    \n",
    "plt.plot(sp,mae)\n",
    "plt.title('% Missing vs Error')\n",
    "plt.xlabel('Missing Data (%)')\n",
    "plt.ylabel('Mean Absolute Error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, more the missing data, greater the recovery error. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For Experimental data set 2\n",
    "\n",
    "This data set is a small subset of the Netflix version, to demonstrate the input, output of the algorithm. \n",
    "It is a matrix of 20 x 5 dimension. \n",
    "\n",
    "X_incomplete: Is the incomplete data with missing entries (0 indicated missing value) \n",
    "\n",
    "X_complete: This is the original complete data\n",
    "\n",
    "X_pred: The prediction matrix after running our algorithm and filling the matrix. \n",
    "\n",
    "RMSE Error = $\\left( \\frac{\\sum_{i=1}^{N} (X_i - Y_i)^2 }{N}\\right)^{1/2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_incomplete:\n",
      "[[2. 5. 3. 0. 0.]\n",
      " [3. 5. 0. 4. 3.]\n",
      " [2. 0. 3. 3. 1.]\n",
      " [4. 0. 4. 5. 2.]\n",
      " [3. 4. 0. 0. 4.]\n",
      " [1. 0. 4. 5. 5.]\n",
      " [2. 5. 0. 0. 1.]\n",
      " [3. 0. 5. 4. 3.]\n",
      " [0. 5. 3. 3. 3.]\n",
      " [2. 0. 0. 3. 3.]\n",
      " [3. 4. 3. 3. 3.]\n",
      " [1. 5. 3. 0. 1.]\n",
      " [4. 5. 3. 4. 3.]\n",
      " [1. 4. 0. 5. 2.]\n",
      " [1. 5. 3. 3. 5.]\n",
      " [3. 5. 3. 4. 3.]\n",
      " [3. 0. 0. 4. 2.]\n",
      " [3. 5. 3. 5. 1.]\n",
      " [2. 4. 5. 5. 0.]\n",
      " [2. 5. 4. 4. 2.]]\n",
      "Missing %: 19.0\n",
      "\n",
      "X_complete:\n",
      "[[2. 5. 3. 4. 3.]\n",
      " [3. 5. 3. 4. 3.]\n",
      " [2. 4. 3. 3. 1.]\n",
      " [4. 4. 4. 5. 2.]\n",
      " [3. 4. 4. 4. 4.]\n",
      " [1. 5. 4. 5. 5.]\n",
      " [2. 5. 4. 5. 1.]\n",
      " [3. 4. 5. 4. 3.]\n",
      " [3. 5. 3. 3. 3.]\n",
      " [2. 5. 3. 3. 3.]\n",
      " [3. 4. 3. 3. 3.]\n",
      " [1. 5. 3. 5. 1.]\n",
      " [4. 5. 3. 4. 3.]\n",
      " [1. 4. 3. 5. 2.]\n",
      " [1. 5. 3. 3. 5.]\n",
      " [3. 5. 3. 4. 3.]\n",
      " [3. 5. 4. 4. 2.]\n",
      " [3. 5. 3. 5. 1.]\n",
      " [2. 4. 5. 5. 3.]\n",
      " [2. 5. 4. 4. 2.]]\n",
      "X_pred:\n",
      "[[2.         5.         3.         3.94430072 1.5291021 ]\n",
      " [3.         5.         3.11560886 4.         3.        ]\n",
      " [2.         4.99001402 3.         3.         1.        ]\n",
      " [4.         4.21254903 4.         5.         2.        ]\n",
      " [3.         4.         3.18503123 3.64288289 4.        ]\n",
      " [1.         4.99965652 4.         5.         5.        ]\n",
      " [2.         5.         3.16507367 4.01202234 1.        ]\n",
      " [3.         4.21314223 5.         4.         3.        ]\n",
      " [2.99360835 5.         3.         3.         3.        ]\n",
      " [2.         4.6344124  3.16745754 3.         3.        ]\n",
      " [3.         4.         3.         3.         3.        ]\n",
      " [1.         5.         3.         4.00011263 1.        ]\n",
      " [4.         5.         3.         4.         3.        ]\n",
      " [1.         4.         4.50283802 5.         2.        ]\n",
      " [1.         5.         3.         3.         5.        ]\n",
      " [3.         5.         3.         4.         3.        ]\n",
      " [3.         4.40676172 4.03609658 4.         2.        ]\n",
      " [3.         5.         3.         5.         1.        ]\n",
      " [2.         4.         5.         5.         2.30793523]\n",
      " [2.         5.         4.         4.         2.        ]]\n",
      "RMSE Error: 0.3157193993270509\n"
     ]
    }
   ],
   "source": [
    "Test()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing on the NetFlix dataset \n",
    "The dataset provides movie rating for 1200 different movies by 1200 different Netflix users. \n",
    "Hence it is a 1200 x 1200 matrix. (This is a subset of the orignal Netflix Prize challenge 2009 which has rating information for millions of users on thousands of movies)\n",
    "\n",
    "The \"netflix_incomplete.txt\" dataset is the incomplete version of the this dataset (with around 23% missing entries). And \"netflix_complete.txt\" is the complete version to be used for measuring the error. \n",
    "\n",
    "RMSE Error = $\\left( \\frac{\\sum_{i=1}^{N} (X_i - Y_i)^2 }{N}\\right)^{1/2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing %: 22.793889\n",
      "RMSE Error: 0.4804908505400682\n"
     ]
    }
   ],
   "source": [
    "X = np.loadtxt(\"netflix_incomplete.txt\")\n",
    "print('Missing %:',np.round(missing(X),6))\n",
    "run_matrix_completion()   # Takes about 1-2 mins to run"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}