green_minimization_fides.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "view-in-github"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Vesalon/xla-nmf/blob/main/notebooks/green_minimization_fides.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0Jh9wPimh4OS",
        "outputId": "ee69e013-a515-4594-9cce-d416f7816eb1"
      },
      "outputs": [],
      "source": [
        "#!pip install fides"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "2v0MLj8T5txU"
      },
      "outputs": [],
      "source": [
        "import jax.numpy as jnp\n",
        "import jax\n",
        "jax.config.update(\"jax_enable_x64\", True)\n",
        "jax.config.update(\"jax_debug_nans\", True)\n",
        "import os\n",
        "import fides\n",
        "import pandas as pd\n",
        "from sklearn.cluster import KMeans\n",
        "from sklearn.metrics import silhouette_score, silhouette_samples\n",
        "import numpy as np\n",
        "import math\n",
        "import time\n",
        "import logging\n",
        "import graphviz\n",
        "from functools import partial\n",
        "import sys\n",
        "from matplotlib import pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "ZNE81o8ltkGg"
      },
      "outputs": [],
      "source": [
        "#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
        "#       set true variables for generated observations       #\n",
        "#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#\n",
        "\n",
        "### parameters ###\n",
        "\n",
        "# The amplitudes of the original sources\n",
        "\n",
        "# 2 sources \n",
        "# As = np.array([0.5, 0.7])\n",
        "# # positions of the source\n",
        "# Xs = np.array([\n",
        "#     [-0.1, -0.2],\n",
        "#     [-0.9, -0.8]\n",
        "# ])\n",
        "\n",
        "# 3 sources\n",
        "# As = np.array([0.25, 0.6,0.2])\n",
        "# # positions of the source\n",
        "# Xs = np.array([\n",
        "#     [-0.36, -0.2],\n",
        "#     [-0.08, -0.2],\n",
        "#     [0.15,-0.16]\n",
        "# ])\n",
        "\n",
        "# #4 sources\n",
        "# As = np.array([0.4, 0.3,0.2,0.1])\n",
        "# # positions of the source\n",
        "# Xs = np.array([\n",
        "#   [0.1,0.5],\n",
        "#   [-0.1,0.3],\n",
        "#   [0.1,0.0],\n",
        "#   [-0.1,-0.2]\n",
        "# ])\n",
        "\n",
        "# # 5 sources\n",
        "As = np.array([0.4, 0.3,0.2,0.5,0.7])\n",
        "# positions of the source\n",
        "Xs = np.array([\n",
        "  [-0.25,0.45],\n",
        "  [-0.08,0.2],\n",
        "  [-0.15,0.25],\n",
        "  [0.1,-0.1],\n",
        "  [0.3,0.3]\n",
        "])\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "# Advection velocity\n",
        "u_scaler = 0.05\n",
        "# U = tf_c(np.array([0.05, 0]))\n",
        "# Dispersivity\n",
        "D = np.array([0.005, 0.00125])\n",
        "\n",
        "\n",
        "### fixed known values ###\n",
        "\n",
        "# the initial time the sources begin the release\n",
        "Ts = -10.\n",
        "\n",
        "# positions of the detectors\n",
        "Xd = np.array([\n",
        "    [0, 0],\n",
        "    [-0.5, -0.5],\n",
        "    [0.5, 0.5],\n",
        "    [0.5, -0.5]\n",
        "])\n",
        "\n",
        "\n",
        "Xd = np.array([\n",
        "[ -0.2,0.6],\n",
        "[0,0.6],\n",
        "[0.2,0.6],\n",
        "[0.4,0.6],\n",
        "[0.6,0.6],\n",
        "[-0.2,0.4],\n",
        "[0,0.4],\n",
        "[0.2,0.4],\n",
        "[0.4,0.4],\n",
        "[0.6,0.4],\n",
        "[-0.2,0.2],\n",
        "[0,0.2],\n",
        "[0.2,0.2],\n",
        "[0.4,0.2],\n",
        "[0.6,0.2],\n",
        "[-0.2,0.0],\n",
        "[0,0.0],\n",
        "[0.2,0.0],\n",
        "[0.4,0.0],\n",
        "[0.6,0.0],\n",
        "[-0.2,-0.2],\n",
        "[0,-0.2],\n",
        "[0.2,-0.2],\n",
        "[0.4,-0.2],\n",
        "[0.6,-0.2],\n",
        "[-0.2,-0.4],\n",
        "[0,-0.4],\n",
        "[0.2,-0.4],\n",
        "[0.4,-0.4],\n",
        "[0.6,-0.4]\n",
        "\n",
        "])\n",
        "\n",
        "\n",
        "# time values of data\n",
        "t = np.linspace(0, 20, 80)\n",
        "\n",
        "# noise to add to \n",
        "noise_level = 1e-4"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "VN6R7ZO060_R"
      },
      "outputs": [],
      "source": [
        "# Helpers - input shaping, generating bounds and initial points\n",
        "\n",
        "def gen_flat_bounds(num_sources):\n",
        "    upper_bounds = (\n",
        "        1.5*np.ones(num_sources),\n",
        "        np.ones((num_sources, 2)),\n",
        "        np.ones(1),\n",
        "        np.ones(2),\n",
        "    )\n",
        "    # instead of 0, use epsilon of 1e-10\n",
        "    lower_bounds = (\n",
        "        1e-10*np.ones(num_sources),\n",
        "        -1*np.ones((num_sources, 2)),\n",
        "        1e-10*np.ones(1),\n",
        "        1e-10*np.ones(2),\n",
        "    )\n",
        "    upper_bounds, _ = flatten_inputs(*upper_bounds)\n",
        "    lower_bounds, _ = flatten_inputs(*lower_bounds)\n",
        "    return upper_bounds, lower_bounds\n",
        "\n",
        "def gen_init(num_sources):\n",
        "    return (\n",
        "        np.random.rand(num_sources),\n",
        "        2*(0.5 - np.random.rand(num_sources, 2)),\n",
        "        np.random.rand(1),\n",
        "        np.random.rand(2),\n",
        "    )\n",
        "\n",
        "def num_source2shapes(num_sources):\n",
        "    sample_input = gen_init(num_sources)\n",
        "    _, input_shapes = flatten_inputs(*sample_input)\n",
        "    return input_shapes\n",
        "\n",
        "def flatten_inputs(*inputs):\n",
        "    flat = []\n",
        "    shapes = []\n",
        "    for input in inputs:\n",
        "        shapes.append(input.shape)\n",
        "        flat.extend(input.flatten())\n",
        "    return np.array(flat), tuple(shapes)\n",
        "\n",
        "def unflatten_inputs(flat_inputs, shapes):\n",
        "    i = 0\n",
        "    inputs = []\n",
        "    for shape in shapes:\n",
        "        size = np.prod(shape)\n",
        "        inputs.append(np.reshape(flat_inputs[i:i+size], shape))\n",
        "        i+=size\n",
        "    return inputs\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "0QbjbBXFUMIk"
      },
      "outputs": [],
      "source": [
        "# green model functions using JAX\n",
        "\n",
        "@partial(jax.jit, static_argnames=['length'])\n",
        "def identity_vec(length):\n",
        "    # returns vector that looks like [1 0 .. 0]\n",
        "    return jnp.squeeze(jax.lax.slice(jnp.identity(length), (0, 0), (1, length)))\n",
        "\n",
        "@jax.jit\n",
        "def green(a_s, x_s, t_s, x_d, t, U, D):\n",
        "    # make matrix nt x dims\n",
        "    t_offset = t-t_s\n",
        "    xx_s = jax.lax.broadcast_in_dim(x_s, (x_s.shape[0], t.shape[0]), (0,)).T\n",
        "    xx_d = jax.lax.broadcast_in_dim(x_d, (x_d.shape[0], t.shape[0]), (0,)).T\n",
        "    tt = jax.lax.broadcast_in_dim(t_offset, (x_s.shape[0], t.shape[0]), (1,)).T\n",
        "\n",
        "    xx = xx_d - xx_s - U*tt\n",
        "    gaussians = (-xx**2)\n",
        "    gaussians = gaussians/(4*D*tt)\n",
        "    gaussians = jnp.exp(gaussians)\n",
        "\n",
        "    dividend = jnp.prod(D)\n",
        "    dividend = jnp.sqrt(dividend)\n",
        "    dividend = dividend*(t_offset)\n",
        "    dividend = 4*math.pi*dividend\n",
        "\n",
        "    res = a_s*jnp.prod(gaussians, axis=1)/dividend\n",
        "    return res\n",
        "\n",
        "green_detector_mapped = jax.vmap(\n",
        "    green,\n",
        "    in_axes=(None, None, None, 0, None, None, None)\n",
        ")\n",
        "green_matrix = jax.vmap(\n",
        "    green_detector_mapped,\n",
        "    in_axes=(0, 0, None, None, None, None, None)\n",
        ")\n",
        "\n",
        "@jax.jit\n",
        "def reconstruct_general(A_s, X_s, T_s, X_d, t, U, D):\n",
        "    U = U * identity_vec(D.shape[0])\n",
        "    return jnp.transpose(jnp.sum(green_matrix(A_s, X_s, T_s, X_d, t, U, D), axis=0))\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "tbHtiieyRRBT"
      },
      "outputs": [],
      "source": [
        "# set up partial functions to set up non-trainable values as constants for jax compilation and autodiff\n",
        "\n",
        "def make_model_fn(T_s, X_d, t): # constants\n",
        "    def model_fn(A_s, X_s, U, D): # parameters/model variables\n",
        "        return reconstruct_general(A_s, X_s, T_s, X_d, t, U, D) # all inputs into original function\n",
        "    return jax.jit(model_fn)\n",
        "\n",
        "# need to set up flattening so that residual function has 1d array as both input and output (R^m -> R^n)\n",
        "\n",
        "def make_flat_resids_fn(model_fn, observations, num_sources):\n",
        "    input_shapes = num_source2shapes(num_sources)\n",
        "    unflatten_func = jax.jit(partial(unflatten_inputs, shapes=input_shapes))\n",
        "    def resid(flat_input):\n",
        "        inputs = unflatten_func(flat_input)\n",
        "        return jnp.ravel(model_fn(*inputs) - observations)\n",
        "    return jax.jit(resid)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "cO1L-a637d4S"
      },
      "outputs": [],
      "source": [
        "# Analysis functions - to be applied to the results dataframe returned by NMF runs\n",
        "\n",
        "def filter_quantile(res, col_to_filter='fval', quantile=0.25):\n",
        "    return res[\n",
        "        res[col_to_filter] < res[col_to_filter].quantile(q=quantile)\n",
        "    ]\n",
        "\n",
        "def extract_point(sol, num_sources):\n",
        "    # extract source amplitudes and positions to use as points for clustering\n",
        "    points = []\n",
        "    input_shapes = num_source2shapes(num_sources)\n",
        "    sol = unflatten_inputs(sol, input_shapes)\n",
        "    a = sol[0] # amplitudes\n",
        "    x = sol[1] # locations\n",
        "    for p in range(len(a)):\n",
        "        point = [a[p]]\n",
        "        point.extend(x[p])\n",
        "        points.append(point)\n",
        "    return points\n",
        "\n",
        "def extract_points(res):\n",
        "    res['points'] = res['sol'].apply(lambda x: extract_point(x, res['num_sources'].iloc[0]))\n",
        "    return res\n",
        "\n",
        "def AIC(sill_avg, recon, num_resids, num_sources, sill_cutoff=0.7):\n",
        "    aic_scores = []\n",
        "    if sill_avg > sill_cutoff:\n",
        "        aic = 2*num_sources + num_resids*math.log(recon/num_resids)\n",
        "    else:\n",
        "        aic = np.inf\n",
        "    return aic\n",
        "\n",
        "def result_analysis(res, num_resids):\n",
        "    num_sources = res['num_sources'].iloc[0]\n",
        "    points = res['points'].sum()\n",
        "    if num_sources == 1:\n",
        "        centers = pd.DataFrame(points).mean().tolist()\n",
        "        avgsill = 1\n",
        "        min_Sil = 1\n",
        "    else:\n",
        "        # kmeans = KMeans(\n",
        "        #     init='random',\n",
        "        #     n_clusters = num_sources,\n",
        "        #     n_init=10,\n",
        "        #     max_iter=300\n",
        "        # )\n",
        "        # kmeans.fit(points)\n",
        "        # centers = kmeans.cluster_centers_\n",
        "        # sill = silhouette_score(points, kmeans.labels_)\n",
        "        import clusts as clst\n",
        "        from numpy.random import SeedSequence\n",
        "        from numpy.random import Generator, PCG64DXSM, SeedSequence\n",
        "        Wall = np.asarray(points).T\n",
        "        Hall = np.random.rand (Wall.shape[1], 10)  \n",
        "        seed_seq     = SeedSequence()\n",
        "        cluster_generators   = seed_seq.spawn(num_sources)\n",
        "        cluster_rand_seq = cluster_generators[num_sources-1]\n",
        "        print(\"Clustering \",num_sources)\n",
        "        centers, random_exposures, processSTE,  exposureSTE, avgSilhouetteCoefficients, clusterSilhouetteCoefficients=clst.cluster_converge_outerloop(Wall, Hall, num_sources, dist=\"cosine\" ,cluster_rand_seq=cluster_rand_seq,n_cpu=-1)\n",
        "        avgsill =  avgSilhouetteCoefficients\n",
        "        min_Sil= clusterSilhouetteCoefficients.min()\n",
        "\n",
        "    recon = res['normF'].mean()\n",
        "    aic_score = AIC(avgsill, recon, num_resids, num_sources)\n",
        "\n",
        "    return pd.DataFrame({\n",
        "        'num_source': [num_sources],\n",
        "        'aic_score': [aic_score],\n",
        "        'avg_sillhouette_score': [avgsill],\n",
        "        'min_sillhouette_score': [min_Sil],\n",
        "        'reconstruction_loss': [recon],\n",
        "        'centers': [centers]\n",
        "    })\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "78TuI3sQuJiF"
      },
      "source": [
        "# NMF run"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "104CUyVIsXb_",
        "outputId": "1a3f7524-de0e-44a1-c09b-dfb2b9e3ca8e"
      },
      "outputs": [],
      "source": [
        "# generate observations from target variables with small noise\n",
        "\n",
        "model_fn = make_model_fn(Ts, Xd, t)\n",
        "observations = model_fn(As, Xs, np.array([u_scaler]), D)\n",
        "noise = noise_level/2 - noise_level*np.random.rand(*observations.shape)\n",
        "observations += noise\n",
        "\n",
        "AA = 0 # some normalization factor to be used for AIC calculation later\n",
        "for i in range(observations.shape[1]):\n",
        "    AA += np.sum(observations[:, i]**2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "id": "OPKh5QsuCqnO"
      },
      "outputs": [],
      "source": [
        "def setup_optimizer(num_sources):\n",
        "    resid = jax.jit(make_flat_resids_fn(model_fn, observations, num_sources))\n",
        "    resid_jac = jax.jit(jax.jacfwd(resid)) # explicitly set forward mode autodiff\n",
        "\n",
        "    # return fval, grad, and (gauss_newton) hessian using resid and resid_jac functions\n",
        "    def obj(x):\n",
        "        r = resid(x)\n",
        "        jac_r = resid_jac(x)\n",
        "        loss = 0.5*jnp.sum(jnp.square(r))\n",
        "        grad = jnp.matmul(r.T, jac_r)\n",
        "        hess = jnp.matmul(jac_r.T, jac_r)\n",
        "        return float(loss), np.asarray(grad), np.asarray(hess)\n",
        "\n",
        "    ub, lb = gen_flat_bounds(num_sources)\n",
        "    # fides optimizer object\n",
        "    opt = fides.Optimizer(\n",
        "        obj,\n",
        "        ub=ub,\n",
        "        lb=lb,\n",
        "        options={\n",
        "            'maxiter': 400,\n",
        "            'xtol': 1e-6,\n",
        "            'stepback_strategy': 'truncate',\n",
        "            },\n",
        "        verbose = logging.ERROR\n",
        "    )\n",
        "\n",
        "    return opt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Lz8ChNAqOqUk",
        "outputId": "41ab35cf-1033-4e4e-c164-92752deb5aac"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 1 SOURCE, TOOK 112.95297932624817 SECONDS\n",
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 2 SOURCE, TOOK 118.17064595222473 SECONDS\n",
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 3 SOURCE, TOOK 187.03818321228027 SECONDS\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/opt/anaconda3/envs/hnmf/lib/python3.10/site-packages/fides/steps.py:194: RuntimeWarning: overflow encountered in divide\n",
            "  (self.ub[nonzero] - self.x[nonzero]) / self.s[nonzero],\n",
            "/opt/anaconda3/envs/hnmf/lib/python3.10/site-packages/fides/steps.py:195: RuntimeWarning: overflow encountered in divide\n",
            "  (self.lb[nonzero] - self.x[nonzero]) / self.s[nonzero],\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 4 SOURCE, TOOK 212.60963201522827 SECONDS\n",
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 5 SOURCE, TOOK 740.962620973587 SECONDS\n",
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 6 SOURCE, TOOK 612.1941528320312 SECONDS\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/opt/anaconda3/envs/hnmf/lib/python3.10/site-packages/fides/steps.py:194: RuntimeWarning: overflow encountered in divide\n",
            "  (self.ub[nonzero] - self.x[nonzero]) / self.s[nonzero],\n",
            "/opt/anaconda3/envs/hnmf/lib/python3.10/site-packages/fides/steps.py:195: RuntimeWarning: overflow encountered in divide\n",
            "  (self.lb[nonzero] - self.x[nonzero]) / self.s[nonzero],\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "DONE WITH RUNNING 1000 SIMULATIONS FOR 7 SOURCE, TOOK 822.4340102672577 SECONDS\n"
          ]
        }
      ],
      "source": [
        "# run NMF for up to max_num_sources\n",
        "\n",
        "nsim=1000\n",
        "max_num_sources = 7\n",
        "OUTPUTDIR = \"Results/\"\n",
        "sill_avg = []\n",
        "recon = []\n",
        "all_results = pd.DataFrame(columns=['fval', 'sol', 'grad', 'hess', 'normF', 'num_sources'])\n",
        "errors = []\n",
        "\n",
        "os.makedirs(OUTPUTDIR,exist_ok=True)\n",
        "\n",
        "for num_sources in range(1, max_num_sources+1):\n",
        "# for num_sources in range(4, 5):\n",
        "    ### define optimization object ###\n",
        "    t1 = time.time()\n",
        "    opt = setup_optimizer(num_sources)\n",
        "    ### run minimization on nsim random inits ###\n",
        "\n",
        "    results = []\n",
        "    for i in range(nsim):\n",
        "        flat_init, _ = flatten_inputs(*gen_init(num_sources))\n",
        "        try:\n",
        "            res = opt.minimize(np.asarray(flat_init))\n",
        "            results.append(res)\n",
        "        # except RuntimeError as exc:\n",
        "        #     print('hit error in fides, skipping iteration and printing error')\n",
        "        #     print(exc)\n",
        "        #     print('\\n\\n')\n",
        "        # except FloatingPointError as fl_exc:\n",
        "        #     print('hit error in jax, skipping iteration and printing error')\n",
        "        #     print(fl_exc)\n",
        "        #     print('\\n\\n')\n",
        "        except:\n",
        "            the_type, the_value, the_traceback = sys.exc_info()\n",
        "            errors.append((the_type, the_value, the_traceback))\n",
        "            print(the_type)\n",
        "    res = pd.DataFrame(columns=['fval', 'sol', 'grad', 'hess'], data=results)\n",
        "    # norm from matlab HNMF code\n",
        "    res['normF'] = np.sqrt((res['fval']/AA))*100\n",
        "    res['num_sources'] = num_sources\n",
        "    all_results = pd.concat([all_results, res])\n",
        "    all_results.to_csv(OUTPUTDIR+\"Results_optimization_sources_\"+str(num_sources)+\".csv\",sep='\\t',index=False)\n",
        "    t2 = time.time()\n",
        "    print(f'DONE WITH RUNNING {len(results)} SIMULATIONS FOR {num_sources} SOURCE, TOOK {t2-t1} SECONDS')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 242
        },
        "id": "sjoJR1gRgb5A",
        "outputId": "32b67bfd-5bd3-43f3-e9a9-98dd11e1b4a3"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Clustering  2\n",
            "Clustering  3\n",
            "Clustering  4\n",
            "Clustering  5\n",
            "Clustering  6\n",
            "Clustering  7\n"
          ]
        }
      ],
      "source": [
        "# run analysis to find which number of sources gives the best fit\n",
        "\n",
        "all_results = all_results.groupby('num_sources', group_keys=False).apply(extract_points)\n",
        "good_res = all_results.groupby('num_sources', group_keys=False).apply(filter_quantile)\n",
        "Forclusts= good_res.groupby('num_sources', group_keys=False).apply(lambda x: result_analysis(x, observations.size))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [],
      "source": [
        "def plot_results_stability_plot(numsources,avgsilsall,minsillsall,recs_all):\n",
        "    fig, ax1 = plt.subplots()\n",
        "    color = 'tab:blue'\n",
        "    ax1.set_xlabel('Number of Sources')\n",
        "    ax1.set_ylabel('Reconstruction Error [%]', color=color)  # we already handled the x-label with ax1\n",
        "    ax1.plot(numsources, recs_all, color=color,marker = 'o')\n",
        "    ax1.tick_params(axis='y', labelcolor=color)\n",
        "    ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis\n",
        "    color = 'tab:red'\n",
        "    ax2.set_ylabel('Silhouette Score', color=color)\n",
        "    ax2.plot(numsources, avgsilsall, color=color,label='Avg_sil',marker = 'o')\n",
        "    ax2.tick_params(axis='y', labelcolor=color)\n",
        "    ax2.plot(numsources, minsillsall, color='pink', label='Min_sill', marker='o')\n",
        "    ax2.legend(loc='upper right')\n",
        "    fig.tight_layout()  # otherwise the right y-label is slightly clipped\n",
        "    plt.savefig(OUTPUTDIR+\"Selection_plot.png\")\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>num_source</th>\n",
              "      <th>aic_score</th>\n",
              "      <th>avg_sillhouette_score</th>\n",
              "      <th>min_sillhouette_score</th>\n",
              "      <th>reconstruction_loss</th>\n",
              "      <th>centers</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>1</td>\n",
              "      <td>-1.323368e+04</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>9.663509</td>\n",
              "      <td>[1.4999999999910378, 0.01696364050661236, 0.25...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>2</td>\n",
              "      <td>-1.573198e+04</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>3.409504</td>\n",
              "      <td>[[1.2536621463448845, 0.7656119818886321], [-0...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>3</td>\n",
              "      <td>-1.802937e+04</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>1.307986</td>\n",
              "      <td>[[0.9233488924396693, 0.4616706819968697, 0.51...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>4</td>\n",
              "      <td>-2.839151e+04</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>0.017424</td>\n",
              "      <td>[[0.5046768621353795, 0.4114857910392516, 0.70...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>5</td>\n",
              "      <td>-3.229694e+04</td>\n",
              "      <td>0.801949</td>\n",
              "      <td>0.106190</td>\n",
              "      <td>0.003420</td>\n",
              "      <td>[[0.12473556885091912, 0.6999190290707978, 0.3...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>6</td>\n",
              "      <td>inf</td>\n",
              "      <td>0.452912</td>\n",
              "      <td>-0.311254</td>\n",
              "      <td>0.004206</td>\n",
              "      <td>[[0.158936921129094, 0.396836192837318, 0.6403...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>7</td>\n",
              "      <td>inf</td>\n",
              "      <td>0.435720</td>\n",
              "      <td>-0.066653</td>\n",
              "      <td>0.005215</td>\n",
              "      <td>[[0.5993869645099362, 0.4550893382269338, 0.13...</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   num_source     aic_score  avg_sillhouette_score  min_sillhouette_score  \\\n",
              "0           1 -1.323368e+04               1.000000               1.000000   \n",
              "0           2 -1.573198e+04               1.000000               1.000000   \n",
              "0           3 -1.802937e+04               1.000000               1.000000   \n",
              "0           4 -2.839151e+04               1.000000               1.000000   \n",
              "0           5 -3.229694e+04               0.801949               0.106190   \n",
              "0           6           inf               0.452912              -0.311254   \n",
              "0           7           inf               0.435720              -0.066653   \n",
              "\n",
              "   reconstruction_loss                                            centers  \n",
              "0             9.663509  [1.4999999999910378, 0.01696364050661236, 0.25...  \n",
              "0             3.409504  [[1.2536621463448845, 0.7656119818886321], [-0...  \n",
              "0             1.307986  [[0.9233488924396693, 0.4616706819968697, 0.51...  \n",
              "0             0.017424  [[0.5046768621353795, 0.4114857910392516, 0.70...  \n",
              "0             0.003420  [[0.12473556885091912, 0.6999190290707978, 0.3...  \n",
              "0             0.004206  [[0.158936921129094, 0.396836192837318, 0.6403...  \n",
              "0             0.005215  [[0.5993869645099362, 0.4550893382269338, 0.13...  "
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Forclusts=good_res.groupby('num_sources', group_keys=False).apply(lambda x: result_analysis(x, observations.size))\n",
        "Forclusts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHVCAYAAAB8NLYkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACa/UlEQVR4nOzdd3hTdfvH8Xf26N6L0rJBRDZSBPdWHDhA/TlxT9zi3rjFvcfjoz7gVhy4EZSyl7JXW+jebZqdnN8fhUqhQFPanqS5X9d1Ls3JycmHQpu736kBFIQQQgghRMjTqh1ACCGEEEK0DynshBBCCCG6CCnshBBCCCG6CCnshBBCCCG6CCnshBBCCCG6CCnshBBCCCG6CCnshBBCCCG6CL3aATrD0KFDKS0tVTuGEEIIIQ5ASkoKy5cvVztGUOvyhd3QoUNZtmyZ2jGEEEII0Q6GDRsmxd0+dPnCbmdL3bBhw6TVTgghhAhRKSkpLFu2TD7L90PVws7UbSDRh56FMaUX+qgEyr54FMfGBc2uiTnsfCIHn4DWHIm7eANVP7+Gp6Ig4PcqLS2lqKiovaILIYQQQgQdVSdPaIxmPGVbqPrl9Rafjz70LKJHnkHVL69T8sEt+BqqST73ETRGSycnFUIIIYQIfqoWds4tS6mZ9yGODbktPh814nRqc2fi2JCLpyKfiu+eQ2swETHgiE5OKoQQQggR/IJ2uRN9TAr6yHgcW3cZIOnz4tz2D6aMAXt9ndFoJCoqqumIjIzshLRCCCGEEOoL2skTusg4APz2mmbnfQ016GOS9/q6qVOn8uCDD3ZgMiGEEF2NXq8nLS0NrTZo2zu6NEVRqKiowG63qx0l5AVtYddEUZo/1mj2PLeLadOm8dxzzzU9TktLY/369R2VTgghRIhLTk7m0UcfxWw2qx0l7M2ZM4f33nsPZR+f82Lfgraw89mqAdBGxOFrqG46r7PG4Guo2evr3G43bre76XFUVFSHZRRCCBHaNBoNl19+OTabjWeeeQaXy6V2pLCk1+vp378/5557LgDvvvuuyolCV9AWdt7aUry2KizZQ/GUbWk8qdVjzjyY6jnvq5pNCCFE1xAbG0v//v159dVX2bBhg9pxwtrmzZsBmDhxIjNmzJBu2TZStbDTGMzo49KaHutjUjAk98DvsOGrL6d+ydfE5JyDp7oIb3URMTnn4Pe4aFj7h4qphRBCdBU7e3XKyspUTiIA1q1bB0BiYiIFBYGvWStULuyMqX1IPX9a0+P4Y64AwPb3L1R+P526hZ+j0ZuIP/4adOZIXEXrKfvkfhS3Q63IQgghuhCNRgOAz+dTOYkA8Hq9wL9/LyJwqhZ2rm1/k//kqfu8pvavj6n96+NOSiSEEEIIEbpkXrcQQgghRBcRtJMngp1ep+eqY0+nZ1IyW8rLeOOXr/H6vGrHClp6rZbJI0bTMymJLeXlvLNkAV6/X+1YQUm+VoGRr5cIBlpguMVKkl5HudfHUoedcP1XuHXrVqZPn84LL7wANK5Rd8YZZ/D111+rnCw8SGHXBo9OnMxtkyZiiotrOvfs1VfxzIyZ3DvzHRWTBacHjzuJO6+5CXN8QtO56VWVPPnaCzz48w8qJgs+8rUKjHy9RDA4NjKSqckppBkMTeeKPR6mlZXyi83Woe+dk5PDvHnz+PnnnznppJM69L1aa+TIkTQ0NKgdg+EWC5fFJzDQbCJZb+CGwu38up+/jxEWC3cmp9DbaKTM6+Xdqipm1tZ0TuB2Il2xAXp04mTuvvoqjLGxzc4bY2O4++qreHTiZHWCBakHjzuJ++9+CFNcfLPzprg47r/7IR48Ljh+EAUD+VoFRr5eIhgcGxnJ9PQMUvTN20lS9Hqmp2dwbAdva3nZZZfx0ksvMXbsWDIzMzv0vVqroqICh0P9SY5WrZb1LiePlpa26voMg4HXu2Wy1G7nrPw83qyq5O6UFI6LDK31cKXFLgB6nZ7bJk0E9pyxo9FoURSF288/jwUr5ktXEKDXaJh6/S3A3r9ed113M0vyt+IN81XG5WsVmP1/vfzcefWNPPrrj/K9KAJmaeWMTC1wT3JK4//v9hqtRoNfUbg7OYXchoZWdcs6AvzetlqtnHvuuYwcOZLU1FQuueQSHnnkEebPn88ff/zB1KlTm65NTEykqKiI448/njlz5pCamsrbb7/N0UcfTUlJCffccw+PP/54sy7UfXnggQe47LLLSElJobKyks8++4ybbroJ2LMrVi3zGhqYF0DL4cSYWIo9Hp4ob1z6ZovbzcFmM5fGx/Ozrb6jYrY7KewCcNWxpzfrft2dRqPBGBPDrDc+6MRUoUuj0WCKjWPWm/9VO0rQk69VYDQaLeaERCaPGM0bi+arHUeEEItGw9K+/drlXlqNhlSDgcWtvN/wDesDKu4mTpzI+vXr2bBhAx9++CEvvfQSjzzyCB999BG33357s8Ju4sSJlJaW8scfjevAfvDBByQmJnLkkUfi8Xh47rnnSE7e+z7suzrrrLO4+eabmTRpEqtXryY1NZXBgwe3OveBioyMbLarlMvlarbjVFsNsViYb29eCP7Z0MCEmFj0QKiMopfCLgA9k1r3j95rt+P3ejo4TfDT6g3ordb9XidfL/laBaq1X6+eSUmdkEYIdUyePJkPP/wQgNmzZxMZGckxxxzDzJkzef755xk7dix//vknAOeffz4ff/wxiqLQr18/jjvuOEaMGMHSpUsBuPzyy9m0aVOr3rd79+6UlJTwyy+/4PV62bZtG4sXL+6YP2QLdt///cEHH+Shhx464Psm6vVUNjRfz7DS68Og0RCr01ERImsdSmEXgC3lrVuZfMoLL/LKj593cJrgd9WoMbz+1P6b4q9/cGrYt6rI1yowrf16bSkv74Q0oitxKArDN6zf/4U0Ds5/M7P7fq+7clsBS1sx5iyQ1rq+ffsyatQoJkyYADQusDxz5kwuu+wyLrjgAn7++WcuuOAC/vzzT7KzsxkzZgzXXHMNAP369cPj8bBs2bKm+23evJmqqqpWvfenn37KlClT2LJlC7Nnz+b7779n1qxZnbbIc79+/SguLm563J77+yo0/zsIxXWSZfJEAN745Wtc1dUoSsujJRTFj7O6ijd+kSndAO8sWYCzqnLfX6/KCt5ZsqCTkwUf+VoFpjVfL3dVJeUb1nRyMtEVOBSlVcd8u51ijwf/Xgoyv6JQ7PEw325v1f0CMXnyZAwGA4WFhXg8HjweD9dccw0TJkwgNjaWjz76iLPPPhu9Xs/555/PP//8w6pVq4C97+rQ2t0etm/fTr9+/bjuuutwOBy8+uqrzJ07F72+c9qKbDYb9fX1TUd7dMMCVHi9JO72Z4jX6fAoCjUh0loHUtgFxOvz8syMmYBmjw+Uxscanp3xiaxnt4PX7+fJ115gX1+vJ19/UQa3I1+rQO3766UAGso++ZBHU1KZlpqGVSM/6kT78wPTyhpnXO5e3O18PK2stN3Xs9PpdFx00UXccsstDBkypOkYPHgw+fn5XHDBBXz11VeYzWZOPPFEzj///KYuW2jcj9VgMDB06NCmc7169SJuH2PId+d0Opk1axY33XQTRx55JGPGjGHQoEHt+ufsbCscDsZYI5qdOywigtVOZ8iMrwPpig3YznXqdl/HzlVTw7MzPpF17Hazcy2x3dcac1VV8eTrL8paY7uQr1Vg9vb18jrsTHv+SYqWLuL6xEROj4lhsMXCrUWFrG3HLhshAH6x2ZhSVLjHOnalXm+HrWN36qmnEhcXxzvvvENdXV2z5z777DMmT57MK6+8wtdff80jjzzCgAED+Pjjf7fmXL9+PT///DNvvvkm11xzDR6Ph2effRa73b7jF6N9u/jii9HpdCxcuBC73c6FF16I3W4nPz+/3f+sB8Kq0dDdaGx6nGEw0N9kotbno9jr5ebEJJL1eqaWNHbrzqyt4fy4OO5ISuaz2hqGmC2cFRPLbUVFav0R2kzpykd6erqiKIqSnp7ervfV6/TKdSecpbx402PKXVPfVVJOnqL6nzWYD71Wq1w1aozy5CmnK1eNGqPotVrVMwXrIV+rtn29fpj6oKLMWazMf+WdpueGWSzKrz17KWv69VdW9OmrXBAbp3peOYLryMrKUj744AMlKyvrgO6jBWWkxaqcHBWljLRYFW0HZv7mm2+Ub7/9tsXnhg4dqiiKogwdOlQ56aSTFEVRlDlz5uxxXWpqqvLdd98pDodD2bp1qzJp0iSlpKREufLKK/f7/qeffrqSm5ur1NTUKPX19cr8+fOVo48+uun5rVu3KjfddFPTY0VRlNNPP/2A/z4C/TwfabEqa/r13+N4LDVNAZTHUtOU9zO7N3vNCItF+SwrW1nRp6/yU4+eysSYWNX/jbbhUD1Ahx4dVdjtPEyZBytZd36rdLvxYwWNfADLIYdaR1pCoqLMWaz4fluoJO1SwMVotcrL6RlNP9RfTM9QYqRYlmPH0V6FXagfGRkZiqIozQq0YPv76OjP865yyMCTA+TavgafvRadJRpz5sFqxxEibBVXVrBk/Rq0Wi2njB7bdL7W7+f6okIeLy3F7fdzbFQUX2T3YKjFomJaIdR11FFHMX78eLKzs8nJyWHGjBls3bqVuXPnqh1NHCAp7A6U4se+sXGmoqVvjsphhAhvs+bPA2D8mHF7PPdhTTXnFeST73aTZjDwn8zuXBWfID8ERVgyGAw8/vjjrF69mi+//JLy8nKOPPJIvF4v559/frNZp7se//zzj9rRxX7I5Il24NiQS9TgE7D2zaH6lzdpbA0VQnS2WfPn8dClV3H8yEMxGY24dlsGYa3LxVl5eTyQmsL46BhuSkpilNXKncVFIbP4qBDt4aefftrrLNZvvvmGhQsXtvicxyMLpAc7KezagSN/JX6XHX1UIsa0PriLN6gdSYiwtHzjeraXl9ItKYUjhwznx0W5e1xjV/zcWVxMboOde1NSyImI4IvsHkwtLuYve+v3lRSiq7LZbNg6YDav6BzSC9EefB4cW5YAYO0j3bFCqOnb3MYtlMbn7Nkdu6uv6mo5Jz+PdU4niXo9b2VmcnNikvy2K4QIaVLYtZOd4+ysfUernESI8LavcXa72+p2c15BPv+rrgbgioQEPuieRbresJ9XCiFEcJLCrp04Ni9G8XowJGSiT+imdhwhwtZvy5ZgdzrpnpLK4N5993u9S1F4pKyUmwoLqfP5GGKx8EV2NsdFRnVCWiGEaF9S2LUTxe3Amb8SkO5YIdTkdLv4eUnjwO/9dcfu6mdbPRPy8ljhcBCt0/FCRgb3JadgCsVdwIUQYUsKu3Zk39g4UNsqy54IoapZua3vjt1VkdfDRQX5vFVZCcB5cXH8r3sWPXbZlkgIIYKZFHbtyL5xIYrix5TWF11UotpxhAhb3+2YQDFqwEBSd9lHtjW8wPMV5Vy+rYAKr5f+ZjOfZmVzRnRMByQVXYVWq+WIIcOYdPTxHDFkGFpt8Hy8/v777zz//POd8l6KonD66acDkJWVhaIoDB48GIAjjjgCRVGIiZHvpY4UPP/yugC/vQbX9rUAWPvIJAoh1FJSVcmitasBOCVn7H6ubtl8u50JeVvJbWjAqtXyeFoaT6alYdXIj03R3JnjjiJvxjfMmf4G/7v/MeZMf4O8Gd9w5rijOuw933vvPRRF4bXXXtvjuVdeeQVFUXjvvfcAmDBhAvfdd1+HZdlVamoqP/zwQ6e8l2iZ/IRqZ9IdK0RwaJodG8A4u91V+HxcsX0b08vL8SoK46Nj+Dw7mwEmU3vFFCHuzHFH8dnDT5KRlNTsfEZiEp89/GSHFncFBQVMmjQJs9ncdM5kMnHeeeeRn5/fdK66urrT1qUrLS3FvdvC4KJzSWHXzhwbGgs7U+bBaC3RKqcRInztHGd33IhDMRvbXoj5gTerKrl4WwHFHg9ZRiP/657F/8XGtVNSEWysZnOrjkirlRdvvA1Q0O7WkqvVakFReOGGW4m0Wlt1v0AtW7aMgoICJkyY0HRuwoQJbNu2jeXLlzed270rduvWrUydOpV33nmHuro68vPzueKKK1r1ngaDgZdeeomioiIcDgdbt27lrrvuanp+165YoQ5Zi7OdeWtLcZduxpjSC0vvUTT8/YvakYQISys3baCgtITuKakcPWwE3y/464Dut9zhYELeVh5JTePYqCjuTklhtNXKPSXF1Pr97ZRaqM1qNtMwe1673Eur1ZKZnEL993+06vqIE8dhdzoDeo/33nuPSy+9lI8//hiAyy67jHfffZcjjzxyn6+79dZbue+++3j88cc5++yzee2115g7dy7r16/f5+tuvPFGTjvtNM4991wKCgrIzMwkMzMzoMyiY0mLXQew72i1k2VPhFBX0y4UAc6O3Ztav58biwp5tLQEt9/P0VFRfJHdg+EWS7vcX4hA/fe//2Xs2LFkZWXRvXt3DjvsMD788MP9vu7777/ntddeY/PmzTz55JNUVFTstxgE6N69Oxs3buTPP/+koKCAv/76ixkzZrTDn0S0F2mx6wD2DbnEjvs/LD2GojGYUTyB/QYmhGgfs+bP49ozzubUnLFc0473/bimhuUOB8+mZ5BtNPJ+ZnderqjgrapKpO0utNmdTiJObN0vAuMOGcLsp17a73Un3nED81ataNV7B6qyspLvvvuOiy++GI1Gw3fffUfljuV69mXVqlXNHpeUlJCcnLzf173//vv8/PPPrF+/ntmzZ/Ptt9/y888/B5xbdBxpsesAnop8PNVFaPRGLD2Hqx1HiLD1+4ol2Bx2uiWlMLRPv3a991qXi7Pz8vi6thadRsNNSUm83S2TRJ2uXd9HdD6709mq4+cli9hWVop/L13xfr+fgtISfl6yqFX3a6t3332XSy65hIsvvph33323Va/xeDzNHiuK0qolWpYvX06PHj247777sFgsfPLJJ3z66adtyi06hhR2HUS6Y4VQn8vt/ncXinbqjt2VXfEztaSYu4qLsPv9jI6I4MvsHoy1RrT7e4ng4/f7uemlZ0Gj2aO48/v9oIEpLz+318KvvcyePRuj0YjRaOTHH3/s0PcCqK+v55NPPuHKK69k4sSJnH322cTFyWSiYCGFXQdx7Fj2xNJrBGilx1sItTQte9IBhd1O39TVcXZeHuucThL0et7MzOTWpCQZ6xIGvpz3O2fffyeFFeXNzm8vL+Ps++/iy3m/d3gGv9/PgAEDGDBgQIcXkVOmTGHixIn069ePPn36cM4551BcXExNTU2Hvq9oPfm500Fchevx2qrQR8ZjzjoE59ZlakcSIix9t+Av/H4/I/odRFpCIsWVFR3yPnkeN5MK8rk9KZkL4uKYHJ/ACIuV24qLKNyt20t0LV/O+52v//qDcYcMIS0+keKqCuatWtHhRdau6uvrO+V9bDYbd955J3369MHn87F48WJOPvlkFEXplPcXraN05SM9PV1RFEVJT0/v9PeOP/46JevOb5X4469T/esghxzhfOS++q6izFmsXHHqmZ3yfsdERiq5vfsoa/r1Vxb27qOcEBml+tdAjpaPrKws5YMPPlCysrJUzyLHvv8+1Pw8D6VDumI70M5dKCx9DgU06oYRIox1Rnfsrn612TgrbyvLHXaidDqez8jggZQUTBr5OSCE6FhS2HUgZ/4q/E4b+sh4TBntOyNPCNF63/w1F4Bjh4/E0knbgRV5vVxcUMCblRX4FYWJsXHMzMqil9HYKe8vRFtMnTqV+vr6Fo/vv/9e7XiiFWSMXUfye3FsXkLEwCOx9MnBVbhO7URChKV/tm4mr6SI7NR0jhk2im9z22dngf3xAtMrKlhot/NkWjp9TWY+ycrmsdJSvqir7ZQMQgTi9ddf55NPPmnxOYfD0clpRFtIi10Hs2+YD4C17xiVkwgR3jq7O3ZXuXY7E/K28ldDAxatlkfT0ngqLY2IVqwbJkRnqq6uZvPmzS0eRUVFascTrSA/VTqYY+tS/B4Xhrg0DEnZascRImztLOxOzRmLRoWxbhU+H1du38Zz5WV4FYVTo2P4PCubgabAN38X7WfnbE6dLCwdFPT6xo5EmWXbdlLYdTDF48KZtxwAa19ZrFgItfyxchn19gbSE5MY1re/KhkU4O2qKi4qKKDI46G70cjHWVlcKIu7qmbnMiGt2U5LdLz+/Ru/NysqOmZZonAgY+w6gX1DLtY+o7H2yaH2r/+pHUeIsOT2ePhx8QLOPuIYxueMY+n6taplWeF0MCFvK4+kpnFcVBRTk1MYbY3gnpJianw+1XKFo5qaGtatW8e5555LVVUVLpdL7UhhSa/X079/f84991zmzJmD3W5XO1LI0tD4S2SXlZ6eTmFhIRkZGaqND9Cao+h2w4dotDoKX5+Mt7ZUlRxChLuLTjiF/0x9kGUb1jH8ygvVjgPApNhY7kxKxqTVUuLxcHtxEUtlkHqnSkpK4rHHHsNslm5xtc2ZM4f33nuvxa7YYPg8DwXSYtcJ/M56XNv+wZw1GEvfHOoXf6V2JCHC0vc7dqEY1rc/GUnJFJaXqR2JGTU1rHA4eDY9nR5GE+9ndufVygreqKyk8/YtCG/l5eVce+21pKamylg7lSiKQkVFhbTUtQMp7DqJfUMu5qzBWPtIYSeEWipqa8hd/TeHDRrMqTljeeObL9SOBMA6l4tz8vK4JyWVM2NiuCExiVFWK3cUFVPu86odLyx4vV62b9+udgwhDphMnugk9o0LADB1G4DWGqtuGCHC2Kwda9iNz+n8ZU/2xa4o3FNSzJ3FRdj9fg61RvBldjbjIiLUjiaECCFS2HUSX30FruINaDRarH0OVTuOEGFr57InxwwfiTUIx1TNqqvjrLytrHU6idfreaNbJrcnJWFQO5gQIiRIYdeJ7Bsa94619pFlT4RQy5q8LWwpKsRsNHHs8FFqx2lRvsfDeQX5fFhdBcCl8Qn8t3sW3QxS3gkh9k0Ku05k39hY2JmzBqMxWlVOI0T4UnMXitZyKwqPl5VxQ+F2an0+DrFY+DwrmxOjotSOJoQIYlLYdSJv5XY8ldvQ6A1Yeo1QO44QYWvnOLtTc8apsgtFIH612ZiQt5WldjtROh3PpWfwYEoq5iDPLYRQhxR2naypO1Z2oRBCNXNXLqPWZiM1PoER/QaoHWe/ir1eLtlWwOuVFfgVhXNjY5mZlU1vo1HtaEKIICOFXSfbWdhZegwHnYyXEUINHq+XHxc3fi+OH3O4ymlaxwe8WFHBFdu3UeH10sdkYmZWNmfFxKgdTQgRRKSw62Tuko1468rRmqxYsoeoHUeIsBUK4+xakmu3c2beVv5ssGHRankkNY1n0tKJ0MqPcyGEFHaq2LmmnXTHCqGe7xfOx+fzMaR3XzKTU9SOE5BKn4+rtm/n2fIyvIrCydHRfJ6VzUBT8C3fIoToXFLYqcCxszu296Ggkb8CIdRQVVfL/NWrgMZJFKFGAd6pquLCgnwKPR66G418nJXFxXFxakcTQqhIqgoVOLf9g89Rh84ag6nbQWrHESJshWp37K5WOp2clbeVn+rrMWg03JmcwqsZ3YiVPU+FCEtS2KlB8ePYtAiQ7lgh1LSzsDt66AgiLBaV07Rdnd/PlKJCHiotweX3c2RkJF9mZTMihP9MQoi2kcJOJbILhRDqW1eQx6bCbZiMRo4bHvpb/c2sqWFiQT6bXS5SDAbey+zOdQmJ8oNeiDAi3+8qceYtx+92oo9JxpjSS+04QoStrtAdu6sNLhfn5ufxZW0NOo2G6xITeTezO8l6PdD4Q3+kxcrJUVGMtFjlQ0CILkavdoBwpXjdOLYuJaLfYVj65uAu3ax2JCHC0qz587j5nPM5ZfRhaDQaFEVRO9IBcygK95SUsMBu5/6UFEZZrXyZlc3MmhpOj4khbZc9Z4s9HqaVlfKLzaZiYiE6zqTYWC6LiydJr2eT280TZaUsdTj2ev2pUdFcFh9PltGIze/nzwYbT5WVUev3d2LqtpNf1lTkkF0ohFDdvFXLqbHVkxKfwKgBA9WO065m1dVxdl4ea5xO4vR6rk5MJFXf/Pf5FL2e6ekZHBsZqVJKITrOiVFRTE1O4Y2qSs7Kz2Op3c4b3TJJ07fcrjXMYmFaWhpf1NZyWt5Wbi4q5GCzhUdS0zo5edtJYaci++bFKD4vxsQs9HHpascRIix5fT5mL9qxC0UILnuyP/keDxcU5NPg9wHssTeudsfjqckp8oEgupxL4uL5vLaGz2tr2eJ280R5GcUeD5NiW14WaLDZQqHHw4c11RR6PCxzOPikpoaB5tBZI1K+j1WkuBpwFjSuoyWtdkKop6uNs9vdIWYLEdq9L3+i1WhIMxgYbrF2Yioh2iYyMpKoqKimw7iXPZMNwEFmM381NDQ7P9/ewJC9zBhf7nCQqtdzeEQEAAk6HcdHRTG3IXSGKkhhpzK7dMcKobofFs7H6/NySK8+dE9JVTtOu0vSt25Nu9ZeJ4Sa1q9fT11dXdMxderUFq+L1enRazRUen3Nzld6fSTu5d/6CqeDO4qLeTY9nZV9+zGvdx/qfT4eKy1t9z9HR5HCTmWOjQtQFD+m9P7oIhPUjiNEWKqur+Ovvxtbz7tiq135bh9sB3qdEGrq168f0dHRTce0adP2eb1C8wlRGmBvc6R6GY3cnZLMaxWVnJOfxxXbtpFhNPBACP3CJ4WdynwN1biK1gNg6RP662gJEapm5e7oju2C4+yWOuwUezz49/Jp5lcUij0eljrsnZxMiMDZbDbq6+ubDrfb3eJ1NT4vXkUhcbeJEvF6HZW+ln+JuSI+geUOB+9WV7HB5eIvewOPlJZyVmwsiSGym4sUdkFAZscKob6d4+yOHDKcyC421swPTCtr7EpqqbjTAK9WVhAaizkI0ToeYI3TyRhrRLPzY6wRrNjLcidmrRb/bt8ivh3fMxo0Lbwi+AR3YafREjvu/8i46m0yb/mc9KveJmbMJAiRL25r7RxnZ84chNYUsZ+rhRAdYcO2fDZsy8dkNHL8yK7Xev6LzcaUokJKvd5m5z2Kgkaj4cyYGFnYVHQ571dXcXZsLBOiY+hpNHJnUjJpBgMza6oBuDkxiWm7LGUyx2bj2KgoJsbG0s1gYKjFwt3JKaxyOCj3eff2NkElqL+Po0efTeSQk6j87nncFQWY0vqQcNJN+F126pd+o3a8duOtKcZdnocxKRtL71E0rP5d7UhChKVZ8//k1olZjB9zOF/M7Xrfh7/YbPxmszHcYiVJr6Pc66PM62FmVjbDLFZuTUrmyfIytWMK0W5m19cTqyvlmsREknQ6NrrdXLV9G0U7fsFJ1OubLdj9VV0tEVotF8TGcUdSMvV+Hwvtdp4tL1frjxCwoC7sTOn9cWxaiGPLEgDsdWVYBxyOMbX3Xl9jNBoxmUxNjyNDZNFN+4ZcjEnZWPvkSGEnhEpm5c7l1okXcMrow9BqtfhDZKX5QPiBxbuNpbu7pJiXM7pxcXw8K50OZtfXqxNOiA4wo6aGGTU1LT53T0nxHuc+qqnmox0teqEoqLtiXdvXYM4a3LR4ryGpB+ZuBzUVei2ZOnVqs2nQ69ev76y4B6SpO7bnMDR6036uFkJ0hL/+Xkl1fR1JsXEcOuBgteN0mt9sNt6qrATgkdRUeuxlXTAhRPAL6sKubuFnNKyZS/oVr9P9tq9Iu/QF6pZ8g33t3L2+Ztq0ac2mQffr168TE7edp2wL3poStAYz5h7D1I4jRFjy+nz8sHA+0DWXPdmXFyvKWWhvIEKrY3p6BhZN1xrLLES4COrCzjrgcCIGHknFrGcofv8mKr97nuhRZxJx8NF7fY3b7W42DdoWQhtb2zcuAGR2rBBq6uq7UOyND7i9qIgyr4c+JhMPpYbOul1CiH8FdWEXd+Sl1C74DPvauXgq8mlY/Tv1i78mZvQ5akfrEDu7Yy29R8E+tv8RQnSc2Yty8fq8HNyjF9mp4bWHc4XPx61FRXgVhVOjY5gUG6t2JCFEgIK6sNMYTKA0H7ysKH7QBHXsNnMVrsXXUIPOHIk5M3zG9wgRTGps9cxbtQIIv1Y7gKUOB8/vmAE4NTmFQSG0+bkQIsgLO8emRcSMmYil5wh00clY+uQQPfKMppatLkfxY9+0EJDuWCHUFK7dsTu9V13FT/X1GDQank/PIEYb1B8VQohdBPV3a9Uvb2Bf/xfxx19L+uWvEXfUZdhW/EDNvA/VjtZhmrpj++TQ1RZiFiJU7Czsjhg8jOiI8Fw0/N6SYvLdbtINBp5KS5efRkKEiKBex05xO6j+9S2qf31L7Sidxpm/Ar/Ljj4qAWNaH9zFG9SOJETY2VS4jXUFefTvns0JI3P4dM4vakfqdDa/n5uKCvlf9yzGRUZydUICr+1YEkUIEbyCusUuLPm8Tev0SXesEOoJ9+5YgA0uFw+XlgBwXUIih1nDs/VSiFAihV0Q2tkda+0jhZ0QatlZ2J186Bh0uvCdpf51XR2f1FSj1Wh4Ki2NNH1Qd/QIEfaksAtCji1LULweDAndMCRkqh1HiLA0f/UqqupqSYiJJeegQWrHUdXjZWWsdjqJ0+t5Lj0Dw/5fIoRQiRR2QUhxO3DkrwDAIt2xQqjC5/PxfZjuQrE7t6IwpaiQWp+PwRYLdySnqB1JCLEXUtgFKcfO7ti+Y1ROIkT4knF2/yr0eLizuAiAC+LiOCUqWuVEQoiWSGEXpOybFqL4fZhSe6OLTlI7jhBhafai+Xi8XgZk9aBXRje146hubkMDb1RWAPBQaiq9jEaVEwkhdieFXZDy22txbV8DyCQKIdRS19DA3JXLARifI612AC9VVJDb0IBVq2V6egbWLroTkBChSr4jg5h9487uWCnshFDLrFzpjt2VH7i9uIgSj4deJhMPp6aqHUkIsQsp7IKYfcMCAEzdDkJrkfEsQqhh1vy5AIw7ZCgxkZEqpwkOVT4ftxQV4VEUTo6O5oLYOLUjCSF2kMIuiPnqynCVbEKj1WHpfajacYQIS1uKClmTtwWDXs+Jo6T1fKcVTgfPlJcBcEdyMoPNZpUTCSFACrug59jY2Gon3bFCqKdpdmzO4SonCS7/ra5mdl0dBo2G59IziAvjhZyFCBZS2AW5nbtQWLKHoDFaVE4jRHjaOc7upENzwnoXipbcW1LCFpeLNIOBp9PS5UNFCJXJ92CQ81Tk46kqQqM3YukxTO04QoSl3NV/U1FbQ3x0DIcdPFjtOEHFrviZUlSI3e9nTEQE1yYkqh1JiLAmhV0IkNmxQqjL7/fz/YK/AJkd25JNbjcPlpQAcG1iIodHRKicSIjwJYVdCGjqju01ErSyAbcQavh3nJ0Udi35tr6O/1VXA/BkWjrpetlRVgg1SGEXAtxF6/HWV6I1RWDOOkTtOEKEpR8XL8Dt8dCvexZ9unVXO05QeqK8jFUOBzE6Hc+np2PQaNSOJETYkcIuJCi7zI6VvWOFUEO9vYE/Vi4DpDt2bzyKwi1FhdT4fAyyWJianKx2JCHCjhR2IaJpnF2fQ0G28BFCFU3dsVLY7VWR18sdRUX4FYVJsXGMj5bF1YXoTFIhhAhnwd/4nDZ0EXGY0vurHUeIsLSzsBs7aDCxkVEqpwlef9obeL2yEoAHU1LpYzSpnEiI8CGFXajw+3BsWgTI7Fgh1JJXUsQ/Wzej1+k56VAZFrEvr1ZW8GeDDYtWywsZGURo5eNGiM4g32khZGd3rEUKOyFUI92xreMH7iguptjjIdto5NHUVLUjCREWpLALIc6ty/B7XBhiUzEk9VA7jhBhaWdhd9KhY9DLLhT7VOPzcXNRIR5F4YSoaC6Ki1M7khBdnhR2IUTxuHBubZyVJ92xQqhj4dp/KK+pJjYyirGDhqgdJ+itcjp5sqwUgNuSkhlmka0RhehIUtiFGNmFQgh1+f1+vpNdKALycU0N39XVoddoeC49nQRp6RSiw0hhF2Icmxah+H0Yk3ugj0lRO44QYWnW/LmAFHaBeKCkmM0uF8l6A0+npSOlnRAdQwq7EON32nAW/AOARRYrFkIVPy1eiMvtpk+37vTrnqV2nJBgVxRuKirE7vczOiKCGxIT1Y4kRJckhV0IcjR1x45WOYkQ4cnmsDNnxVJA9o4NxBa3m3tLigG4MiGRoyIiVU4kRNcjhV0I2jnOzpQxAG1ErLphhAhTs3Jl2ZO2mF1fz4fVVQBMS0ujm8GgciIhuhYp7EKQr74SV9F6NBot1t7SaieEGnYue3LYwYOJj45ROU1oebqsjBUOB9E6HdPTMzBqNGpHEqLLkMIuRNk3yOxYIdRUUFrCqs0b0el0sgtFgDzALUWFVHm9HGQ2c0+yTAQTor1IYReidnbHmrMOQWO0qpxGiPDUtAuFjLMLWInXy+3FRfgVhXNiYzlTWj2FaBf61lwUd/TlAd+4dv4M/E5bwK8TreOtKsRdUYAxsTuWXiOxr/1D7UhChJ1ZufO458LLOHFUDga9Ho/Xq3akkJJrt/NyZQU3JiZxX0oKa11O1rlcascSIqS1qsUuasRpGNP6Ykzu2aojatipaE0RHZ097DmkO1YIVS1au5rSqkpiIiMZd8hQteOEpDcqK5lrs2HWapmenkGUVjqShDgQrWqxAyj/8jH89tpWXZs55ZM2BxKtZ9+YS8yYiVh6DkejN6J43WpHEiKsKIrCt7l/MvmU0xk/Zhy/LVusdqSQowB3FhfxWXYPuhuNPJaaxo1FhWrHEiJktepXo8rvp+N3NbT6ppU/voLPXtPWTKKV3CWb8NaVoTVaMGcNUTuOEGGpadkTGWfXZrV+P1MKC3H7/RwbFcVlcfFqRxIiZLWqsGv45zfwtX7siH3tHygeGSfRGewbFgDSHSuEWn5eshCn20WvjG4MyOqhdpyQtdrlZFpZGQA3JyUxwmJROZEQoemABjMYErOIHHoKUcPHY0zt3V6ZRAB2zo619B4FGhmbIkRnszud/LZsCSCLFR+ombU1fFNbi06j4bn0DBJ1sqOsEIFqcyUQOfQUUiY9irn7wZi7H0LKpMeIzjm3PbOJVnBtW43PXovOGoOp20C14wgRlpqWPZHC7oA9VFrCRpeLRL2eZ9MzkNJOiMC0urDTRSU0exw9/FSK3rmOiq+fpPzLxyj95AGiR5ze7gHFfih+HJsXAdIdK4Ravs39E4CcgwaRECPrsR0Ih6JwU2EhNp+PkVYrU5KS1I4kQtyk2Fh+6tGT5X368mlWNsP3081v0Gi4KTGRX3r2YkWfvszu0ZMJIbTOYqsLu5SJjxE1/LSmxz5HHZYew0CnR2O0YMkegq+Vs2ZF+7Kv37nsiWwvJoQatpeXsnzjenQ6HScfepjacUJensfNvSUlAEyOT+CYyEiVE4lQdWJUFFOTU3ijqpKz8vNYarfzRrdM0vR7XxTk+bR0RlsjuK+kmJO3buX24iK2uENn1YlWF3bFH9yMIaEbqRc+hyG5J1WzXyZq5Bl0v/ULMm+agbX/OCq/e74js4q9cOQtx+92oI9OlrGOQqhEumPb10+2et6vqgLg8dQ0sgwGlROJUHRJXDyf19bweW0tW9xunigvo9jjYVJsXIvXj7VGMMJq5ert28i12ynyevjb6WSF09HJyduu1evYKW4HVT+9iiljAImn3IwzbwWlH90JWi1otCgBLIci2pnPg2PLUiL6j8XaJwd3ySa1EwkRdmbNn8f9F1/OCSNHyy4U7eS58jIGmc0Mt1p5Pj2D8wvycSqK2rGEyiIjI4mKimp67HK5cLfQomYADjKbeauqstn5+fYGhuylO/aoyEhWO51cFp/AadHROBSF3231vFhRgStE/u0FPHnCVbiW4v9Mwe9qIO2SFzB3GyhFXRCw79iFwiLj7IRQxdINaymurCA6IpIjBg9TO06X4AVuLS6iwuulv9nMfSkpakcSQWD9+vXU1dU1HVOnTm3xulidHr1GQ6XX1+x8pddHor7laTndjAaGWSz0MZm4saiQJ8pKOT4qOqT+7bW+sNNoiRxyEnHHXkXkwcdQm/sJZZ8/TPSoCSSefhfaiNiOSyn2y7F5MYrPgzGxO/r4DLXjCBF2du5CAdId257KvF5uKyrCpyicGRPLWTI5Jez169eP6OjopmPatGn7vF6heUubBthb45sWDQpwR3ERfzudzG1o4MmyUs6IjsGk0bTPH6CDtbqwSzjpBqKGj0fxOIkYdCxxR1+Ot7qI0hl349y6jLT/e4bIISd1ZFaxD4rbjjN/FQDWPtJqJ4QaZJxdx1jksPNiRTkA9yanMMBkUjmRUJPNZqO+vr7paKkbFqDG58WrKCTuNlEiXq+j0udr8TXlXi9lXi82v7/p3Ba3G61GQ+o+JlwEk1YXdpY+OZR/9Tg1f/yHspn3Yuk1suk526qfKP7vrZgzZR01Ne3sjpVlT4RQxy9LF+JwOemRlsHAHr3UjtOlvF1Vxe82GyatlhfSM4jWyoLsYt88wBqnkzHWiGbnx1gjWOFoeTLEcoedJL0e6y6tc9lGIz5FoSRExs22+jvD31CDJXsYaPWYswbjd9Q3f95eS8WsZ9o9oGg9+6YFKIofU3o/dJEJ+3+BEKJdOVwufl22GJC9Y9ubAkwtLmKb2003o5FpaemERseYUNP71VWcHRvLhOgYehqN3JmUTJrBwMyaagBuTkxiWmpa0/Xf1dVR4/PxWFoavYxGhlss3JaUzBe1tV1v8kTVL68Tk3MO3W/9nPgTrqP6tzc7MpdoA39DDa7CdQBY+siadkKoQbpjO06d38+UokJcfj9HRUZyeXy82pFEkJtdX8+0slKuSUzki6xsRlitXLV9G0U7Wt8S9XrSdllKx64oXL59G1FaHZ9kZfNUWjpzbDYeLytV648QMA0QUAmqtUTjd9R1UJz2l56eTmFhIRkZGRQVFakdp8NFjzqTuKMm48hbQdnMe9WOI0TYSU9MovCz7/H7/aROOJHyHS0Dov2cFRPDI6lp+HZ8CC+029WOJDpBuH2et1XAgxRCqagLRzvH2Zm7D0JrjtrP1UKI9lZUUc7S9WvRarWcPFp2oegIn9fW8mVtDTqNhmfS0kkOkUHtQnSGVhV2qRdPR2uK2P+FO6Rc8KSM8VKJt6YEd9lWNFodlt4j9/8CIUS7m5W7oztWxtl1mEdKS1nndJKg1/NsWnrrV9sXootrVWFnTOmJIbknhqTsVh3G5J5odPJtppam2bGy7IkQqtg5zu74kYdilK2wOoRTUZhSVEi9z8dwq5VbkpLVjiREUGh19ZUy6VFo7eJ8ITJzpKuyb8glduz5mHsMQ2MwoXhcakcSIqws27COwvIyMpKSOXLIcH5avEDtSF1SgcfD3SXFvJTRjUvi41nhcPCTrX7/LxSiC2tVYVf4+uSAb+yrr9z/RaJDeMq34qkpwRCbirnHMBw7WvCEEJ3n29w/ueq0CYwfM04Kuw70q83GO1WVTI5P4LG0VDbkucjztLxgrRDhoFVdsb668oAPFP/+byw6jEO6Y4VQlYyz6zzTy8tZbLcTodXxQkYGlhDZ+kmIjiBLd3dRO8fZWXqPAm3Lmx0LITrOr0sXY3c6yUpNY1DP3mrH6dJ8wK1FhVR4vfQxmXggJVXtSEKoRgq7LspVtA5fQzU6cyTm7oPUjiNE2HG6XfyydBEgixV3hgqfj5uLCvEqCqfFxDAxJlbtSEKoQgq7rkrxY9+4EJDuWCHUsnN27GljDlc5SXhY6nAwvbwcgKnJyRxsNqucSIjOF1hhp9Fiyjw4oDXthHrsG3d0x/YZDbKrohCd7tsd4+wOPehgUuJlbc/O8G51Fb/U12PUapmenkGMVtovRPAztuO40MD+xSt+Us59GK05st0CiI7jzF+J32VHH5WAMb2v2nGECDslVZUsXrcGgFNkF4pOc09JMQVuN+kGA0+mpcuvtSIoaYCrExL4vWcvlvTpS7cda17ekJDIhJiYNt834F9l3OV56GNlYGpI8HlxbF4MgLXvGJXDCBGednbHyji7zlPv93NTUSFOv5/DIyO5KkFaS0XwuTohgTOiY3i2vBzPLuv/bnS7OPsAxogGXNjVzP0vcUddhqXXSHQRcWiMlmaHCC72DfMBGWcnhFpmzZ8LwHHDD8VkNKqcJnysd7l4uLQUgOsTEsmxWlVOJERzp0fH8GBpCd/W1+Hn38JuvctFjwP4WRHwvl/J5z4EQNJZ9zXfYUKjAUWh4OnT2xxGtD/HlqUoXjeG+HQMiVl4KvLVjiREWFmxaQPbykrJTE7h6KEj+GHhfLUjhY2v6moZarFwTmwsz6Slc1Z+HiVer9qxhAAgWa+nwL3nYtpawHAAY+4CLuxK/3d3m99MdD7F48SRtwJr71FY++ZQK4WdEJ3u29x5XHP62YwfM04Ku072WFkpA81mDjKbeT49g4sK8vGoHUoIYJPbxXCrlaK6umbnT4iKYq3T2eb7BtwV69r2zz4PEXyaFivuM1rlJEKEp53j7E6VXSg6nVtRmFJUSK3Px2CLhduTk9WOJAQAr1ZUck9yCpPj49Gg4bjIKB5KSeXK+AReraxo830DbrED0JgiiDzkeAwJ3QDwVBRgW/Uzitve5iB7o4tMIPbIS7D0HI5Gb8RbVUTlDy/gLt3c7u/VVTk2LUTx+zCl9kYXnYyvrkztSEKEld+WL6HB4SAzOYUhvfuyYtMGtSOFle0eD1OLi3m1Wzf+Ly6e5Q4HP9TXqx1LhLk5DTZuLSriyh2Te65PTGSN08m1hdvJtbe9ngq4xc6Y2puMq94ieuTp6MxR6CzRRI88g4yr3sKY0qvNQVoMZ4og9f+eAr+Xsk8fpOjta6n+/R38roZ2fZ+uzu+ow7W9cckFa19ptROis7ncbn5e2rhg+HhZrFgVcxpsvLmjFeSR1DR6yUQWoSIdcG1CApvdLi7eVsCIjRsYvnEDF24rYP4BFHXQhsIu7ugrcGxaSOFrl1H+1eOUf/kYha9PxrF5MXHHXHFAYXYXPfpsvHUVVH7/Au7iDfjqynDmr8RbU9Ku7xMOdnbHyuxYIdQhy56o76WKChY0NGDdsXixtR0XhRUiED7gsvgEtB2wymKbWuxqF3wOiv/fk4qfuoWfY0xt342uLb0PxV2ykcTT76Lb9R+SdskLRA4+Yd/5jEaioqKajshIWUwZ/t2FwtTtILSWaJXTCBF+vlvwFwAj+x9EWkKiymnCkw+4vbiIMq+HXiYTD6WmqR1JhLHchgZGdcAyPAEXdorbjj46aY/zuugkFLejXULtZIhNJWroyXiriyj95H7ql/9A3DFXEjHw6L2+ZurUqdTV1TUd69evb9dMocpXV46rZBMarQ5rn0PVjiNE2CmtqmThmsYJZqeMHqtymvBV6fNxS1ERXkXhlOhozo+NVTuSCFPzGhq4OSmJ25OSODkqiqMiIpsdbRVwYdewdh4JJ92Itf84dFGJ6KISsA44nIQTb6Bh7dw2B2mRRoOrdDM1cz/AU7YF28rZ2Fb+SNTQk/f6kmnTphEdHd109OvXr30zhTBH0+xY6Y4VQg2zcqU7Nhgsczh4prxxEtmdySkcYjarnEiEo/tTUkjQ6bg4Lp6n0tJ5KSOj6XgxI6PN9w14Vmz17+8CkHjKLbBzc2W/j/rl31P9x/ttDtISn60aT0VBs3Oeym1Y++19z0W32417lwX/oqKi2jVTKLNvzCX28AuxZA9FY7S0ewurEGLfZs2fx6OTr+HY4aMwG0043S61I4WtD6qrGWqxcEJUNM+nZ3B2fh7VPp/asUQYOXhDx/QoBlbYabSYMvpT+9fH1Pzxn8Y9YzUavNXFKN72/wHlKlyDIb5bs3OG+Ay8slxHm3gqCvBUFWKIz8DScwT2dfPUjiREWFm1eSMFpSV0T0nlmOEj+S73T7UjhbV7S0roazLRw2jiqbR0rtq+Df/+XyZEUAusK1bxk3Luw2hNESheF56KfDzleR1S1AHULf4aU3o/okefgz42DeuAI4gcfCK2Zd91yPuFg6a9Y/tKd6wQamiaHSuLFauuwe9nSmERDr+fwyIiuFYmtYhONsJi4ZWMbszu0ZMfevTk5YwMhlssB3TPgMfYucvzGlvqOoG7ZCPlXz5GxEFHkD75FWIPm0T1b2/RsGZOp7x/V9S0C0XPEaBr0/rUQogDsHOc3ak5MoEiGGx0u3iwtHEJrasTEhgXEaFyIhEuxkdH805md5x+Px9WV/NRdTUuv8K7md05Jartq1cE/MleM/e/xB11GTXzPsRdsgm/p/l+Zu09bsuxeTGOzYvb9Z7hzF28EW99JfqoBMxZQ3BuWaJ2JCHCypwVS7E57GQkJTOsb3+WbVindqSwN6uujqEWC5Ni43gyLZ2z87ZS5PWqHUt0cVfGJ/BseRkfVFc3nfuwppqLnXFcnZDAd/V1+3j13gVc2CWf+xAASWfdB4ry7xMaDSgKBU+f3qYgorMoODbmEjXsVKx9c6SwE6KTudxuflq8kAmHH8X4MeOksAsS08rKGGgyM8hi4fn0DP5vWwGeXT/jhGhnmQYDc2y2Pc7/brMxJXHPZeVaK+DCrvR/d7f5zURwsG/YUdj1PpQqzSvNF5sWQnS4WfPnNRZ2OeN46P231I4jAI+icHNREZ9nZzPIYuGupGQeKStVO5bowkq8XkZbIyiorWl2frQ1gpIDaDEOrLDT6og97Hwqf3wZb3VRm99UqMu57R98jnp0EbGYMgbg2r5a7UhChJXvFvyJ3+9neL8BpCcmUVRRrnYkARR5PdxZXMTr3TI5Ly6OFQ4HJV4vSXod5V4fSx12mTUr2s171VXcnZxMf7OJ5Y7GYWzDLBbOiI5hWlnbV/8IbPKE34chKavNbyaChN/XNG5RZscK0fnKa6pZuLZxFwqZRBFc5jY08GpFBQDT0tL4T/fuPJOewX+6d+fnnr04VrapFO1kZk0NtxYX0cdoYmpyClOTU+htNHFrcRGf7NaKF4jAd5745zciDzmuzW8ogsPO2bHWPqNVTiJEeGpa9kR2oQg6G1xOFEVBq2m+QXuKXs/09Awp7kS7+dVm48JtBYzZtJExmzZy4bYCfmth3F0gAl/vQqcnavDxmLOH4i7ZiOJpvoZd9W9vH1Ag0TmcW5fh9zjRx6ZiSO6Bp2yr2pGECCuz5s/j8Suu45hhI7GYTDhcsgtFMNDSuM2YAmh2f06jwa8oTE1O4TebTbplxQE52GxGC6xyNl9d5BCzGZ8Cq13Oll+4HwEXdsbELNwlmwEwxO2+l5nMIAoViteFc+tyrH1zsPYdQ60UdkJ0qn+2biavpIjs1HSOHX4os+a3817bok2GW6ykGQx7fV6r0ZBmMPBJ9yy2e73U+3zU+X3U+fzU7/hvS49dMsNW7Obe5BTerapiFc0LuGS9nsvjE5hUkN+m+wY+K3aGzIrtKuwbcncUdjnU/vmR2nGECDuz5s/jhgkTGT9mnBR2QSJJr2vVdQdZLBwUwH3dfj91fj91Ph/1O/5b599Z/O34f18LheGO86Gwi62WxsJYJpu0Tm+TiTUttMqtdbnoZTK2+b7tuvWA1hqD317bnrcUHcixeRGKz4sxKRt9bBremmK1IwkRVnYWdqfmjEWj0aBIq47qyr2tK6Feq6ig3OclWqsjWqclasd/o7U6onb8N1qnI0qrRafRYNRqSdRqSdS37WO3YZfWv3rf7oXhLud3LxT9fhr8HV9eHRsZydTklGatncUeD9PKSvnlAMeMdVVuRSFBp2O7x9PsfJJOj/cAfhS0+l9Y5i2fUfjaZfgdjSshJ5/zEJXfT8fX0LhistYaS7fr/iMLFIcQv9OGc9vfWLKHYu2bQ92iL9SOJERY+WPlMurtDaQlJDK87wCWrF+jdqSwt9Rhp9jjIUWv32PyBIBfUSj1enmlsqLVrVERWi3RWi1ROh3RWi3ROt1uBeCOwnDHc7uej9DqdtxDR4RWRxp77ybeG5+i7NJKuGsBuHv38W6tiTv+697PLxzHRkYyPX33oVn/TjaZUlQoxV0L5jc0cHNSMtcXbse2o/iO0mqZkpREbkNDm+/b6sJOozc27i6xgyljQOO5Zhft+U0ggpt9Qy6W7KFYpLATotO5PR5+XLyAs484hvFjxklhFwT8wLSyUqanZ+DfbWasf0eBM62sNKAuxoYdrWbFbVh0VgdNBeEerYK7FIrNCkPtv88Zd7QYxup0xOpa1828O6ffv9fC0Ob3MSk2Dg2g2a0GkMkm+/ZUeRkfZHbnl569WLujS7a/yUyl18tdxW1fK7h9d4GXboSQ49i4EI6/FnPGAHQRcU0tsEKIzjFr/rymwu6B995QO44AfrHZmFJUuEfXYqnX2+ldiz6gxuejxucDPPu7fA8mjWY/3cU7H+9ZGEZptWg1GsxaLWatlqQ2dCPvnGwy3GJlscMe8Ou7sjKvlzPztnJqdDT9TGZcip8va2v5vq6OA9mpuH0LOxFyfLZKXEXrMKX3x9JnNLYVP6gdSYiw8v2Cv/D7/Qzt049uSSlsL5dtrILBLzYbv9lsIT8ZwKUolPu8lLdh9oWGnd3Iey8MDzGbGduKdf1aOykl3DgUhU9ra4H2m5/Q+sJOUaRFrouyb8jFlN4fa98cKeyE6GQVtTXkrv6bwwYN5tScsbz+zedqRxI7+CGsW5kUwOb3Y/P7KdpLE9JIi7VVhV1rJ6WEg+4GA5FaLWt2WbtytNXKVQkJWDRafrPZeLOqss33b31hp9GQfuWbTcWdxmgm7ZIX/p3FJePrQpZ9Qy5xR16KufshaEwRKK62D9oUQgRuVu48Dhs0mPFjxklhJ0JKayebLA3jAnl3tyUls9HlairsMgwGXsnoxlKHnQ0uB1ckxONQ/Py3um1Do1pd2FV+P71NbyCCn7e6CHdFPsbELKy9RtKwZo7akYQIK7Pmz+OJK6/n6GEjsJrN2J1tW3FeiM7WEZNNurqBZjPv7NIid2pUNPluN1du3w7AepeL/4uN6/jCruGf39r0BiI0ODbkYkzMwtI3Rwo7ITrZmrwtbCkqpGd6BseNOJSv//xD7UhCtFowTTYJBXE6HaW7zJA+1Grl912+Rovtdu5ISm7z/bUHlE50GfYNuQBYegzfcxkbIUSHmzV/HgDjc8apnESIwP1is3Hcls1cXFDAbUWFXFxQwHFbNktR14Jan69phrGGxha8VU5H0/MGjWaPfYoDIYWdAMBduhlvbRlaoxlz9lC14wgRdmblNhZ2p+zYhUKIULNzssn39fUsDqIZxJNiY/mpR0+W9+nLp1nZDLdYWvW6oRYLq/r244us7HbNs9hh5+qERFL1ei6Ki0OrgUX2f8cg9jKaKPQEvrTNTlLYiSb2jY2tdta+OSonESL8zF25jFqbjdT4BEb2D2QXUiHE3pwYFcXU5BTeqKrkrPw8ltrtvNEtk7T9rMkXqdUyLTWNBfb2n0z4QkUFvYxGfu7Zi1uSknmmvBzHLquOnBYdzUJ72yebSGEnmjR1x/YeBRr5pyFEZ/J4vfy4uPF7cPwY6Y4Voj1cEhfP57U1fF5byxa3myfKyyj2eJgUG7fP1z2Yksp3dXWsdLT/RKZCj4dTt27h7Pw8jt28mZk1Nc2ef7mygtcrK9p8f/n0Fk1c29fgs9eis0RjzjxY7ThChB0ZZyfE/kVGRhIVFdV0GI0tjws3AAeZzfy1276r8+0NDNlHd+yZ0TFkGgy8egDF1f54aZz9Wu7bc4HA9S4Xtf62d2QHvvOERkvkoGMwZw1Ga43dYyxI6Yx72hxGqEzx49i0kMhDjsfSNwdnwSq1EwkRVr5fOB+fz8fg3n3pnpJKQWmJ2pGECDrr169v9vjBBx/koYce2uO6WJ0evUZD5W6LI1d6fSRGtLwTRpbBwM1JSVxYkE+oLqkccItd3DFXEnfMlaDR4qnIx122tdkhQtvO7lhrn9EqJxEi/FTV1TJ/deMvVKdKq50QLerXrx/R0dFNx7Rp0/Z5vULzXbM0tLyRlhZ4Ki2dVyoqyD+AyQtqC7jFLmLAOMq/fhLnliUdkUeozJG3Ar/Ljj46CWNqH9wlG9WOJERYmTV/HuMOGcr4MeN49atP1Y4jRNCx2WzU19fv97oanxevopC420SJeL2OSt+e7XERWi2DLBYGmM3ck5ICNBZ7Wo2GVX37ccX2bQc0qaGzBNxip/i8eKuLOiKLCAY+D46tSwGZHSuEGnaOsztqyHAiLVaV0wgRujzAGqeTMdaIZufHWCNY4XDscb3N7+e0rVuYkLe16ZhZU8MWl4sJeVtZ1cJrglHAhV3d4q+IHnFaR2QRQaKpO1YKOyE63bqCPDYVbsNkNHLciEPVjiNESHu/uoqzY2OZEB1DT6ORO5OSSTMYmFnTuF3XzYlJTEtNA0ABNrndzY4qnw+3orDJ7W62JEl7GW6x8GRaGh93zyJ5R8vi+OhohrVyrb2WBNwVa+52EObugzD3HIGnogB2m9FR/tXjbQ4jgoNj8xIUnwdDQib6hG54K7erHUmIsDJr/jxuPud8xo8Zx5fzflc7jhAha3Z9PbG6Uq5JTCRJp2Oj281V27dRtGNLr0S9vtk2aJ3puMgonkhL49u6OgaYTBh3TEaN0Gq5Mj6Bqwvb9tkbcGHnd9maFrIVXZPituPMX4ml5wisfXKoq5RxPkJ0pp2F3SmjD0Or1eI/gKUPhAh3M2pqmLHbWnE73VNSvM/XvlJZwSsdtOzJ1QkJPFRawjd1dZwcHdV0foXDwbUJiW2+b8CFXeX3L7T5zUTosG/IbSzs+uZQt0AKOyE607xVy6mx1ZMcF8+o/gNZsOZvtSMJIdpZttHIEnvLY/2itG1fZrjNr9RaojFlHIQpYwBaS3SbA4jgZN+4EEXxY0rriy6q7b85CCEC5/X5+GHhfEB2oRCiqyr3esky7tkNPMxiYXtn7hWrMZhIOOkmul3/X1IueIKUC56k23UfkHDSjWj0pjYHEcHFb6/BtX0tIGvaCaGGpl0opLATokv6pLaGqckpHGI2oyiQrNdzalQ0tycl878dkzvaIuCu2LijL8eUeTBlnz+Ma/saAEzdBhJ/7JXEHT2Zqp9ebXMYEVzsG3MxZw7E2jeH+mXfqh1HiLDyw8L5eH1eBvXsTVZqGvn7GQskhAgt71ZVEaXV8l5md0waDR9kdsetKLxXVcXHexkT2BoBt9hZ+46h8ocXcW5ZiuJ2oLgdOLcsoXL2S1j7HdbmICL4OHYse2LKPBitOWo/Vwsh2lONrZ4//14JyN6xQnRVL1RUcNimjUzMz+e8gnzGbtrISwc4WaNNXbE++55NhH57rXTFdjHe2lLcpVvQaHVYeo9SO44QYUe6Y4Xouh5NTcWq0eJUFFa7nPztdGJXFCwaDY+mprb5vgEXdq7C9cSOvQB0/w740+iNxBx2Hq6idW0OIoLTzqVtZLFiITrfzsLuyCHDidpt9XwhRGg7PToGs1azx3mTVstp0TFtvm/AY+yqf32T5HMfotu17+MuywMUjMk9UHweyj65v81BRHCyb8glduwFmLOHojGYUDwutSMJETY2bi9gfUE+/bpncfzIQ/n8j9/UjiSEOEARWi0aQLPj/1277GihAw6PiKBqt80fAhFwYeepyKfozSuJGHgUhvhuoNFgXzuXhjVzULzuNgcRwclTnoenuhhDXBqWHsOxb5ivdiQhwsqs3Hn0657F+JxxUtgJ0QUs7N0HhcYtzL7v0XOP5xXglYq2j7MLuLADULxubCt/bPObitBi35hLzKgJWPrmSGEnRCebNX8et038P07JGSu7UAjRBVyyrQANGt7LzOSmokJqfb6m5zyKQpHHS3lHt9hZeo/CsWUp+H37HUTv2LSozWFEcHJsaCzsrL1GUqnVg7/t/+CEEIH565+VVNXVkhgTS87AQfy1Y6asECI0LXE07jZx/JbNTXvW7i5Nr6d4L8/tT6sKu6QJ97L95Qvx22tJmnDv3i9UFAqePr1NQUTwchWuw2erRhcZhznrEJxbl6kdSYiw4fP5+GFhLhccdyLjc8ZJYSdEF/Fjz14csXkTVbu02AHEaLX81LMXgzasb9N9W1XYFTx1Wov/L8KFgn3TAqKGnIS1T44UdkJ0slm5cxsLuzHjuOvNl9WOI4RoB3vOh220+4SKQAU8xi5i4NE0rJsLu/f/avVEDDichtUyuLcrsm/IJWrISVj6HAo/vUrj8E4hRGeYvSgXj9fLQdk96ZmewZaiQrUjCSHa6I6kZKDxU/SGxEQc/l1mxWrgELOFdS5nm+8f8Dp2CSffhNa053pKWqOFhJNvanMQEdyc+avwO23oI+MxZfRTO44QYaXWZmPequUAjB9zuMpphBAHYoDZxACzCQ3Qx2RqejzAbKKH0ch6l5N7itu+hWDgs2I1GmihiVAXlYjfZW9zEBHk/F4cm5cQMfBILH1ycBXKYtRCdKZZ8+dx9LCRjM8Zxwuf/U/tOEKINrp02zYAHktN5fGyMhraeaZ7qwu7tEteQFEUUBRSJj2G4v93sJ9Gq0Mfk4Jj69J2DSeCi31jLhEDj8TaN4eaOe+pHUeIsDIrdx7PX38Lhw8eSnREBHUNDWpHEkIcgHtKSgDobjCQaTCyxGE/oLF1O7W6sLNvXACAMaUnjq3LUDz/9v8qPi/e2lLs62WNs67MsWUpiteNIS4dQ2IWnop8tSMJETY2F25nbf5WBmT14MRRY/jk95/VjiSEOAAxWi3Pp2cwympFAU7auoXtHg+PpKRS7/fzVHlZm+7b6sKu9q/Gpn9vbSkNa1uYPCG6PMXjxLF1OdY+h2Ltm0OtFHZCdKpZ8+cxIKsH48eMk8JOiBB3V3IKHhSO2bKZb3v0aDr/Q30ddyWn8FR52+4b8OQJT8U2jMl7boFhTOuLMbV321KIkGHfmAuAtW+OykmECD+z5s8D4ORDx6DT6VROI4Q4EIdFRPBceTmluy1EnO/xkG4wtPm+ARd28cddgz46aY/zuqgE4o+7ps1BRGhwbFqE4vdhTOmFPiZF7ThChJXcNX9TWVtDfHQMYwYeonYcIcQBsGg1OFuYOBGn0+E+gLF2ARd2hsRM3CWb9jjvKd2CISGzzUFEaPA76nBtWw2Apc9oldMIEV58Ph/fL2wcyzx+zDiV0wghDsQSu4PTomOaHis0Llp8WVw8C+1tnxwVcGGn+DxoI+L2OK+LjENRZHPqcPBvd+wYlZMIEX52dsdKYSdEaHumvIxzY2N5I6MbBjTclpTEN9k9GGG18lx5GwfY0YbCzpm3grgjLkZjtDad05giiD38Ypx5y9scRIQO+4bGGdKmbgPQWmPVDSNEmPlxceMuFP27Z9M7Q3pJhAhVm91uzszbyt9OJ/PtdiwaLT/b6pmQt5VtHk+b7xvwAsXVv71DyvlP0O2ad3GXbgHAkNITf0M1Fd8+2+YgInT46stxFW/AlNYXa59Dsa38Ue1IQoSNuoYG/li5jGOHj2L8mHE8/+nHakcSQrRRhc/Hy5UV7XrPgAs7n62S4veuJ+KgIzEm98DvdWP7+xca1v4BuyxaLLo2+4bcHYVdjhR2QnSyWfPnSWEnRIgbbrHs8/mlDkeb7hv4lmKA4nHJh3mYs2/MJe6IizFnDUZjtKC42/YPUAgRuFnz5/HCDbcy7pAhxEZGUWOrVzuSECJA/8nsvse5XefCDtqwvk33Dbiwixh49D6fb1j9W5uCiNDirdyOp3IbhoRMLL1GYl87V+1IQoSNrcWFrN66hYE9enLiqBxm/PaT2pGEEAEavWljs8d6jYYBJhM3JibxQkXbJ08EXNjFHXtls8carQ6NwQQ+L36PSwq7MGLfsICYnEysfUZLYSdEJ5uVO4+BPXoyfsw4KeyECEG2Ftawy7XbcZeXcVdyCufk57XpvgEXdttfmLTnTeLSiT/+WuoWfdGmECI02TfmEpNzDpaeI0BnAF/bZ/EIIQIza/5c7jr/Yk46dAx6nQ6vT8Y4C9EVVPt89DAa2/z6gJc7aYm3uoiaP94n/pgr2uN2IkS4izfira9Aa7JiyRqsdhwhwsqCNf9QXlNNXFQ0hw2S7z8hQk1fk6nZ0c9kYqw1gvtTUlnncrb5vm2aPNESxe9HF5nQXrcTIUHBvmEB0cNPxdI3B8eWJWoHEiJs+P1+vl/wFxefeCrjcw7njxXL1I4khAjAF1nZTbtN7Gql08G9JSVtvm/AhZ2l96jdzmjQRcYTNexUXIVr2hxEhCbHxlyih5+Ktc9oqn58BWT3ESE6zazceY2F3Zix3PbadLXjCCECcPyWzc0e+4Eqn++A9omFNhR2SRPubX5CUfDb63AWrKT6t3cOKIwIPc5t/+Bz1KOzxmDqdhCubf+oHUmIsPHT4oW4PR76ZmbRNzOLDdvy1Y4khGilIq+3Q+4bcGFX8NRpHZFDhCq/D8emRUQOOgZr3xwp7IToRPX2BuasWMrxI0czfsw4np0phZ0QoWSExcKl8Qn0MhpRgM1uF+9VVbV5cWIIdPKEVkf6VW9jSJD9CcW/7BvmA2Dtk6NyEiHCz6z58wAYnzNO5SRCiECMj47mnczuOP1+Pqyu5qPqalx+hXczu3NKVHSb7xtYYef3odEZUDiw/l/RtTjzluN3O9HHJGNM6aV2HCHCyqzcxsLusEGHEHcAHwZCiM51ZXwCz5aXcWtxER/WVPNhTTW3FhfxXHkZVye0fTJqwMud1C+dRcyhZ4OmXVZKCUj06HPIuvNb4mRZlaCieN04ti4FwNJXWu2E6Ez5JcX8vWUTep2ekw4do3YcIUQrZRoMzLHZ9jj/u81GN4OhzfcNeIydKb0f5qzBmLOH4qnIR3E3X2ul/KvH2xxmX4ypfYgafALusq0dcn9xYBwbconodxjWPjnUzvtQ7ThChJVZ8+cxqGdvxo8Zx8e/zFY7jhCiFUq8XkZbIyiorWl2frQ1gpIDmFgRcGHnd9mwb/irzW/YFhqDmcTxt1E5+yVixuy588WujEYjJpOp6XFkZGRHxxOAY/NiFJ8XY1IW+rh0vNVFakcSImzMmj+Pu//vUk4clSO7UAgRIt6rruLu5GT6m00s3zFZYpjFwhnRMUwrK2vzfQMu7Cq/f6HNb9ZW8cddg2PzYpz5K/db2E2dOpUHH3ywc4KJJn5XA86CVVh6DMPaJ4e6RZ+rHUmIsLFo3WrKqqtIjotn3CFD+X25LBYuRLCbWVNDhdfLJXHxnLhjfOxml4tbi4v4rYUu2tYKeKBcyqTH0Jgi9jivMVpImfRYm4PsjXXA4RhTe1H9x39adf20adOIjo5uOvr169fumUTL7BtyAbDKODshOpXf7+e7BY09KePHyOxYIULFrzYbF24rYMymjYzZtJELtxUcUFEHbSjsTN0HodHt2dCn0RsxdRt4QGF2p4tKJP6YK6iY9WyrN5h3u93U19c3HbYD/AKJ1nNsWgiAKaM/ush4ldMIEV6alj2Rwk6IkGIAUvR60nY72qrVrzQkZf/7/wnd8UfU//ukRoul53B8tso2B2mJMbU3uog40i6Zvstb6TBlDiRq2KkUPHOmbGEVRHy2KlyF6zBl9MfS+1BsK35QO5IQYeOnJQtwud30zsikf/ds1hXkqR1JCLEPWQYDj6amMcRiaXZeAyjAoA3r23TfVhd2aZe+CIoCikLKeXt2uSoeN1W/vNGmEHvjzF9J0TvXNTuXcPJNeCq3U7fwcynqgpB9Qy6mjP5Y+42Rwk6ITtTgcPD7iqWcOCqH8WPGSWEnRJB7LC0NnwLXFm6n3OtttxWCW13YFb4+GdCQcfXblHxwCz577b9P+ryNj9u50FLcDjwVzbfIUTwu/M76Pc6L4GDfmEvcUZdizhyE1hSB39WgdiQhwsas+fOaCrunZ/xX7ThCBIVJsbFcFhdPkl7PJrebJ8pK97pl17GRkUyKjaO/yYRRo2GT280rFRX8ZW//z7L+JjPn5Oex1e1u1/u2eoydr64cX10ZBU+dhrtk047HO46Gamk9EwB4q4twl+ej0emx9B6ldhwhwsq3O3ahGDPwEOKjY1ROI4T6ToyKYmpyCm9UVXJWfh5L7Xbe6Ja51zFsIyxW5tsbuLpwO+fk57HI3sCr3boxYJdl1NrLZpeLOJ2u3e8b8OSJiIOPxtJzRNPj2CMvJfOmGaT839PoopPaNVxLSv83lepf3+rw9xFtJ3vHCqGOgtISVm7agE6n42TZhUIILomL5/PaGj6vrWWL280T5WUUezxMio1r8fonyst4t6qKf5xO8j0epldUkO92c2Q7rYkbodU2Hc9VlHNrUhIjLVZidjm/82irgF8ZM/pcFG9js6ExvT9Rw06hes57+O11xB8tW32Jf5c9Mfcchkbf/r/lCCH2bufesTI7VnRVkZGRREVFNR1Go7HF6wzAQWYzfzU070adb2/YY8LC3mhoLMZq22nR74W9+7Bgx/F2t0wGmy28m5nJn7uc33lNWwU8n1YXnYinuhgAa5/R2NfPx7byR1yFa0k5b1qbg4iuw1O2BW9tKfqYFMw9huLYuEDtSEKEjVnz53HvhZM5cVQOBr0ezwFsTSREMFq/vvls0QcffJCHHnpoj+tidXr0Gg2V3uZFWaXXR2JE67pAL42Lx6LVMru+fv8Xt8Il2wra5T77EnBhp7idaC1R+OrLsfQYSt3irxvPe91o9C1XzSL82DfkEj3yDKx9cqSwE6ITLV63hpKqSlLjEzh88DB+XbpI7UhCtKt+/fpRXFzc9Njlcu3zemW3+aYaGhf52J+To6K4NjGRGwq3U9VOLXZL9jJpoz0FXNg581aQcNKNuEs3o4/LwLF5MQCGxO54a9u+t5noWnYWdpbeo0CrA7/sXSlEZ1AUhe9y/2TyKaczPmecFHaiy7HZbNS3ogWtxufFqygk7jZRIl6vo3I/hdqJUVE8kprGzUWF5NrtB5R3V30DmISxYT8F694Evlfsz68RO+5C9NGJlH/1OH5n4xfXmNqbhrV/tCmE6HpchWvxNdSgi4jFnHkwzvyVakcSImzMyp3XWNiNGceUl59VO44QqvAAa5xOxlgj+HWXXajGWCP2uW3XyVFRPJqaxu3FRcxtaN9lTr7IykahsdVwXzplgeKmN3M1UP3L63ucr/3z4zYFEF2U4se+aSFRg0/A0idHCjshOtHPSxbidLvomZ7BQdk9WZO3Re1IQqji/eoqnkxLZ7XTyQqng3NiYkkzGJhZUw3AzYlJJOv1TC1p7No9OSqKaWnpTCsrZaXDQeKO5UicioLNf+DLuh2/ZfMB32N/2rQZmcYUgSmtLzprDGh2nVir0LD693aKJkKdY0MuUYNPwNp3NNW/vAHttq62EGJf7E4nvy1bwsmjD2P8mHFS2ImwNbu+nlhdKdckJpKk07HR7eaq7dso2jGpKFGvJ81gaLr+3Ng4DBoN96ekcn9KatP5L2truaekeI/7B6qoEyYzBVzYWXqNInH8bWgMJhS3o/lHtSKFnfiXI38lfpcdfVQixrQ+uIs3qB1JiLAxa/68xsIuZxxPfvwfteMIoZoZNTXMqKlp8bndi7WOnrV6VEQk8xpseHf8/7783rD37uJ9Cbiwizt6Mra/f6bmjw9QvG0b2CfChM+DY8sSIgYcjrVvjhR2QnSib3P/5DUgZ+AgEmNiqaitUTuSEGHvpYwMDt+8iSqfj5cyMvZ6XaeOsdNFJlC/5Bsp6kSr2DcuaCzs+uRQ84e0GgjRWbaXl7J843qG9unHyaMP44Mfv1M7khBh7+BdirWD21i47U/AO084ty7DmNb2FZFFeHFsXozi9WBI6EbEwUdjHXA4psxBu43NFEJ0hFnzZRcKIcJNwC12ji2LiTvyMgwJ3fGU56H4mw8EdGySNZPEvxS3A3dFPqbU3iSeckvTeW9dOVW/voljx/ZjQoj2N2v+PO6/+HJOGDkao8GA2+NRO5IQYe8Qs5kYnY55uyylclp0NNcnJGLRavnVZuOxslI8rVlFuQUBF3bxJ94AQMxhk/Z8UlEoePr0NgURXZOlbw7GlF57nNdFJZB0xt2Uf/W4FHdCdJClG9ZSXFlBWkIiRwwexs9LFqodSYiwd11CIosc9qbCro/RxCOpaXxVW8sWt4vL4hMo93p5pbKiTfcPuLAreOq0Nr2RCEMaLfHHXNnyUxotiuIn/pgrKdy4EJQDXx9ICNGcoih8m/snV5x6BuPHjJPCTogg0N9s5sWKf4u2k6OjWOVw8EBpCQAlHi/XJya2ubCTgU6iw5i6DUQfnYRG0/Ia2xqNFn10EqZuAzs5mRDhQ8bZCRFcorVaKn3/DmMbYbHy5y7dsn87naQa2rTMMNDGBYpNmQcTPWoChoRuoICncht1i77AtX11m4OIrkcXGdeu1wkhAvfL0oU4XE6yU9M5uEcv/tna8SvfCyH2rtLnI8NgoMTrxQAcZDY3a52L0GrxHsB6/gG32EUcdCQpEx9F8TipXzqL+mXfonhdpEx6FOuAI9qeRHQ5Plt1u14nhAicw+Xi12WLAWm1EyIY/Nlg45akJIZbLNyclIzT72ep3d70fD+TiQK3u833D7iwixkzkeo571HxzVONhd3Sb6j45imq//gPsS1NqBBhy7V9Nd66cpS9jJ9TFAVvXbm09ArRwaQ7Vojg8UJFBX4F/pPZnbNjYri/tIRd56tPiIlhvr1hr6/fn4C7YvUxqS0uaeLYuJC4wy9qcxDRBSl+qn59k6Qz7kZR/GhaWLuufuWPMnFCiA72be6fABw64GCS4+Ipq65SOZEQ4ava5+PCbQVEarXY/X52/wS8uagQu7/tn4sBt9h568sxZw3e47w5azDeurbN4BBdl2NDLuVfPY6vvrLZeb/HhUajIWbUBIxpfVVKJ0R4KKooZ+n6tWi1Wk4ZfZjacYQQgK2Fog6g1u/nQFacDLjFrm7Rl8QfexXGlJ64CteComDqNpDIQcdQ9cubBxBFdFWODbkUblyIqdtAdJFx+GzVuIo3kHzWfViyh5B8zkOUfnwXnop8taMK0WXNyp3H8H4DGD9mHO/9MEvtOEKIDqKhca/ZgFj65BA96gwMCZnAjlmxC7/AsSn41khKT0+nsLCQjIwMioqK1I4jdqExmEme+AjmjAH4bNWUfHwn3mr5OxKiIwzr25+lb/6XBoeDhNOPxXUAg7OFUIN8nrdOm5Y7cWzMxbFRdgsQB0bxOCn/9EFSzpuGMaUnKRMfpeSjO/DVS5e+EO1t2YZ1FJaXkZGUzFFDRjB70Xy1IwkhOkDAY+yMqX1aHBNlTOuLMbV3u4QS4cPvaqD0k/vwVG5HH5NMysRH0Vpj1I4lRJe0cxLFVadNYNLRx3PEkGFotbJOvRBdScDf0fHHXYM+OmmP87qoBOKPu6ZdQonw4rfXUjrzXry1ZRgSupFy7iNoTRFqxxKiy6moqwHgjLFH8L/7H2PO9DfIm/ENZ447St1gQoh2E3BhZ0jMxF2yaY/zntItTWPuhAiUr76C0pn34rNVY0zpSdI5D6IxmNWOJUSXcea4o5h6waUoSvNh1RmJSXz28JNS3AnRRQRc2Ck+D9qIPbeA0kXG7XUhWiFaw1td1FjcOeoxZwwgacK9oDOoHUuIkKfVannhhlsBZY+9m7VaLSgK06+/RbplhegCAv4uduatIO6Ii9EYrU3nNKYIYg+/GGfe8nYNJ8KPpyKfsk8fwO92YMkeQtJpd4BWp3YsIULauEOGkJmcgraFRcKhsbjrnpLKuEOGdG4wIUS7C7iwq/7tHXRRiXS75l1SJj1OyqTHybj6HXQRsVT/9k5HZBRhxl28gbLPH0bxurH2zSHhpJtoXJlHCNEWafGJ7XqdECJ4Bbzcic9WSfF71xNx0JEYk3vg97qx/f0LDWv/AL+vIzKKMOQq+Jvyr6aRdOY9RB58NIrbQdXPr6kdS4iQVFzVuiWEWnudECJ4tWkdO8Xjwrbyx/bOIkQzjs2LqfjuORLH30bUsFPwu+3U/PEftWMJEXLmrVrBtrJSMhKTWhxH5/f72V5exrxVKzo/nBCiXbVppGzEwKNIueBJMq79D7odS59EjTgdS+9D2zWcEPa1c6n68RUAYkafQ/Toc1ROJETo8fv93PTSs6DR4G9hc3GNRsOUl59r8TkhRGgJuLCLHHIScUdfjmPLUrTmSDQ7BuP6nTaiR5ze7gGFsK38kaod4zfjjriYyKGnqJxIiNDz5bzfOfv+OymsKN/jObfXw99b9lzGSggRegIu7KKGj6dy9kvU5X7SbEydu2QThqSsdg0nxE71i7+k5q//AZBw/DVEDDxa5URChJ4v5/1O9qTTOHLKVZz38D0cOeUqflq8AJPByJu33a12PCFEOwh4jJ0+JgV36eY9zis+jywoKzpU7Z8foTVFED3iNBJOvgm/2yF7FgsRIL/fzx8rljU9zi8p4Z/3ZnDU0BFcfsoZvP3dV+qFE0IcsIBb7Ly1pRiTe+5x3tJzOJ7Kbe0SSoi9qf71LWyrfkaj1ZF02h2Ys4eoHUmIkJZXUsS97zTOOH/6mhtJS5AlT4QIZQEXdnWLviD++Guw9h8HGg3G9H5E55xL7OEXUbfoi47IKMQuFCpnv0TDuj/R6A0knXkvpowBaocSIqS9+MVMFq75h9jIKF6+6Q614wghDkDAhV3D379Q+9f/iDvyUjQGU+NSFENOovrXN7GvndsRGYVoTvFTMesZHFuWoDWaST7nQQwttCILIVrH7/dz+dOP4vF6mXD4UUw4XPaNFSJUaQBlv1fthdYS3Th93l4LgC4yAZ+tsr2ytYv09HQKCwvJyMigqKhI7TiiHWn0JpLPfQhz5sH47LWUfHQn3qrtascSImQ9dOlV3H/x5RRXVnDQxedSY6tXO5IQTeTzvHUOaMdnv6MOv70WbUQsccdeRfqVb7ZXLiH2S/G6KPvsIVzFG9FZY0iZ9Ci66GS1YwkRsh778F3W5G0hLSGRZ665Se04Qog2aHVhpzFFkHjqbXS74SMyrv0PUcPHAxpixl5AxlVvY0rvR+UPL3RgVCH2pLgdlH36AO6KfPRRiY3FXUSc2rGECEluj4fLn34Mv9/P5FNO5+hhI9WOJIQIUKsLu7jDL8KUORDbP7/id9qIO/pyks9+AHO3gyj79EFKPrhFxtgJVfgddZTNvA9PdTGGuHSSJz6K1hyldiwhQlLu6lW8+vVnALx5691YTCaVEwkhAtHqws7SaySV30+n5vd3Kf/8YdBo8FQXUjrjHlzb/unIjELsl89WRdnMe/HWV2JMyiL5nIfQGC1qxxIiJE198xUKSkvoldGNhy+7Wu04QogAtLqw00XGN61T560tRfG6sa38qcOCCREob20pZTPvxWevxZTel+Sz7kejN6odS4iQY3PYufq5aQDcfPZ5jOh3kMqJhBCt1frJExoNis/772O/H8Xj7IBIQrSdp3IbZZ/cj9/VgLn7IBLPmAragDdYESLs/bBwPh/9PBudTsfbt9+DXqdTO5IQohVa/4mn0ZBwys3g9TQ+1BuJP+E6FHfz4q78q8fbNaAQgXKXbqbs04dInvgw1l4jSTz1FipmPQOKX+1oQoSUKS8/ywmjRjO4d1/uOO8iHv/wPbUjCSH2o9Utdg3//IrfXoPf3YDf3UDDmt/x2SqbHu88hAgGrsI1lH/5OIrPQ8SAw4k/4Tq1IwkRcipqa7jppWcBuP+iy+nXPUvlREKI/TmgBYpDgSxoGN6sfceQePqdaLQ66hZ/RfVvb6sdSYiQ890T0zl59GHMW7WcI266CkXp0h8bIkjJ53nrHNACxUIEO/uG+U3rK0aPPIOYw85TOZEQoefq56ZRb29g3CFDuWr8BLXjCCH2QQo70eU1/PMbVT+/DkDs2AuIGnG6yomECC3bykqZ+tYrADx51fVkJMkOL0IEKynsRFioX/Yt1XM/ACD+mCuIPOQ4lRMJEVpe/eoz5v+ziuiISF67+S614wgh9kIKOxE26nI/oXbh5wDEn3gD1v5jVU4kROhQFIXLn34Ul9vN+DHjmHi0/HIkRDCSwk6ElZo571G//Ac0Gi2Jp96GuecItSMJETLW5m/lsR1Lnrx4w23ER8eonEgIsTsp7ETYqfr5NRpWz0Gj05N0xlRMmQerHUmIkPHEx+/z95ZNJMfF8/x1N6sdRwixGynsRPhR/FR8/zz2jQvRGkwkn3U/xtQ+aqcSIiR4vF4uf/pR/H4/F51wCieMylE7khD7NCk2lp969GR5n758mpXNcMu+9xEfYbHwaVY2y/v05ccePZkYE9s5QduJFHYiPPl9lH/9BI68lWhNVpLPfQhDoiy+KkRrLFq7mhc+nwHAG7dMJWI/H5RCqOXEqCimJqfwRlUlZ+XnsdRu541umaTpW954K8Ng4PVumSy12zkrP483qyq5OyWF4yKjOjl520lhJ8KXz0P5F4/gKlqHzhJN8sRH0MemqZ1KiJBw7zuvsbW4kKzUNB6bfK3acYRo0SVx8XxeW8PntbVscbt5oryMYo+HSbFxLV4/MSaWYo+HJ8rL2OJ283ltLV/U1nBpfHwnJ287KexEWFM8Tso+fRB32Vb0kfGkTHoUXVSC2rGECHp2p5Ornp0GwA0TzuXQg2SsqugckZGRREVFNR1Go7HF6wzAQWYzfzU03+50vr2BIXtpZR5isTDf3vz6PxsaGGg203IbX/CRwk6EPb/TRunM+/BUFaKPSSFl4qNorTLbT4j9+XnJQt6f/S1arZa3b78Xw166t4RoT+vXr6eurq7pmDp1aovXxer06DUaKr2+ZucrvT4S9boWX5Oo17d4vUGjIVbX8muCjRR2QgB+ew2lM+7FW1eGISGTlHMfRmOKUDuWEEHvlleep7SqkoN79GLqBZeoHUeEgX79+hEdHd10TJs2bZ/XKzTf21gD7Gu74z2u17Q1qTqksBNiB199OaUz7sXXUI0xpRfJZz+AxmBSO5YQQa26vo4bXnwGgHv+7zIOyu6pciLR1dlsNurr65sOt9vd4nU1Pi9eRSFxt5bkeL2OSp+vxddUeL17Xq/T4VEUavbymmAjhZ0Qu/BWF1E68z58ThvmbgeRdOY9oJPuJSH25dM5v/DNX3MxGgy8ffs9aLXy0SLU5wHWOJ2MsTbvfRljjWCFw9Hia1Y4HHtcf1hEBKudTrwdFbSdyXefELvxlOdR9ukD+N0OLD2GkTT+DtDIt4oQ+3Lt9Cepa7CRM/AQrjvjHLXjCAHA+9VVnB0by4ToGHoajdyZlEyawcDMmmoAbk5MYlrqv6shzKytIc1g4I6kZHoajUyIjuGsmFjeq6pS648QMPm0EqIF7qL1lH/xKIrXjbXfGBJOnkLjyAwhREsKy8u4442XAHj8imvpnpKqciIhYHZ9PdPKSrkmMZEvsrIZYbVy1fZtFHkb298S9XrSDIam6ws9Hq7evo1RVitfZGVzdUICj5eW8rOtXq0/QsA0wD6GEIa+9PR0CgsLycjIoKioSO04IsRYeo8i6cx70Gh11C39lupfXlc7khBBS6PR8Pvzr3PEkGHMXpTLSXfcqHYk0YXI53nrBHWLXfToc0i96Dkyp3xCt+s/JOnMe9DHZ6gdS4QRx6ZFVHz7HIriJ3r4qcQefpHakYQIWoqicMUzj+F0uzhxVA7/d9xJakcSIuwEdWFnzjyY+mXfUfLhbZTOvA+0OlLOfURmKopOZV/7B1U/vgpATM65RB96lsqJhAheG7cX8ND7bwMw/fpbSNrLCv9CiI4R1IVd2acP0PDPr3gqCvCUb6Xy++noY5IxpvRWO5oIM7aVs6n+/V0A4o68lMihJ6ucSIjg9czM/7J843oSYmKZfv0tascRIqwEdWG3O+2OBWP9TtterzEajc22GomMjOyseKKLq1v0BbXzGzc+Tzj+WiIOOlLdQEIEKa/Px+VPP4rP5+P8Y0/k5NGHqR1JiLARUoVd3NGX49y2Gk9F/l6vmTp1arOtRtavX9+JCUVXVzPvQ+qWfANAwik3Y+kzWuVEQgSnZRvW8ewnHwHw+i1TibLKTi5CdIaQKezij7saY3I2FbOe2ud106ZNa7bVSL9+/TopoQgX1b++he3vn9FodSSddifmrMFqRxIiKD34/ptsKtxGZnIK0664Tu04QoSFkCjs4o69CkvvQyn939346iv3ea3b7W621YjNtvduWyHaRqHyh5doWP8XGr2BpAn3Ycror3YoIYKOw+XiymceB+C6M8/hsEHyS1BXptVqOWLIMCYdfTxHDBkmO5CoJOi/6nHHXo217xhKZ9yDt7ZU7ThCNFL8VMx6GsfWZWiNZpLPfhBDcg+1UwkRdH5fvoS3v/sKgLdvvxeT0ahuINEhzhx3FHkzvmHO9Df43/2PMWf6G+TN+IYzxx2ldrSwE9SFXfxx1xA58EgqZj2N321HGxGLNiIWjV5+MIgg4PNS/uVjOLetRmuOJOXcR2SdRSFacPtrL1JcWUH/7tnce+FlascR7ezMcUfx2cNPkpGU1Ox8RmISnz38pBR3nSyod57IuvPbFs9XfPc8Df/82qp7yErVoqNpjFZSznscU2pvvHXllHx0B766crVjCRFUzhx3FF888hQer5fhV17I31s2qR1JtAOtVkvejG/ISEpC28Ke2n6/n+3lZfQ473T8fv8BvZd8nrdOUBd27UH+IYjOoLVEk3rBkxgSMvFUFVHy8R34G2rUjiVEUPnsoSc564ijWbxuDTnXXYbP51M7kjhARwwZxpzpb+z3uiOnXMUfK5Yd0HvJ53nrBHVXrBChwu+oo3TGvXhqSjDEp5My8VG0ZllDUYhdXf/CU1TX1zGy/0HcdNYkteOIdpAWn9iu14kDJ4WdEO3EZ6ukbMY9eOsrMSZlk3zOQ2iMFrVjCRE0Sqoque21FwB45LKr6ZkuY1JDXXFVRbteJw6cFHZCtCNvbSlln9yHz16LKb0fyRPuk8k+Quzi3e+/4deli7Gazbxx691qxxEHaMGaf2hwOvb6vN/vp6C0hHmrVnReqDAnhZ0Q7cxTUUDZpw/gd9kxZx1C4ul3gVandiwhgsaVzz6Gw+Xk2OGjuOTE8WrHEW0UYbHw9aPPEGG2oCgKfqX5kH2/3w8amPLycwc8cUK0nhR2QnQAd8kmyj57CL/HhbX3KBJPuQVamDEmRDjaUlTI/e82Drh/7roppMQnqJxIBCoxJpbfnnuNE0blYHPYeeSDdygsL2t2zfbyMs6+/y6+nPe7SinDk8yKFaIDmXsOJ3nCvWh0BupX/kjV7JfUjiREUNDpdCx49V1G9DuIT+f8wrkPTlU7kmil7imp/PT0y/TrnkVFbQ0n33kTi9etQavVMu6QIaTFJ1JcVcG8VSvataVOPs9bR5oQhOhAzi1LqZj1DIrfR9TgE4g9ShZnFQLA5/Nx+dOP4fV5OefIYzl97BFqRxKtcHCPXsx/+R36dc8iv6SYsTdczuJ1a4DGrtc/Vixjxm8/8ceKZdL9qhIp7IToYPb1f1G5o6UuZtQEYsbIMg9CAKzctIGn/vdfAF6dcicxkbJEUDA7bNBg5r74JhlJyfyzdTNjrp/M+oJ8tWOJ3UhhJ0QnaPj7F6p+eROA2HH/R9Tw01ROJERwePiDt1lfkE96YhJPXXWj2nHEXpyaM46fn3mZuKho/vx7BeNuuIKiCtlhJxhJYSdEJ6lf+g018z4EIP7YK4kYdJzKiYRQn8vt5opnHgPgyvFncsSQYSonEru75MTxfPnoU1hMZmbNn8fxt11Pja1e7VhiL6SwE6IT1c6fQe2iLwBIOPF6rP0OUzmREOqbt2o5r339GQBv3XYPZqNJ5URipzvPv5j37rofvU7Pu99/w5n33Y7D5VI7ltgHKeyE6GQ1v79L/YrZaLQ6EsffhrnncLUjCaG6O994me3lpfTp1p0HL7lC7ThhT6PR8Nx1N/PEldcD8MTH/2HyU4/I/r4hQAo7IVRQ9dOrNKz5A43OQNIZd2PqNlDtSEKoqt7ewDXPPQnArRMvYGiffionCl8GvZ7/3v0wN59zPgA3v/wcU998WeVUorWksBNCDYqfiu+ew75pEVqDieSzH8CY2lvtVEKo6tvcecz47Sf0Oj1v334vOp3s2NLZIiwWvnn8OS447kQ8Xi//99j9TP/sf2rHEgGQwk4Itfh9VHz9BM78lWhNVpLPfRhDYnfQaDFlDsI64HBMmYNkxwoRVm566Vmq6moZ1rc/t557gdpxwkpCTAy/PvsqJ47KocHhYPzdt/DRzz+oHUsESHaeEEJlGqOFlImPYkrvh89Rj+LzoI+Mb3reW1dO1a9v4tiQq2JKITrPRSecwn+mPojD5eSQy85nU+E2tSN1ed1TUvnx6Zfo3z2bitoaTrlrCovWrlY7VjPyed460hQghMoUt4OyTx/AU1uGzhKFLiKu2fO6qASSzrgbS98clRIK0bk++PE7flyUi8Vk5q3b70Gj0agdqUs7KLsnf738Nv27Z1NQWsLYGy4PuqJOtJ4UdkIEAb/LjkarQ1GUPT7ENBotoBB/zJXSLSvCxlXPTaPB4eDIIcO5/JQz1I7TZY05+BD+fOktuiWlsHrrFtlNoguQTwkhgoCp20D0UQl7bZnQaLToo5Nk9qwIG/klxdzzzmsAPH31jaQnJqmcqOs5JWcsvzz7CnFR0fz190rG3XgFheVlascSB0gKOyGCgC4ybv8XBXCdEF3BS1/MZMGav4mJjOSVKXeoHadLufjEU/nq0aexmMx8mzuP4267jur6OrVjiXYghZ0QQcBnq27X64ToCvx+P5c//Rhuj4czxh7J2Ucco3akLuH2SRfy/l0PoNfpeX/2t5x5r+wm0ZVIYSdEEHBtX423rhxF8e/1GkVRMKX1kXF2Iqys3rqZaR+9D8DLN91OXFS0uoFCmEaj4ZlrpvDU1TcC8OTH/+HSJx7CK7tJdCnyCSFEMFD8VP36JqDZo7hTFH/TpIq4oy4j9f+ewZCYpU5OIVTw+EfvsSZvCynxCTx77RS144QkvU7Hf6Y+yK0TG9cGvPXV6dwlu0l0SVLYCREkHBtyKf/qcXz1lc3O++orKf/qcSq+fwG/04YpvS9pl0wn5rDzQatXKa0Qncft8TD5qUfx+/1cetJ4jhk+Su1IIcVqNvPN489x4fEn4/F6ufCx+3nuk4/UjiU6iCxQLESw0WgxdRuILjIOn60a1/bVsKMVTxcZT/xx12DdsaaduzyPyu9fwF2yUc3EQnSKF2+8jRsmTGRLUSGDLpuE3elUO1LQi4+O4bsnnmf0QYNocDg4+4G7mL1ovtqx2kQ+z1tHWuyECDaKH9e2v7GvnYtr299NRR2Az1ZF+ZePUf71E/gaajAmZZN64TPEHnUZGr1JxdBCdLy733qVgtISeqZn8PClV6sdJ+hlJqfw50tvMfqgQVTW1nDMrdeGbFEnWk8KOyFCkH3dnxS9cy22f35Do9URM2oCaZe+hCnzYLWjCdFhbA47Vz07DYApZ09iRL+DVE4UvBp3k3iHAVk92FZWytgbrmDhmn/UjiU6gRR2QoQov6OOyu+eo+zTB/HWlWOITyf1/CeIP/5aNEaL2vGE6BCzF83nw59/QKfT8c4d92LQyzjT3eUMPIR5L75JZnIKa/K2MOa6yawryFM7lugkUtgJEeIcW5ZQ9M611C//AYCooSeTPvlVzD1HqJxMiI4x5eVnKa+p5pBefbjjvIvUjhNUTh59GL88+wrx0THkrl7FuBuvZHt5qdqxRCeSwk6ILkBxO6j66RVKPp6Kp7oIfXQSKec8SMIpt6A1R6kdT4h2VVlby00vPQvAfRdOpn/3bHUDBYkLjz+Zrx97BqvZzHe5f3LsrddRVVerdizRyaSwE6ILcW37m+J3b6B20Rcofh+RBx9N+uWvYe0/Vu1oQrSr//36I9/l/onJaOSt2+/Z6z7L4eK2if/HB3c/hF6n5z+zv+WMe2+TWcNhSgo7IboYxeui5vd3Kfnwdtzl+egiYkk6/S6SzrwHXYTsNSu6jmuef4J6ewNjBw3hmtPPVjuOKjQaDU9fcxNPX3MTAE/P+C+XPvmw7CYRxqSwE6KLchdvoPj9m6j582MUnxdr3xzSLn+NiEHHqh1NiHaxrayUu958BYAnrryOzOQUlRN1Lr1Ox/t3PcBtE/8PgNtee4E7Xn8RRenSy9OK/ZDCToiuzO+l9q+PKX7/JlzFG9CZI0k8eQrJ5z6MLjpZ7XRCHLDXvv6Mv/5eSZQ1gtduvkvtOJ3Gajbz1aPPcNEJp+D1ebno8Qd4duaHascSQUAKOyHCgKcin5L/3kb17+/g97iw9BhG+uRXiBp2Ko0b0AgRmhRF4fKnH8XldnNKzlgmHX282pE6XHx0DL88+wqn5IzF7nRy+j238d+fvlc7VkiL1mp5IjWNhb37sLB3H55ITSNKu/cSSQ/ckpjEV9nZLOnTlzm9ejEtNY0knfrL70hhJ0S4UPzULfqS4vduwLntH7RGC/HHXU3KBU+gj89QO50QbbauII9H//su0LjtWEJMjMqJOk63pBTmvfgmOQMPoaqulmNuuZbvF/yldqyQ93RaOv3NJq7cvo0rt2+jv9nEE2npe73erNVykNnM65WVnJ2Xx42FhWQbjbzSLTh+lipd+UhPT1cURVHS09NVzyKHHMFzaJTIIScpmVM+UbLu/FbpfusXSvShZytotEGQTQ45Aj8Mer2y6t3/KcqcxcoHdz+kep6OOPp3z1YKPvlWUeYsVgo++VYZkNVD9UydeXTU53lPo1FZ06+/cojZ3HTuELNZWdOvv5JtMLb6PgfveE2aXq/q10la7IQISwq2FT9Q9M51OLYsQaM3EnfkJaRe9ByG5B5qhxMiYB6vl8lPPYLP5+PC40/mhFE5akdqV6MPGsSfL71FZnIKa/O3Mua6yazN36p2LFVERkYSFRXVdBiNxgO63xCzhTqfj1W7LA+zyumkzudjqKX1u/hEabX4FYU6v3//F3cgKeyECGO++nLKPn2Qim+fw+eox5Tam7SLnid23P9BEIwVESIQi9et4YXPZwDwxi1TibRYVU7UPk46dAy/PvcqCTGxLFjzN2NvuCKsd5NYv349dXV1TcfUqVMP6H6Jej1VLSwPU+XzkajXteoeRo2Gm5OS+K6+jgaVCzv5yS2EoGH1bzi2LiP+uKuJ6D+WmDGTsPQdQ+UPL+IuWqd2PCFa7b53X+eMsUfSMz2Dxy6/pmmHilD1f8edxHt33Y9ep+f7BX9xzoN3hf3Cw/369aO4uLjpscvlavG66xISuS4xcZ/3OicvDwAFZY/nNNDC2T3pgWfT0tGi4eFS9QtuKeyEEAD47TVUfP0EDWtzSDjuWoyJ3Un9v6eoXzqLmrkfoHha/uEpRDCxO51c+ezj/PLsK1x/5rn879efWLDmb7Vjtckt517As9dOAeC/P33PZbLwMAA2m436+vr9XvdRTTXf19ft85pCj4d+JhMJLfRQxOl0VHr3/fXWA8+lZ5BhMHDptgLVW+tAumKFELtxbMil6J1rsP39MxqNlugRp5N22SuYswarHU2IVvl16SLe+2EWWq2Wd+64F6PBoHakgD151Q1NRd2zMz/i4mkPSlEXoBqfj61u9z4Pt6KwwukgWqdjkNnc9NpDzGaidTqWOxx7vf/Ooi7LaGTy9m3UBkFRB1LYCSFa4HfaqPz+BUo/uR9vbRmG2FRSJj1G/Ik3oDFFqB1PiP269dXplFRVclB2T+6+4FK147SaXqfjvbse4I7zLgLg9tde4LbXpstuEh1oi9vNPJuNh1JSOcRs5hCzmYdSU/ndZiPP42667tvsHhwTGQmADpiensFAs5k7iovQAYk6HYk6HWr/GiGFnRBir5xbl1H07nXULf0WgKjBJ5A++VUsvUepnEyIfauur+OGF54GYOoFlzCwRy+VE+2fxWTiy0ef5pITT8Xr83LJEw/xjOwm0SnuKC5io9vF290yebtbJhtcLu4qLmp2TU+TiSht42SKFL2Bo6OiSDMY+DK7B3N792k6hqg8aae1YwNDVnp6OoWFhWRkZFBUVLT/FwghWmTqNpCEk27EsGMx44Y1f1D1yxv4HfsewyKEmr569BlOH3sEC9b8zWHXX44/SLrLdhcXFc23055nzMGHYHc6OfehqXyX+6fasYKKfJ63jrTYCSFaxbV9NcXv3UDtgk9R/D4iDjqC9MtfwzrgcLWjCbFX105/klqbjdEHDeL6M89VO06LMpKSmffiW4w5uHE3iWNvvU6KOtFmUtgJIVpN8bqp+eM/lHxwC+6yreisMSSddgdJZ92PLjJB7XhC7KGoopw73ngRgMcvv5as1DSVEzXXv3s2819+h4E9erK9vJRxN15J7upVascSIUwKOyFEwNylmyn+z83UzP0viteDtfco0i9/lcjBJ6gdTYg9vPXtV/yxYhkRFgtv3HJgi9m2p1EDBvLnS2/RPSWVdQV5jLluMmvytqgdS4Q4KeyEEG3j91KbO5Pi92/EVbQOrSmChBNvIGXSY+hjU9VOJ0QTRVG44pnHcLpdnDAqhwuPP1ntSJw4agy/PfcaCTGxLFzzD2NvuJxtZeovbitCnxR2QogD4qncRsmHd1D161v4PU7MWYNJu+xlokacDhr5ESOCw8btBTz4/lsAPH/dzSTFxqmW5YLjTuKbx58lwmJh9qJcjr7lGipra1XLI7oW+akrhDhwip/6JV9T/M51OPJWojWYiT/mClIveApDQqba6YQA4NmZH7J843oSYmJ58cbbVMlw8znn8+E9D2PQ6/nw5x847e5bwn6LMNG+pLATQrQbb20pZTPvoXL2S/hdDZgy+pN2yYvEjJkE2tZtpi1ER/H6fEx+6hG8Pi+Tjj6eU3PGder7P3Hl9Tx33c0APPfJR1z0+AN4vN5OzSC6PinshBDtzrbyR4reuRb7pkVo9AZix/0faRdPx5gS/IvEiq5t+cb1PDvzIwBeu+VOoiM6ficVnU7Hu3fez53nXwzAnW+8xK2vym4SomNIYSeE6BC++krKP3+Y8m+ewmevxZjcg9SLniP2iIvR6I1qxxNh7MH332JT4Ta6JaXwxJXXd+h7WUwmvnzkaS49aTxen5dLn3iYp/73QYe+pwhvUtgJITqUfe1cit6+hoY1f6DR6ogZfQ5pl7yIqdtAtaOJMOV0u7ji6ccAuOb0sxk7aEiHvE9cVDQ/PfMy48eMw+Fycua9d/D+7Fkd8l5C7CSFnRCiw/kddVTMepqyzx/GW1+JIaEbqRc8SdyxV6MxWtSOJ8LQnBVLeevbrwB4+/Z7MRnbtxU5IymZuS++ydhBQ6iur+O4267n29x57foeQrRECjshRKdxbFpE0TvXUr/yRwCih59K+mWvYO4xTOVkIhzd/voLFFWU0697FvddOLnd7tuvexZ/vfQ2B/foRWF5GeNuvJK//l7ZbvcXYl+ksBNCdCrF1UDV7JconXEPnpoS9DHJpJz7MAknT0FrjlQ7nggjtTYb101/CoA7zruIQ3r1OeB7Nu4m8TZZqWmNu0lcP5nVWzcf8H2FaC0p7IQQqnDmr6T43euoW/I1iuInctCxpE9+DUvfHLWjiTDy1Z9z+OyPXzHo9bxzx73odG1fluf4kaP57bnXSIyJZdHa1Yy94XIKSkvaMa0Q+yeFnRBCNYrHRfWvb1H60R14Krehi4wj+cx7SDxjKtqIWLXjiTBx/QtPU11fx4h+BzHl7PPadI/zjz2Rb6c9T4TFwo+ym4RQkRR2QgjVuQrX8f/t3XtYVHX+B/D3zHBREVS8gGDegcUwMNDaSm1VNq1fYP00ts3CSq3MTLMEtFLXRZ8yMlL6uZli9nNN1y2l3FVDRRRcUyEUb91A5SIaIgwwwMB89g90ZABl8HYOw/v1PJ9Hzplzec/3OTgfzpwzkxc/HcWpX0JMNXDyeRAeL/4fnO4eoXQ0agUKLhZi1icfAQD+8vxL6Ovh2az1Xx/3NNa9vRD2dnb4e+I2PD7nDZQZDLchKVHT2NgRkTrUGHFp7/8j//MZqDz3M3RtndHlf95At/HzoXPuqnQ6snHx//4GiYe/R7s2bfDprLlWr7do8qv4aNobAICPNq3HhOh3+W0SpCg2dkSkKsbzWTj3xSwUJa2BVFehbd8geLwYh/aDHgWgUToe2bCXYhahvKICIwMH44VHQ667rE6nw2dvvY2oZyYCACI/XY6Zyz/kt0mQ4tjYEZH6mGpQcmAT8uKnoyLnGLSO7dD5j1Ph9vQi2HXyUDod2ahf83LxbvzfAAAxU2fA3bVzo8u1dXTEV395Hy8+Foqay98/+97fP7+TUYmuSQPApv+88PDwQG5uLjw9PZGXl6d0HCJqNg2c730MHYeHQ+vQFiZjJYr3rUPJwc2AmACNFo497oaufSfUlBahMudY7XxqHMfrunQ6Hf7zyWoE+QzApj278Om+g/Bw90TeuVzs3LUZLm3b4pvFH+KhgQEwVFYgbMFcfJOarHRsdbjNxxZfz63TIhq79oMeRYchT0LX3hVVv51B0c6VtQeMFXggENkGnUs3dB49DW0vf5hxZf6PKDuWBJchT8DO5eo1eNUlF3Bx56cw/Lhfqaiq1db793AdOYXj1YR7+nnh8KdfwK7eR58UGIyoLNejZ2dXXCrV4/GoN7Dv6A/KhFSZO3Fs8fXcOqp/K7bd74bCdeRkFO/fiLw101GZc4wXUxO1QjUl53F+47v47V8foaaiFI7dvdFp5GTonLtYLKdz7oyuY+fw8/Dqaev9e3QdO4fjZYUBQSOg0zZ8eezWxg49O7uiqKwMQ1+bzKbuMh5b6mKndICmuAwei9Ij36H0yA4AQNHOlWjb5144D3oUl5J5TQNRa1N2NBEVWWnwmLISWnvHBo9rNFqICDo/Mg2/VRtr3woSMf8rV34GADFdnRYBIJeXubJO7TyLaYvtCICrj1lsS0y1b4dI3cdN5n2gbo7bTaOF68gpl8fH8gaU2vEywXXkFOT+dKDVvy2r1dkhJvw5CBreqqPRaCAiqNQ44PiZ00rEUx8eW6qj7sZOawcH9/4o/s8mi9mGrHQ4ev6u0VUcHBzg6Hj1P/v27fkVRUS2xq6TZ6NN3RUajQa6dh3gNn7+nQt1E6RB82mqvUjmciNo2TDWbxKvLisWDePV5lNj52DxFll9Go0Wdi5d0f3FOEhV3c9fs+IuZM2tu1O5fmNwnSWt2dgNZXjAsws82jlcZ7MauLezxxNzvkBqbmHTeZqMcf0FrB+TG9j5dbfd1LqX/7FzgF29M3WWu6g9thx73I3Ks0evv026JVTd2OnauUCj1cFUXmQxv6asCDqnxr80PCoqCvPnz78D6YhIKbr2naxazliUD1NlWe2Lo0Zb+0Km0QDQ1JunrX2dqjMNjQYaaK6uo9FeXafOspq6270y3Uy12wCA2mu6lPpQF4fOdym0Z/Xw6NLGuuU6d4VDlfNtTmM7rP2dpZun6sbuigYfC6S59j0fixcvxocffmie7t69O06dOnX7whHRHVdTWtT0QgAK//2xQmcJNA2bvXpNYf1ms7aJ1DZsEOs2n3Ua06uN5dXG0zwNy2bTwa0fXEdMajL1xaQ1MF7IbuSRZtxjd7s+x83q7d5c1p/vHw4EPNvkqj/vT0DBf/Zcd1tXEzWRqcnI11ngZsf7uus3tW2pPbaCX2lyN9b+ztLNU3VjV1NeAjHVQOdk2enr2nVETdmlRtepqqpCVVWVedrZmX9REdmaypxjqC65AJ1z50bPkImYUKMvtPru+Vvvylui5qnGlrhjKs8eg0tQaJPjpf/+q1Z/HdT2nGPICw+De1t7aBt5q9IkgvxyI7Zv+gSmGn7DRGXej3C5b5yKfxdbH3XfFWuqRtW5n9G2d4DF7Da9A1CZe1KZTESkPDHh4s5PAWguX2tW5yExAdDUPt7KmxQzjpfVTDXVmPX52tqf653NujL95tq1bOqu4LGlOupu7ACUHNyM9v5/hNPAYNh17oFOIybBzqUr9D/8S+loRKQgw4/7cWHzItToCy3m1+gLcWHzIn4uWz0cL+t9ueFveGbFKpwzGC3m55cb8cyKVfhyw98USqZOPLbUpeV8QPF9/wudkyuqfjvNDygmoqv4TQrNw/GymlZnh5Ejxlp88wTP1F0Hv3lCFVpEY3czeCAQERG1fHw9t47q34olIiIiIuuwsSMiIiKyEWzsiIiIiGwEGzsiIiIiG8HGjoiIiMhGsLEjIiIishFs7IiIiIhsBBs7IiIiIhvBxo6IiIjIRrCxIyIiIrIRbOyIiIiIbISd0gHuFDc3N6UjEBER0Q3i67h1bL6xu3IgpKWlKZyEiIiIbpabmxvy8vKUjqFaGgCidIjbbdCgQSgoKLjl223fvj1OnToFHx8flJaW3vLt2xqOl/U4Vs3D8Woejpf1OFbNc7vHy83NDenp6bd8u7akVTR2t4uzszNKSkrg4uICvV6vdBzV43hZj2PVPByv5uF4WY9j1TwcL+Xx5gkiIiIiG8HGjoiIiMhGsLG7CZWVlZg/fz4qKyuVjtIicLysx7FqHo5X83C8rMexah6Ol/J4jR0RERGRjeAZOyIiIiIbwcaOiIiIyEawsSMiIiKyEWzsiIiIiGwEG7sbMHToUCQkJCA3NxcigtDQUKUjqVZkZCS+//57lJSUoKCgAF9//TW8vb2VjqVaL7/8MjIyMlBcXIzi4mKkpqZi9OjRSsdqESIjIyEiWLp0qdJRVGnevHkQEYvKz89XOpaqeXh44IsvvsBvv/2GsrIypKen495771U6liplZWU1OL5EBMuXL1c6WqvDxu4GODk5ISMjA9OmTVM6iuoNHz4ccXFxuP/++xEcHAw7Ozvs2LED7dq1UzqaKuXk5CAyMhJBQUEICgrCrl27sGXLFgwYMEDpaKoWFBSEKVOmICMjQ+koqpaZmQl3d3dzDRw4UOlIqtWxY0ekpKTAaDRizJgxGDBgAGbNmoVLly4pHU2VBg8ebHFsjRo1CgDwj3/8Q+FkrZOwbrxEREJDQxXP0VKqS5cuIiIydOhQxbO0lCosLJQXXnhB8RxqLScnJzl16pSMHDlSdu/eLUuXLlU8kxpr3rx5kp6erniOllKLFy+W5ORkxXO01Fq6dKn89NNPiudojcUzdnRHdejQAQBw8eJFhZOon1arRVhYGJycnLB//36l46hWXFwctm7dip07dyodRfW8vLyQm5uLX3/9FevXr0efPn2UjqRaISEhOHToEDZu3IiCggKkpaVh0qRJSsdqEezt7TFhwgSsXr1a6SitluLdZUsunrFrXm3ZsoV/BTdRfn5+otfrxWg0SlFRkYwZM0bxTGqtsLAwOXLkiDg6OgoAnrG7To0ePVqefPJJ8fPzM5/dzM/PF1dXV8WzqbEMBoMYDAaJjo6WgIAAmTJlipSXl8uzzz6reDa11/jx48VoNEr37t0Vz9JKS/EALbrY2Flfy5cvl6ysLPH09FQ8i5rL3t5e+vXrJ4GBgbJo0SI5f/68+Pr6Kp5LbdWjRw85d+6c3HPPPeZ5bOysr3bt2kl+fr7MnDlT8SxqrMrKSklJSbGYFxsbK6mpqYpnU3tt27ZNEhISFM/RikvxAC262NhZVx9//LGcOXNGevfurXiWllbfffedrFixQvEcaqvQ0FARETEajeYSEampqRGj0SharVbxjGqvHTt2yCeffKJ4DjVWdna2rFy50mLeyy+/LDk5OYpnU3P17NlTqqurJSQkRPEsrbXsQHSbLVu2DE888QQefvhhZGdnKx2nxdFoNHB0dFQ6hurs3LkTfn5+FvPi4+Nx8uRJvPfeezCZTAolaxkcHBzg6+uLvXv3Kh1FlVJSUuDj42Mxz9vbG6dPn1YoUcvw/PPP4/z589i6davSUVo1xbvLllZOTk7i7+8v/v7+IiIyY8YM8ff3l7vuukvxbGqruLg4KSoqkmHDhombm5u52rRpo3g2NVZ0dLQ89NBD0qtXL/Hz85O//vWvUl1dLaNGjVI8W0sovhV77VqyZIkMGzZMevfuLUOGDJGEhAQpLi6Wnj17Kp5NjRUUFCRVVVUSFRUl/fr1k6efflpKS0vlz3/+s+LZ1FoajUays7Nl8eLFimdp5aV4gBZXw4cPl8bEx8crnk1tdS3h4eGKZ1NjffbZZ5KVlSUVFRVSUFAg3333HZu6ZhQbu2vX+vXrJTc3VyorKyUnJ0c2bdrEazebqMcee0yOHDkiBoNBjh8/LpMmTVI8k5orODhYRES8vLwUz9KaS3P5ByIiIiJq4fg5dkREREQ2go0dERERkY1gY0dERERkI9jYEREREdkINnZERERENoKNHREREZGNYGNHREREZCPY2BERERHZCDZ2RHTTevXqBRGBv7+/0lHMfHx8sH//fhgMBqSnpysdh4jojmBjR2QD4uPjISKIiIiwmB8aGgqR1vnlMgsWLEBZWRl8fHwwcuTIRpfp2rUrVqxYgdOnT6OiogL5+fnYtm0b7r///juclojo1mBjR2QjDAYDIiIi0LFjR6Wj3DL29vY3vG6/fv2wb98+nDlzBhcvXmx0mX/+85/w9/dHeHg4vL29ERISgqSkJLi6ut7wfq1xM8+LiKgpin9hLYvFurmKj4+XhIQEOX78uLz33nvm+aGhoSK1p+wEgMybN0/S09Mt1n399dclKyvLYltff/21REVFyblz56SoqEjeffdd0el08v7770thYaGcPXtWnn/+efM6vXr1EhGRsLAwSUlJEYPBIJmZmTJ8+HCLffn6+srWrVtFr9fLuXPnZO3atdK5c2fz47t375Zly5ZJTEyMXLhwQZKSkhp9vhqNRt555x05e/asVFRUSHp6ujzyyCPmx+ubN29eg2106NBBRESGDRt23bG96667ZPPmzaLX66W4uFg2bNgg3bp1azBedddZunSp7N69u8nnNWDAAPn222+luLhYSkpKJDk5Wfr27Wteb+LEiXL8+HExGAxy4sQJeeWVV8yP2dvby7JlyyQvL08MBoNkZWVJZGSk4scii8VStnjGjshG1NTUYM6cOXjttdfg6el5U9saMWIEPDw8MGzYMLzxxhtYsGABvv32WxQVFeG+++7DihUrsGLFCvTo0cNivSVLliAmJgaDBg1CamoqEhISzGe/3N3dsWfPHvzwww8ICgrC6NGj4ebmho0bN1psIzw8HNXV1XjwwQfx0ksvNZrv9ddfx6xZs/Dmm2/innvuwfbt25GQkID+/fub95WZmYkPPvgA7u7u+OCDDxpso7S0FHq9HmPHjoWDg8M1x2Lz5s1wdXXF8OHDERwcjH79+mHDhg3NGs/GnpeHhweSk5NRUVGBESNGIDAwEKtXr4adnR0AYNKkSYiOjsbcuXPh6+uLOXPmYOHChXjuuecAANOnT0dISAieeuop+Pj4YMKECcjOzm52LiKyPYp3lywW6+aq7lmj1NRU+eyzzwS48TN2WVlZotFozPNOnDghe/bsMU9rtVrR6/USFhYmwNUzdrNnzzYvo9Pp5MyZM/LWW28JAFmwYIFs27bNYt+enp4iIuLl5SVA7ZmttLS0Jp9vTk6OREVFWcw7cOCALF++3Dydnp7e6Jm6uvXkk09KYWGhlJeXy759+yQ6OloGDhxofnzUqFFiNBqlR48e5nm+vr4iIhIUFNRg7K9UY2fs6j+v6Oho+eWXX8TOzq7RbKdPn5Y//elPFvPmzp0rKSkpAkBiY2MlMTFR8WOPxWKpq3jGjsjGREREIDw8HL6+vje8jWPHjlncdFFQUICjR4+ap00mEwoLC9GtWzeL9fbv32/+uaamBocOHTLnCAwMxB/+8Afo9XpznTx5EkDt9XBXHDp06LrZnJ2d4enpiZSUFIv5KSkpzX7OX331FTw8PBASEoLt27fj4YcfRlpaGsLDwwEAvr6+OHv2LHJycszrnDhxAkVFRc3eV/3nFRAQgL1796K6urrBsl26dEHPnj2xatUqi/F6++23zWO1Zs0aBAQE4NSpU4iNjUVwcHCz8hCRbbJTOgAR3Vp79+7F9u3bsWjRIqxZs8biMZPJBI1GYzGvsQv5jUajxbSINDpPq236b8MrDaJWq8U333zT4M5dAMjPzzf/XFZW1uQ26273Co1Gc0N3AFdWViIxMRGJiYlYuHAhVq5ciQULFuDzzz+/5jbrzrd2TOs/L4PBcM1MV8Z18uTJOHDggMVjNTU1AID09HT06dMHY8aMwahRo7Bx40YkJiZi/PjxVjxrIrJVPGNHZIMiIyPx+OOP44EHHrCYf+HCBbi7u1vMCwgIuGX7rfsxITqdDoGBgeazcmlpabj77ruRnZ2NX375xaLKy8ut3oder0dubi4eeughi/kPPPAATpw4cdPP4fjx43BycjL/3LNnT4trCX19fdGxY0fzvi5cuIDu3btbbMOaMT1y5AiGDh1qvqaurvPnzyMnJwd9+/ZtMFZ1r6PT6/XYuHEjpkyZgrCwMIwbNw6dOnW6gWdNRLaCjR2RDcrMzMS6devw2muvWcxPSkpC165dMXv2bPTt2xdTp07FmDFjbtl+X331VYwdOxY+Pj6Ii4tDp06dsHr1agBAXFwcXF1dsX79egwePBh9+vRBcHAwVq1aZdWZv7qWLFmCiIgIPPXUU/D29sbixYsREBCA2NhYq7fh6uqKnTt34plnnsHAgQPRu3dvjBs3DrNnz8aWLVsAAImJiThy5AjWrVuHQYMGYfDgwVi7di2SkpJw+PBhAMCuXbsQFBSEZ599Fv3798f8+fPh5+fX5P6XL18OFxcXfPnllwgMDET//v0xYcIEeHt7AwDmz5+PqKgoTJ8+HV5eXvDz88PEiRMxc+ZMAMCMGTMQFhYGHx8feHl5Yfz48cjPz8elS5eaNZZEZFvY2BHZqHfeeafBW4QnT57E1KlT8eqrryIjIwNDhgxp9I7RGxUZGYmIiAhkZGRg6NChCA0NRWFhIYDat1sffPBB6HQ6bN++HZmZmYiNjUVxcTFMJlOz9vPxxx8jJiYGMTExOHr0KEaPHo2QkBD8/PPPVm+jtLQUBw4cwMyZM5GcnIzMzEzzW7HTpk0zLzd27FgUFRUhOTkZiYmJ+PXXXxEWFmZ+fMeOHVi4cCHef/99HDx4EM7Ozli7dm2T+7948SJGjBiB9u3bY8+ePTh8+DAmT55sfst71apVmDRpEiZOnIijR49iz549mDhxIrKyssz5IyIicOjQIRw8eBC9e/fGo48+2mo/kJqIamlQexcFEREREbVwPGNHREREZCPY2BERERHZCDZ2RERERDaCjR0RERGRjWBjR0RERGQj2NgRERER2Qg2dkREREQ2go0dERERkY1gY0dERERkI9jYEREREdkINnZERERENuK/yekeMAoMXtIAAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 640x480 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "numsources = Forclusts['num_source'].values\n",
        "avgsilsall    = Forclusts['avg_sillhouette_score'].values\n",
        "minsillsall   = Forclusts['min_sillhouette_score'].values\n",
        "recs_all   = Forclusts['reconstruction_loss'].values\n",
        "plot_results_stability_plot(numsources,avgsilsall,minsillsall,recs_all)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.12473557,  0.69991903,  0.38559303,  0.38979156,  0.50006274],\n",
              "       [-0.12817872,  0.30006231, -0.10364522, -0.24311644,  0.09998255],\n",
              "       [ 0.24967361,  0.30001035,  0.215059  ,  0.43861195, -0.09998018]])"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# good_res.groupby('num_sources', group_keys=False).apply(lambda x: result_analysis(x, observations.size))['centers'].iloc[1]\n",
        "\n",
        "Forclusts['centers'].iloc[4]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "X2UoxPrhIO0D",
        "outputId": "62956299-0baa-47f7-811a-9fa76565a2b7"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.4 ,  0.3 ,  0.2 ,  0.5 ,  0.7 ],\n",
              "       [-0.25, -0.08, -0.15,  0.1 ,  0.3 ],\n",
              "       [ 0.45,  0.2 ,  0.25, -0.1 ,  0.3 ]])"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "ground_truth=np.vstack([As,Xs.T])\n",
        "ground_truth"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.7 ,  0.39,  0.12,  0.39,  0.5 ],\n",
              "       [ 0.3 , -0.24, -0.13, -0.1 ,  0.1 ],\n",
              "       [ 0.3 ,  0.44,  0.25,  0.22, -0.1 ]])"
            ]
          },
          "execution_count": 17,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Forclusts['centers'].iloc[4].round(2)[:,[1,3,0,2,4]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "33FA_RTWuC34"
      },
      "source": [
        "# jax graph analysis"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "-v4945FCwO-J",
        "outputId": "0757949d-afea-42c0-c3a1-165927e68e84"
      },
      "outputs": [],
      "source": [
        "# # # import tensorflow as tf\n",
        "# # # from jax.experimental import jax2tf\n",
        "\n",
        "# # # f_tf = jax2tf.convert(green_matrix)\n",
        "# # # ff_tf = tf.function(f_tf)\n",
        "\n",
        "# # # logdir = 'tensorboard-test2/2'\n",
        "# # # writer = tf.summary.create_file_writer(logdir)\n",
        "# # # tf.summary.trace_on(graph=True, profiler=True)\n",
        "# # # ff_tf(As, Xs, Ts, Xd, t, u_scaler, D)\n",
        "# # # with writer.as_default():\n",
        "# # #   tf.summary.trace_export(\n",
        "# # #       name=\"my_func_trace\",\n",
        "# # #       step=0,\n",
        "# # #       profiler_outdir=logdir)\n",
        "\n",
        "# resid = jax.jit(make_flat_resids_fn(model_fn, observations, num_sources))\n",
        "# flat_init, _ = flatten_inputs(*gen_init(num_sources))\n",
        "# z=jax.xla_computation(resid)(flat_init)\n",
        "# graphviz.Source(z.as_hlo_dot_graph())\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "authorship_tag": "ABX9TyPC9nWDImGL6+nVkb9Bo85F",
      "gpuType": "T4",
      "include_colab_link": true,
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.10.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}