ECAI_Exp2.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "6eXzKXL9J76b"
      },
      "outputs": [],
      "source": [
        "%matplotlib inline\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from tqdm.notebook import trange, tqdm\n",
        "\n",
        "ETA_Q = 5\n",
        "ETA_R = 5\n",
        "EPS = 1e-8\n",
        "STEP = 0.2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g4GAQ9f1J76e"
      },
      "source": [
        "# The power of contrastive examples\n",
        "\n",
        "We compare RL, contrastive examples and demonstrations."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RgcJDZxsJ76g"
      },
      "source": [
        "## Environment definition\n",
        "\n",
        "The environment looks like this:\n",
        "\n",
        "```\n",
        " ---- ---- ---- \n",
        "|  1 |  2 |  3 |\n",
        " ---- ---- ---- \n",
        "|  4 |    |  5 |\n",
        " ----      ----    \n",
        "|  6 |    |  7 | \n",
        " ----      ----  \n",
        "|  8 |    |  9 |\n",
        " ---- ---- ---- \n",
        "| 10 | 11 | 12 |\n",
        " ---- ---- ---- \n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "MhM-GZkGJ76g",
        "outputId": "00b5ef94-1f08-483d-9c71-5ebbd714b314"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "States: ['1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12']\n"
          ]
        }
      ],
      "source": [
        "S = list(map(lambda x : str(x), range(1, 13)))\n",
        "\n",
        "print('States:', S)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UlZgXAGoJ76h",
        "outputId": "02043b78-10ab-4b23-9549-40e7e4410c3f"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Actions: ('Up', 'Down', 'Left', 'Right', 'No-Op')\n"
          ]
        }
      ],
      "source": [
        "A = ('Up', 'Down', 'Left', 'Right', 'No-Op')\n",
        "\n",
        "print('Actions:', A)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7zPS3eBhJ76i",
        "outputId": "21e188fd-bbbb-41a0-9dc7-47d11ffb4cd0"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Transition probabilities ok?\n",
            "\t- Action Up: True\n",
            "\t- Action Down: True\n",
            "\t- Action Left: True\n",
            "\t- Action Right: True\n",
            "\t- Action No-Op: True\n"
          ]
        }
      ],
      "source": [
        "#               1   2   3   4   5   6   7   8   9   10  11  12               \n",
        "Pu = np.array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 1\n",
        "               [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 2\n",
        "               [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 3\n",
        "               [.8, 0., 0., .2, 0., 0., 0., 0., 0., 0., 0., 0.],  # 4\n",
        "               [0., 0., .8, 0., .2, 0., 0., 0., 0., 0., 0., 0.],  # 5\n",
        "               [0., 0., 0., .8, 0., .2, 0., 0., 0., 0., 0., 0.],  # 6\n",
        "               [0., 0., 0., 0., .8, 0., .2, 0., 0., 0., 0., 0.],  # 7\n",
        "               [0., 0., 0., 0., 0., .8, 0., .2, 0., 0., 0., 0.],  # 8\n",
        "               [0., 0., 0., 0., 0., 0., .8, 0., .2, 0., 0., 0.],  # 9\n",
        "               [0., 0., 0., 0., 0., 0., 0., .8, 0., .2, 0., 0.],  # 10\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],  # 11\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., .8, 0., 0., .2]]) # 12\n",
        "\n",
        "#               1   2   3   4   5   6   7   8   9   10  11  12               \n",
        "Pd = np.array([[.2, 0., 0., .8, 0., 0., 0., 0., 0., 0., 0., 0.],  # 1\n",
        "               [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 2\n",
        "               [0., 0., .2, 0., .8, 0., 0., 0., 0., 0., 0., 0.],  # 3\n",
        "               [0., 0., 0., .2, 0., .8, 0., 0., 0., 0., 0., 0.],  # 4\n",
        "               [0., 0., 0., 0., .2, 0., .8, 0., 0., 0., 0., 0.],  # 5\n",
        "               [0., 0., 0., 0., 0., .2, 0., .8, 0., 0., 0., 0.],  # 6\n",
        "               [0., 0., 0., 0., 0., 0., .2, 0., .8, 0., 0., 0.],  # 7\n",
        "               [0., 0., 0., 0., 0., 0., 0., .2, 0., .8, 0., 0.],  # 8\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., .2, 0., 0., .8],  # 9\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],  # 10\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],  # 11\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]]) # 12\n",
        "\n",
        "\n",
        "#               1   2   3   4   5   6   7   8   9   10  11  12               \n",
        "Pl = np.array([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 1\n",
        "               [.8, .2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 2\n",
        "               [0., .8, .2, 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 3\n",
        "               [0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],  # 4\n",
        "               [0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],  # 5\n",
        "               [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],  # 6\n",
        "               [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],  # 7\n",
        "               [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],  # 8\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],  # 9\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],  # 10\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., .8, .2, 0.],  # 11\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .8, .2]]) # 12\n",
        "\n",
        "\n",
        "#               1   2   3   4   5   6   7   8   9   10  11  12               \n",
        "Pr = np.array([[.2, .8, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 1\n",
        "               [0., .2, .8, 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 2\n",
        "               [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],  # 3\n",
        "               [0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],  # 4\n",
        "               [0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],  # 5\n",
        "               [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],  # 6\n",
        "               [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],  # 7\n",
        "               [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],  # 8\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],  # 9\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., .2, .8, 0.],  # 10\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .2, .8],  # 11\n",
        "               [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]]) # 12\n",
        "\n",
        "Pn = np.eye(len(S))\n",
        "\n",
        "P = [Pu, Pd, Pl, Pr, Pn]\n",
        "\n",
        "print('Transition probabilities ok?')\n",
        "\n",
        "for aidx in range(len(A)):\n",
        "    print('\\t- Action %s:' % A[aidx], np.all(np.isclose(P[aidx].sum(axis=1), 1.)))\n",
        "    \n",
        "gamma = 0.99"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "utyA_tD3J76k"
      },
      "outputs": [],
      "source": [
        "# Rewards are a convex combination of 7 features/rewards\n",
        "CELL_GOAL = np.zeros((len(S), len(A)))\n",
        "CELL_GOAL[2, :] = 1\n",
        "\n",
        "CELL_6 = np.zeros((len(S), len(A)))\n",
        "CELL_6[7, :] = 1\n",
        "\n",
        "CELL_7 = np.zeros((len(S), len(A)))\n",
        "CELL_7[8, :] = 1\n",
        "\n",
        "\n",
        "RFEAT = (CELL_GOAL, CELL_6, CELL_7)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zcthHE-sJ76k"
      },
      "source": [
        "## MDP solver\n",
        "\n",
        "The next function is used to compute the optimal $Q$-function for MDPs. For computational efficiency, we pre-compute the optimal $Q$-function for all MDPs."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "id": "KOVzIT71J76l"
      },
      "outputs": [],
      "source": [
        "def vi(M):\n",
        "    ''' \n",
        "    Compute optimal Q-function for MDP M.\n",
        "\n",
        "    :param M: 5-tuple corresponding to MDP description: (S, A, P, R, gamma)\n",
        "    :return: np.array with shape (|S|, |A|), corresponding to Q* for M.\n",
        "    '''\n",
        "\n",
        "    S = M[0]\n",
        "    A = M[1]\n",
        "    P = M[2]\n",
        "    R = M[3]\n",
        "    gamma = M[4]\n",
        "    \n",
        "    Q = np.zeros((len(S), len(A)))\n",
        "    quit = False\n",
        "    niter = 0\n",
        "    \n",
        "    while not quit:\n",
        "        Qnew = np.zeros((len(S), len(A)))\n",
        "        V = np.max(Q, axis=1, keepdims=True)\n",
        "        for a in range(len(A)):\n",
        "            Qnew[:, a, None] = R[:, a, None] + gamma * P[a].dot(V)\n",
        "            \n",
        "        if np.linalg.norm(Q - Qnew) < EPS:\n",
        "            quit = True\n",
        "            \n",
        "        Q = Qnew\n",
        "        niter += 1\n",
        "        \n",
        "    # print('Finished after %i iterations.' % niter)\n",
        "    return Q\n",
        "\n",
        "# --\n",
        "\n",
        "def greedy(Q):\n",
        "    ''' \n",
        "    Compute greedy policy with respect to Q.\n",
        "    \n",
        "    :param Q: np.array with shape (|S|, |A|).\n",
        "    :return: np.array with shape (|S|, |A|), corresponding to pi_g(Q).\n",
        "    '''\n",
        "    \n",
        "    pol = (Q == np.max(Q, axis=1, keepdims=True)).astype(int)\n",
        "    pol = pol / np.sum(pol, axis=1, keepdims=True)\n",
        "\n",
        "    return pol\n",
        "\n",
        "# --\n",
        "\n",
        "def boltzmann(Q):\n",
        "    ''' \n",
        "    Computes the Boltzmann policy with respect to Q.\n",
        "\n",
        "    :param Q: np.array with shape (|S|, |A|).\n",
        "    :return: np.array with shape (|S|, |A|) corresponding to pi_B(Q).\n",
        "    '''\n",
        "    \n",
        "    Q = ETA_Q * (Q - np.max(Q, axis=1, keepdims=True))\n",
        "\n",
        "    pol = np.round(np.exp(Q), 3)\n",
        "    pol = pol / np.sum(pol, axis=1, keepdims=True)\n",
        "    \n",
        "    return pol\n",
        "\n",
        "# --"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "id": "M462WAssJ76m"
      },
      "outputs": [],
      "source": [
        "w = np.random.randn(4)\n",
        "r = np.sum([RFEAT[m] * w[m] for m in range(len(RFEAT))], axis=0)\n",
        "M = (S, A, P, r, gamma)\n",
        "Q = vi(M)\n",
        "pol = boltzmann(Q)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zHt6k7pbJ76m"
      },
      "source": [
        "## RL updater\n",
        "\n",
        "We perform maximum likelihood estimation using SGD, where \n",
        "\n",
        "$$r_t\\sim\\mathrm{Normal}(w^\\top\\phi(s_t,a_t);\\eta)$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "id": "gGr9dYM7J76m"
      },
      "outputs": [],
      "source": [
        "# Weight update\n",
        "def rl_update(wold, s, a, r):\n",
        "    ''' \n",
        "    Updates the weights wold given the triplet (s, a, r) and the weight w. \n",
        "    The update is a Stochastic Gradient Descent (SGD) update on the negative log likelihood, assuming the\n",
        "    reward is disturbed by 0-mean Gaussian noise.\n",
        "    \n",
        "    :param wold: np.array with shape (nfeat,)\n",
        "    :param s: int (state index)\n",
        "    :param a: int (action index)\n",
        "    :param r: float (reward)\n",
        "    :return: np.array with shape (nfeat,)\n",
        "    '''\n",
        "    \n",
        "    feat = np.array([RFEAT[m][s, a] for m in range(len(RFEAT))])\n",
        "    \n",
        "    grad = feat * (feat.dot(wold) - r)    \n",
        "    wnew = wold - STEP * grad\n",
        "\n",
        "    return wnew"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-bLtdUSwJ76n"
      },
      "source": [
        "## Demonstration updater\n",
        "\n",
        "We perform maximum likelihood estimation using SGD, where \n",
        "\n",
        "$$\\mathbb{P}[a_t\\mid s_t]=\\frac{\\exp\\{\\eta Q^*_w(s,a)\\}}{\\sum_{a'}\\exp\\{\\eta_Q Q^*_w(s,a')\\}},$$\n",
        "\n",
        "where $Q^*_w$ is the optimal $Q$-function for the MDP with reward $w^\\top\\phi$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "id": "O-ltVWL7J76n"
      },
      "outputs": [],
      "source": [
        "def policy_feat(pol, feat):\n",
        "    '''\n",
        "    Compute the policy-weighted feature matrix. \n",
        "    \"pol\" can be a single action index or a full \n",
        "    |S| x |A| policy matrix.\n",
        "    \n",
        "    :param pol: int (action index) or np.array with shape (|S|, |A|) (policy)\n",
        "    :param feat: tuple (each element is a (|S|, |A|)-shaped np.array)\n",
        "    :return: np.array with shape (|S|, |nfeat|)\n",
        "    '''\n",
        "    \n",
        "    if isinstance(pol, int):\n",
        "        return np.hstack([feat[m][:, pol, None] for m in range(len(feat))])\n",
        "    \n",
        "    else:\n",
        "        return np.hstack([np.sum(feat[m] * pol, axis=1, keepdims=True) for m in range(len(feat))])\n",
        "    \n",
        "# --\n",
        "\n",
        "def policy_trans(pol, P):\n",
        "    '''\n",
        "    Compute the policy-weighted transition probabilities.\n",
        "    \n",
        "    :param pol: np.array with shape (|S|, |A|) (policy)\n",
        "    :param P: tuple (each element is a (|S|, |S|)-shaped np.array)\n",
        "    :return: np.array with shape (|S|, |S|)\n",
        "    '''\n",
        "    \n",
        "    return np.sum([np.diag(pol[:, a]).dot(P[a]) for a in range(len(A))], axis=0)\n",
        "\n",
        "# --\n",
        "\n",
        "def policy_rew(pol, r):\n",
        "    '''\n",
        "    Compute the policy-weighted reward.\n",
        "    \n",
        "    :param pol: np.array with shape (|S|, |A|) (policy)\n",
        "    :param r: np.array with shape (|S|, |A|) (reward)\n",
        "    :return: np.array with shape (|S|,)\n",
        "    '''\n",
        "    \n",
        "    return np.sum(pol * r, axis=1)\n",
        "\n",
        "# --\n",
        "\n",
        "def discounted_freq(pol, P):\n",
        "    '''\n",
        "    Compute the discounted visitation frequencies given a policy pol.\n",
        "    \n",
        "    :param pol: np.array with shape (|S|, |A|) (policy)\n",
        "    :param P: tuple (each element is a (|S|, |S|)-shaped np.array)\n",
        "    :return: np.array with shape (|S|, |S|)\n",
        "    '''\n",
        "    \n",
        "    Ppol = policy_trans(pol, P)\n",
        "    return (1 - gamma) * np.linalg.inv(np.eye(len(S)) - gamma * Ppol)\n",
        "\n",
        "# --\n",
        "    \n",
        "# Weight update\n",
        "def demo_update(wold, s, a, pol=None):\n",
        "    ''' \n",
        "    Updates the weights wold given the pair (s, a) and the weight wold. \n",
        "    The update is a SGD update on the negative log likelihood, assuming the\n",
        "    actions are selected from a Boltzmann distribution (pol, which can be\n",
        "    precomputed).\n",
        "    \n",
        "    :param wold: np.array with shape (|nfeat|,)\n",
        "    :param s: int (state index)\n",
        "    :param a: int (action index)\n",
        "    :param pol: np.array with shape (|S|, |A|) (policy)\n",
        "    :return: np.array with shape (nfeat,)\n",
        "    '''\n",
        "\n",
        "    if pol is None:\n",
        "        rhat = np.sum([RFEAT[m] * w[m] for m in range(len(RFEAT))], axis=0)\n",
        "        M = (S, A, P, rhat, gamma)\n",
        "        Q = vi(M)\n",
        "        pol = boltzmann(Q)\n",
        "\n",
        "    feat = policy_feat(greedy(Q), RFEAT) # (nS, nF)\n",
        "    disc = discounted_freq(greedy(Q), P) # (nS, nS)\n",
        "    \n",
        "    grad = np.zeros(wold.shape)\n",
        "\n",
        "    for a_idx in range(len(A)):\n",
        "              #              (nS, nF)                     (nS, nS) (nS, nS) (nS, nF)\n",
        "        grad += ETA_Q * (policy_feat(a_idx, RFEAT) + gamma * P[a_idx].dot(disc).dot(feat))[s, :] * (pol[s, a_idx] - int(a == a_idx))\n",
        "        \n",
        "    # Gradients are too small, so we push the step-size up\n",
        "    \n",
        "    wnew = wold - 2 * STEP * grad\n",
        "        \n",
        "    return wnew"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vQFCuICrJ76n",
        "outputId": "9661fb32-3b29-4d8d-9436-37260697b032"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "policy_feat(.): (12, 3)\n",
            "policy_feat(.): (12, 3)\n",
            "policy_trans(.): (12, 12)\n",
            "policy_rew(.): (12,)\n",
            "discounted_freq(.): (12, 12)\n"
          ]
        }
      ],
      "source": [
        "print('policy_feat(.):', policy_feat(pol, RFEAT).shape)\n",
        "print('policy_feat(.):', policy_feat(0, RFEAT).shape)\n",
        "print('policy_trans(.):', policy_trans(pol, P).shape)\n",
        "print('policy_rew(.):', policy_rew(pol, r).shape)\n",
        "print('discounted_freq(.):', discounted_freq(pol, P).shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r7k45CqBJ76o"
      },
      "source": [
        "## Contrastive examples updater\n",
        "\n",
        "We consider the following contrastive examples:\n",
        "\n",
        "- In state $s$ action $a$ is better than action $a'$ because it will eventually lead you through state $s'$ and that is good.\n",
        "- In state $s$ action $a$ is worse than action $a'$ because it will eventually lead you through state $s'$ and that is bad.\n",
        "\n",
        "We consider that, given the target reward $r$,\n",
        "\n",
        "- A state $s$ is good if $r^*(s)>0$ for the optimal policy $\\pi^*_r$ given that reward.\n",
        "- A state $s$ is bad if $r^*(s)<0$ for the optimal policy $\\pi^*_r$ given that reward.\n",
        "- An action $a$ is _better_ than an action $a'$ in state $s$ and leading to state $s'$ if the following two conditions are cumulatively met:\n",
        "\n",
        " - $Q^*(s,a)>Q^*(s,a')$ ($a$ is better than $a'$)\n",
        "    * The transition probabilities by first taking action $a$ and then following $\\pi^*$ lead to larger transition probability from $s$ to $s'$ than those same probabilities taking action $a'$.\n",
        "    \n",
        "- An action $a$ is _worse_ than an action $a'$ in state $s$ and leading to state $s'$ if the following two conditions are cumulatively met:\n",
        "\n",
        " - $Q^*(s,a)<Q^*(s,a')$ ($a$ is worse than $a'$)\n",
        "    * The transition probabilities by taking action $a$ in $s$ and following $\\pi^*$ elsewhere lead to larger transition probability from $s$ to $s'$ than those same probabilities taking action $a'$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "id": "MTFJKY7ZJ76o"
      },
      "outputs": [],
      "source": [
        "def sigmoid(z, eta):\n",
        "    ''' \n",
        "    Compute sigmoid function.\n",
        "    \n",
        "    :param z: float\n",
        "    :param eta: float\n",
        "    :return: float\n",
        "    '''\n",
        "    \n",
        "    return 1 / (1 + np.exp(-eta * z))\n",
        "\n",
        "# --"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mQ4Px7V8J76o"
      },
      "source": [
        "We perform maximum likelihood estimation using SGD, where \n",
        "\n",
        "$$\\mathbb{P}[(s,a,a',s',k)\\mid w]=\\sigma(kr_w(s'))\\sigma(k(Q(s,a)-Q(s,a')))\\mathbf{P}_{a}^\\infty(s'\\mid s)$$\n",
        "\n",
        "(see precise definition in the paper)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "id": "Uq4h5tnFJ76p"
      },
      "outputs": [],
      "source": [
        "# Weight update\n",
        "def explain_update(wold, s0, s1, a0, a1, k, Q=None):\n",
        "    '''\n",
        "    Updates the weights wold given the tuple (s0, s1, a0, a1, k) \n",
        "    and the weight wold. The update is a SGD update on the negative \n",
        "    log likelihood, assuming the distribution above.\n",
        "    '''\n",
        "    \n",
        "    # Compute policy\n",
        "    r = np.sum([RFEAT[m] * wold[m] for m in range(len(RFEAT))], axis=0)\n",
        "\n",
        "    if Q is None:\n",
        "        M = (S, A, P, r, gamma)\n",
        "        Q = vi(M)\n",
        "        \n",
        "    pol = boltzmann(Q)\n",
        "\n",
        "    paux = pol.copy()\n",
        "    paux[s0] = np.zeros(len(A))\n",
        "    paux[s0, a0] = 1.0\n",
        "        \n",
        "    # Component 1 - reward\n",
        "    fpol = policy_feat(greedy(Q), RFEAT)      # (nS, nF)\n",
        "    rpol = policy_rew(greedy(Q), r) # (nS,)\n",
        "\n",
        "    grad_R = ETA_R * k * fpol[s1, :] * (sigmoid(k * rpol[s1], ETA_R) - 1) # (nF,)\n",
        "    #print('grad_R:')\n",
        "    #print(grad_R)\n",
        "    \n",
        "    # Component 2 - Q-values\n",
        "    fa0 = policy_feat(a0, RFEAT)   # (nS, nF)\n",
        "    fa1 = policy_feat(a1, RFEAT)   # (nS, nF)\n",
        "    freq = discounted_freq(greedy(Q), P) # (nS, nS)\n",
        "\n",
        "    diff_Q = Q[s0, a0] - Q[s0, a1]       # scalar\n",
        "    diff_P = P[a0][s0, :] - P[a1][s0, :] # (nS,)\n",
        "    diff_F = fa0[s0] - fa1[s0]           # (nF,)\n",
        "    \n",
        "    grad_Q = ETA_Q * k * (diff_F + gamma * diff_P.dot(freq).dot(fpol)) * (sigmoid(k * diff_Q, ETA_Q) - 1) # (nF,)\n",
        "    #print('grad_Q:')\n",
        "    #print(grad_Q)\n",
        "\n",
        "    # Component 3 - Pam\n",
        "    Pam = discounted_freq(paux, P) # (nS, nS)\n",
        "    pf = Pam.dot(fpol)             # (nS, nF)\n",
        "    \n",
        "    grad_P = np.zeros(len(RFEAT))\n",
        "    \n",
        "    for k in range(len(RFEAT)):\n",
        "        # for a single k\n",
        "        grad_aux = np.zeros((len(S), len(S)))\n",
        "\n",
        "        for a in range(len(A)):\n",
        "            for ap in range(len(A)):\n",
        "                fa  = policy_feat(a, RFEAT)  # (nS, nF)\n",
        "                fap = policy_feat(ap, RFEAT) # (nS, nF)\n",
        "                    #  (nS,)            (nS, nS) (nS, nS) (nS,)\n",
        "                dT  = fap[:, k] + gamma * P[ap].dot(freq).dot(fpol[:, k]) # (nS,)\n",
        "                tmp = ETA_Q * np.diag(paux[:, a] * dT * (pol[:, ap] - int(a == ap))).dot(P[a]) # (nS, nS)\n",
        "\n",
        "                grad_aux += tmp\n",
        "\n",
        "        dPk = Pam.dot(grad_aux).dot(Pam)      # (nS, nS)\n",
        "        if np.isclose(Pam[s0, s1], 0) and not np.isclose(dPk[s0, s1], 0):\n",
        "            print('Whoa! Zero!')\n",
        "            print('dPk:', dPk[s0, s1])\n",
        "\n",
        "        elif np.isclose(Pam[s0, s1], 0):\n",
        "            grad_P[k] = 0.                        # scalar\n",
        "\n",
        "        else:\n",
        "            grad_P[k] = dPk[s0, s1] / Pam[s0, s1] # scalar\n",
        "        \n",
        "    #print('grad_P:')\n",
        "    #print(grad_P)\n",
        "    \n",
        "    grad = grad_R + grad_Q + grad_P\n",
        "    \n",
        "    wnew = wold - STEP * grad\n",
        "    \n",
        "    return wnew"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "EHi5ZyhkJ76p",
        "outputId": "6a61ed35-c776-4638-a16d-e1248e07bc70"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "P[a][s, :]: (12,)\n",
            "freq: (12, 12)\n",
            "fpol: (12, 3)\n",
            "P[a][s, :].dot(freq): (12,)\n",
            "P[a][s, :].dot(freq).dot(fpol): (3,)\n"
          ]
        }
      ],
      "source": [
        "freq = discounted_freq(pol, P)\n",
        "fpol = policy_feat(pol, RFEAT)\n",
        "print('P[a][s, :]:', P[0][0, :].shape)\n",
        "print('freq:', freq.shape)\n",
        "print('fpol:', fpol.shape)\n",
        "print('P[a][s, :].dot(freq):', P[0][0, :].dot(freq).shape)\n",
        "print('P[a][s, :].dot(freq).dot(fpol):', P[0][0, :].dot(freq).dot(fpol).shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2b_HimgRJ76p"
      },
      "source": [
        "## Experiments\n",
        "\n",
        "We now generate a random reward, and compare the performance of an RL updater, a demo updater and an contrastive examples updater in randomly selected sample rewards, demos and contrastive examples.\n",
        "\n",
        "We start by defining the samplers."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "id": "EaY9kXR1J76p"
      },
      "outputs": [],
      "source": [
        "def sample_reward(M):\n",
        "    ''' Samples a reward for MDP M. '''\n",
        "\n",
        "    S = M[0]\n",
        "    A = M[1]\n",
        "    R = M[3]\n",
        "    \n",
        "    # Sample a random state\n",
        "    sample_s = np.random.choice(len(S))\n",
        "    \n",
        "    # Sample a random action\n",
        "    sample_a = np.random.choice(len(A))\n",
        "    \n",
        "    # Sample a noisy reward\n",
        "    sample_r = R[sample_s, sample_a] + np.random.randn() / ETA_R\n",
        "    \n",
        "    return (sample_s, sample_a, sample_r)\n",
        "\n",
        "# -- \n",
        "\n",
        "def voice_reward(M, s, a, r):\n",
        "    ''' Converts a reward sample into a natural language string. '''\n",
        "    \n",
        "    return 'The reward in state %s when performing action %s is %.2f.' % (M[0][s], M[1][a], r)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 36
        },
        "id": "NlSSm3f0J76q",
        "outputId": "c093f169-d179-48a6-b9f5-ffddecfde683"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'The reward in state 1 when performing action Down is -0.38.'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 16
        }
      ],
      "source": [
        "voice_reward(M, *sample_reward(M))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "id": "YJNptXPOJ76q"
      },
      "outputs": [],
      "source": [
        "def sample_demo(M):\n",
        "    ''' Samples a demonstration for MDP M. '''\n",
        "\n",
        "    S = M[0]\n",
        "    A = M[1]\n",
        "    \n",
        "    # Sample a random state\n",
        "    sample_s = np.random.choice(len(S))\n",
        "\n",
        "    # Sample a noisy action\n",
        "    Q = vi(M)\n",
        "    pol = boltzmann(Q)[sample_s, :]\n",
        "    sample_a = np.random.choice(len(A), p=pol)\n",
        "    \n",
        "    return (sample_s, sample_a)\n",
        "\n",
        "# --\n",
        "\n",
        "def voice_demo(M, s, a):\n",
        "    ''' Converts a demonstration sample into a natural language string. '''\n",
        "\n",
        "    return 'In state %s you should perform action %s.' % (M[0][s], M[1][a])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 54
        },
        "id": "6G6YXh88J76q",
        "outputId": "c5b02127-e0dd-4d33-ba1c-9ce9b8ccfa95"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'In state 10 you should perform action Down.'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 18
        }
      ],
      "source": [
        "Q = vi(M)\n",
        "pol = boltzmann(Q)\n",
        "print(pol.sum(axis=1))\n",
        "voice_demo(M, *sample_demo(M))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "id": "s_RXa80KJ76q"
      },
      "outputs": [],
      "source": [
        "def sample_explain(M, echo=False):\n",
        "    ''' Samples an contrastive examples for MDP M. '''\n",
        "\n",
        "    S = M[0]\n",
        "    A = M[1]\n",
        "    P = M[2]\n",
        "    R = M[3]\n",
        "    gamma = M[4]\n",
        "    \n",
        "    Q = vi(M)\n",
        "    pol = boltzmann(Q)\n",
        "    rpol = policy_rew(pol, R)\n",
        "\n",
        "    sample_complete = False\n",
        "    \n",
        "    while not sample_complete:\n",
        "        \n",
        "        # This ensures that the search for next state doesn't go on forever\n",
        "        stop_sample = False\n",
        "        \n",
        "        # Sample a random initial state\n",
        "        sample_s0 = np.random.choice(len(S))\n",
        "        \n",
        "        if echo:\n",
        "            print('Sampled state %s' % S[sample_s0])\n",
        "\n",
        "        # Sample a random initial action\n",
        "        sample_a0 = np.random.choice(len(A))\n",
        "        \n",
        "        if echo:\n",
        "            print('Sampled action %s' % A[sample_a0])\n",
        "\n",
        "        paux = pol.copy()\n",
        "        paux[sample_s0] = np.zeros(len(A))\n",
        "        paux[sample_s0, sample_a0] = 1.0\n",
        "        \n",
        "        Pam = discounted_freq(paux, P)\n",
        "\n",
        "        # Sample a random second action (different from the first)\n",
        "        quit = False\n",
        "        while not quit:\n",
        "            sample_a1 = np.random.choice(len(A))\n",
        "            quit = sample_a1 != sample_a0\n",
        "\n",
        "        if echo:\n",
        "            print('Sampled another action: %s' % A[sample_a1])\n",
        "\n",
        "        # Sample a random via state with non-zero reward (to fit the \"good\" or \"bad\" description)\n",
        "\n",
        "        quit = False\n",
        "        while not quit:\n",
        "            paux = np.maximum(rpol * (Q[sample_s0, sample_a0] - Q[sample_s0, sample_a1]) * Pam[sample_s0, :], 0)\n",
        "\n",
        "            # Ensure that there is a next state to sample\n",
        "            if np.isclose(np.sum(paux), 0):\n",
        "                if echo:\n",
        "                    print('Zero support distribution. Restarting sampling.')\n",
        "                stop_sample = True\n",
        "                quit = True\n",
        "            else:\n",
        "                paux = paux / paux.sum()\n",
        "                if np.isclose(paux[sample_s0], 1.):\n",
        "                    if echo:\n",
        "                        print('Departing from sink. Restarting sampling.')\n",
        "                    stop_sample = True\n",
        "                    quit = True\n",
        "                else:\n",
        "                    if echo:\n",
        "                        print(np.round(paux, 3))\n",
        "                    sample_s1 = np.random.choice(len(S), p=paux)\n",
        "                    quit = (sample_s1 != sample_s0)\n",
        "\n",
        "        if not stop_sample:\n",
        "            if echo:\n",
        "                print('Sampled another state: %s' % S[sample_s1])\n",
        "\n",
        "            # Check which one is better\n",
        "            if Q[sample_s0, sample_a0] > Q[sample_s0, sample_a1]:\n",
        "                k = 1\n",
        "            else:\n",
        "                k = -1\n",
        "                \n",
        "            sample_complete = True\n",
        "            \n",
        "    if echo:\n",
        "        print('Q(%s, %s) =' % (S[sample_s0], A[sample_a0]), Q[sample_s0, sample_a0])\n",
        "        print('Q(%s, %s) =' % (S[sample_s0], A[sample_a1]), Q[sample_s0, sample_a1])\n",
        "        print('reward(%s) =' % S[sample_s1], np.sum(R * pol, axis=1)[sample_s1])\n",
        "        print('P_m[%s](%s, %s) =' % (A[sample_a0], S[sample_s0], S[sample_s1]), Pam[sample_s0, sample_s1])\n",
        "        \n",
        "    return (sample_s0, sample_s1, sample_a0, sample_a1, k)\n",
        "\n",
        "# --\n",
        "        \n",
        "def voice_explain(M, s0, s1, a0, a1, k):\n",
        "    ''' Converts an exmplanation sample into a natural language string. '''\n",
        "\n",
        "    REL1 = ('better', 'worse')\n",
        "    REL2 = ('good', 'bad')\n",
        "    \n",
        "    k = 1 - (k + 1) // 2\n",
        "    \n",
        "    return 'In state %s, action %s is %s than %s because it may eventually lead you through state %s and that is %s.' % (M[0][s0], M[1][a0], REL1[k], M[1][a1], M[0][s1], REL2[k])\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 346
        },
        "id": "yQLwYih7J76r",
        "scrolled": false,
        "outputId": "b7a3e06a-f4bc-46e6-d59e-95d0a2d29954"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Sampled state 2\n",
            "Sampled action Down\n",
            "Sampled another action: No-Op\n",
            "Zero support distribution. Restarting sampling.\n",
            "Sampled state 7\n",
            "Sampled action Left\n",
            "Sampled another action: Right\n",
            "Zero support distribution. Restarting sampling.\n",
            "Sampled state 6\n",
            "Sampled action Up\n",
            "Sampled another action: Left\n",
            "[ 0.  0.  1.  0.  0.  0.  0.  0.  0. -0. -0. -0.]\n",
            "Sampled another state: 3\n",
            "Q(6, Up) = 31.52620995857553\n",
            "Q(6, Left) = 31.21094785771591\n",
            "reward(3) = 0.33148853482453977\n",
            "P_m[Up](6, 3) = 0.7787574779843445\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'In state 6, action Up is better than Left because it may eventually lead you through state 3 and that is good.'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 20
        }
      ],
      "source": [
        "voice_explain(M, *sample_explain(M, echo=True))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E5_t25yaJ76r"
      },
      "source": [
        "We now run a comparative study. Every 10 steps, for each of the three approaches, we:\n",
        "\n",
        "* Estimate the reward;\n",
        "* Compute the associated \"optimal\" policy;\n",
        "* Evaluate that policy in the correct MDP.\n",
        "\n",
        "The next function is used in the last step."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "id": "1K4B1KnTJ76r"
      },
      "outputs": [],
      "source": [
        "def evaluate_reward(M, w):\n",
        "    ''' Computes the policy from weights w in the MDP M.'''\n",
        "\n",
        "    S = M[0]\n",
        "    A = M[1]\n",
        "    P = M[2]\n",
        "    R = M[3]\n",
        "    gamma = M[4]\n",
        "    \n",
        "    Raux = np.sum([w[i] * RFEAT[i] for i in range(len(RFEAT))], axis=0)\n",
        "    \n",
        "    Maux = (S, A, P, Raux, gamma)\n",
        "    pol  = greedy(vi(Maux))\n",
        "    \n",
        "    Ppol = policy_trans(pol, P)\n",
        "    Rpol = policy_rew(pol, R)\n",
        "    \n",
        "    Vaux = np.linalg.inv(np.eye(len(S)) - gamma * Ppol).dot(Rpol)\n",
        "\n",
        "    # Compute average value\n",
        "    return Vaux.sum() / len(S)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yzbntx4sJ76r"
      },
      "source": [
        "We can now run our comparative study."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 49,
          "referenced_widgets": [
            "41b50f7bfd3b480a952593c96da0980b",
            "1131e727887c4030ad451a44f11870e9",
            "581f0a2815d44294a7084bfb2ed3ba95",
            "1764f7a1af18499c8ec767b06914f439",
            "021806461ab44c5abdd4629d6d268199",
            "e633f0851f984feea42c0f676462a47a",
            "e79702d34d7b426eb7f6612532abb0e1",
            "86248af124514eb995033fcd604c8fab",
            "39a5ec80a5234a7dbcdf7615df616518",
            "faade32c915f4120895baada845e1c90",
            "3274dd63da9747959237c72b1ab8a738"
          ]
        },
        "id": "KcBq_eJKJ76s",
        "scrolled": false,
        "outputId": "7ceca9b7-217e-4ef7-c5e7-769ce6b4add1"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "  0%|          | 0/30 [00:00<?, ?it/s]"
            ],
            "application/vnd.jupyter.widget-view+json": {
              "version_major": 2,
              "version_minor": 0,
              "model_id": "41b50f7bfd3b480a952593c96da0980b"
            }
          },
          "metadata": {}
        }
      ],
      "source": [
        "# Simulation constants\n",
        "NRUNS = 30\n",
        "NSTEPS = 15\n",
        "\n",
        "# For performance statistics\n",
        "perf_rl   = np.zeros((NRUNS, NSTEPS))\n",
        "perf_demo = np.zeros((NRUNS, NSTEPS))\n",
        "perf_expl = np.zeros((NRUNS, NSTEPS))\n",
        "perf_ref  = np.zeros((NRUNS, NSTEPS))\n",
        "\n",
        "# For reproducibility\n",
        "np.random.seed(40)\n",
        "\n",
        "samples_expl = []\n",
        "samples_demos = []\n",
        "samples_rl = []\n",
        "\n",
        "# Run sims\n",
        "for run in trange(NRUNS):\n",
        "    \n",
        "    # Initialize a random reward\n",
        "    wopt = np.random.randn(len(RFEAT))\n",
        "    ropt = np.sum([wopt[i] * RFEAT[i] for i in range(len(RFEAT))], axis=0)\n",
        "    \n",
        "    #print('Reward:')\n",
        "    #print(np.round(ropt,2))\n",
        "                  \n",
        "    # This is the \"target\" MDP\n",
        "    Mopt = (S, A, P, ropt, gamma)\n",
        "    \n",
        "    # Optimal policy\n",
        "    popt = greedy(vi(Mopt))\n",
        "    \n",
        "    #print('Optimal policy:')\n",
        "    #print(np.round(popt, 2))\n",
        "    \n",
        "    # Weights for 3 approaches\n",
        "    w_rl   = np.zeros(len(RFEAT))\n",
        "    w_demo = np.zeros(len(RFEAT))\n",
        "    w_expl = np.zeros(len(RFEAT))\n",
        "    \n",
        "    # Run test\n",
        "    for t in range(NSTEPS):\n",
        "        #print('\\n[Step %i]' % (t + 1))\n",
        "\n",
        "        # Sample a reward and update\n",
        "        #print('Sampling reward... ')\n",
        "        (s, a, r) = sample_reward(Mopt)\n",
        "        #print('Updating weights for reward... ')\n",
        "        w_rl = rl_update(w_rl, s, a, r)\n",
        "        #print('Weights:')\n",
        "        #print(w_rl)\n",
        "\n",
        "        \n",
        "        # Evaluate\n",
        "        #print('Evaluating reward...')\n",
        "        perf_rl[run, t] = evaluate_reward(Mopt, w_rl)\n",
        "       \n",
        "        #print(voice_reward(Mopt, s, a, r))\n",
        "        \n",
        "        # Sample a demo and update\n",
        "        #print('\\nSampling demo... ')\n",
        "        (s, a) = sample_demo(Mopt)\n",
        "        #print('Updating weights for demo... ')\n",
        "        w_demo = demo_update(w_demo, s, a)\n",
        "        #print('Weights:')\n",
        "        #print(w_demo)\n",
        "\n",
        "        \n",
        "        # Evaluate\n",
        "        #print('Evaluating reward...')\n",
        "        perf_demo[run, t] = evaluate_reward(Mopt, w_demo)\n",
        "        \n",
        "        #print(voice_demo(Mopt, s, a))\n",
        "        \n",
        "        # Sample an contrastive examples and update\n",
        "        #print('\\nSampling contrastive examples... ')\n",
        "        (s0, s1, a0, a1, k) = sample_explain(Mopt)\n",
        "        \n",
        "        #print('Updating weights for contrastive examples... ')\n",
        "        w_expl = explain_update(w_expl, s0, s1, a0, a1, k)\n",
        "        #print('Weights:')\n",
        "        #print(w_expl)\n",
        "        \n",
        "        # Evaluate\n",
        "        #print('Evaluating reward...')\n",
        "        perf_expl[run, t] = evaluate_reward(Mopt, w_expl)\n",
        "        \n",
        "        #print(voice_explain(M, s0, s1, a0, a1, k))\n",
        "        \n",
        "        perf_ref[run, t] = evaluate_reward(Mopt, wopt)\n",
        "        \n",
        "        #print('Done. Step complete.')\n",
        "        \n",
        "# Compute statistics\n",
        "rl_mean = np.mean(perf_rl, axis=0)\n",
        "rl_std = np.std(perf_rl, axis=0) / np.sqrt(NRUNS)\n",
        "\n",
        "demo_mean = np.mean(perf_demo, axis=0)\n",
        "demo_std = np.std(perf_demo, axis=0) / np.sqrt(NRUNS)\n",
        "\n",
        "expl_mean = np.mean(perf_expl, axis=0)\n",
        "expl_std = np.std(perf_expl, axis=0) / np.sqrt(NRUNS)\n",
        "\n",
        "ref_mean = np.mean(perf_ref, axis=0)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "rl_min = np.min(rl_mean, axis=0)\n",
        "rl_range = np.max(rl_mean, axis=0) - np.min(rl_mean, axis=0)\n",
        "rl_norm = (rl_mean - rl_min)/rl_range\n",
        "rl_normstd = np.std(rl_norm, axis=0) / np.sqrt(NRUNS)\n",
        "\n",
        "demo_min = np.min(demo_mean, axis=0)\n",
        "demo_range = np.max(demo_mean, axis=0) - np.min(demo_mean, axis=0)\n",
        "demo_norm = (demo_mean - demo_min)/demo_range\n",
        "demo_normstd = np.std(demo_norm, axis=0) / np.sqrt(NRUNS)\n",
        "\n",
        "expl_min = np.min(expl_mean, axis=0)\n",
        "expl_range = np.max(expl_mean, axis=0) - np.min(expl_mean, axis=0)\n",
        "expl_norm = (expl_mean - expl_min)/expl_range\n",
        "expl_normstd = np.std(expl_norm, axis=0) / np.sqrt(NRUNS)\n",
        "\n",
        "\n",
        "ref_min = np.min(ref_mean)\n",
        "ref_range = np.max(ref_mean) - np.min(ref_mean)\n",
        "ref_norm = (ref_mean - ref_min)/ref_range\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SAFud7XsSAhF",
        "outputId": "ca343ac9-f22b-4a97-d8db-1c51f82a8983"
      },
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-23-a6cb8b7f165e>:19: RuntimeWarning: invalid value encountered in true_divide\n",
            "  ref_norm = (ref_mean - ref_min)/ref_range\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 470
        },
        "id": "d3cHTDXIJ76s",
        "outputId": "391aa1b4-bb71-40b0-ac2b-e7f870a46679"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 400x500 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAHFCAYAAAAUkEVeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACTz0lEQVR4nOydd5yTVdbHf096nWR6oQ1NulSlubICguKiuCgWVBBFX0WliCKLYBdlFRELrmXB3lax4IqLqCDSm4IgoNKZAsxMMmlPve8fd5KZTIFJJnVyv5oPk5v73OdkJrnnnnPPPYcjhBAwGAwGIyVRxVsABoPBYMQPpgQYDAYjhWFKgMFgMFIYpgQYDAYjhWFKgMFgMFIYpgQYDAYjhWFKgMFgMFIYpgQYDAYjhWFKgMFgMFIYpgQYDAYjhYmrEli7di1Gjx6NgoICcByHzz77LPCaKIqYNWsWevToAbPZjIKCAtx00004ceJE0BhlZWUYP3480tLSYLfbccstt8DlcsX4nTAYDEZyElcl4Ha70bNnT7z00kt1XvN4PNi+fTvmzp2L7du349NPP8W+fftw+eWXB/UbP348fv31V6xatQorVqzA2rVrcdttt8XqLTAYDEZSwyVKAjmO47B8+XKMGTOmwT5btmzB+eefj8OHD6N169bYu3cvunbtii1btqBfv34AgJUrV2LUqFE4duwYCgoKYiQ9g8FgJCeaeAsQCg6HAxzHwW63AwA2bNgAu90eUAAAMHz4cKhUKmzatAlXXnllvePwPA+e5wPPFUVBWVkZMjMzwXFcVN8Dg8FgxAJCCCorK1FQUACVqmGnT9IoAZ/Ph1mzZuG6665DWloaAKC4uBg5OTlB/TQaDTIyMlBcXNzgWPPnz8cjjzwSVXkZDAYjETh69ChatmzZ4OtJoQREUcS4ceNACMGSJUuaPN7s2bMxY8aMwHOHw4HWrVvj6NGjAQXDYDAYyYzT6USrVq1gtVrP2C/hlYBfARw+fBjfffdd0CSdl5eH0tLSoP6SJKGsrAx5eXkNjqnX66HX6+u0p6WlMSXAYDCaFWdzcSf0OQG/Ajhw4AC+/fZbZGZmBr0+cOBAVFRUYNu2bYG27777DoqioH///rEWl8FgMJKOuFoCLpcLv//+e+D5wYMHsXPnTmRkZCA/Px9XXXUVtm/fjhUrVkCW5YCfPyMjAzqdDl26dMEll1yCyZMn45VXXoEoirjrrrtw7bXXssggBoPBaARxDRH94YcfcNFFF9VpnzBhAh5++GG0bdu23uu+//57/PWvfwVAD4vddddd+PLLL6FSqTB27FgsXrwYFoul0XI4nU7YbDY4HA7mDmIwGM2Cxs5rCXNOIJ4wJcBgMJobjZ3XEnpPgMFgMBjRhSkBBoPBSGGYEmAwGIwUhikBBoPBSGGYEmAwGIwUJuFPDDPCR5YBSap+Hk4cWEPXNHQIMdz22vep775n61P7+ZkOSob7WkOynan9bK9FS85IEO57Cuf9htqeSmg0gFodpbGjMywj1sgyIIr0IQiAx0P/FcXgfqEqgkgogcb0jcSkHwkl0JgJJxGUQCwnxli930R4r4lKVhZ9RAOmBJIQSWp4wpdl+qXRaunqwWQKbWI+GxwX+qTQmAm+PpnO9vxMfcKZoEJRkOFMWPW9FukJNlrE6/2yU0yA0wkoSvTGZ0ogwfFP9qII8Dzg9VZP+IoSPOFbLNEzGWuSDCuzeLtRGkuqrX5T7f1Ggmj/bpgSSDA8nuoJ3/+zKNIVUc0J32gEzlAngsFgMBoFUwIJhNsNHDtGJ321mk74Wm3DLh0Gg8FoKkwJJBAuF/XpZ2TEWxIGg5EqMIdCgiBJQGUlYDDEWxIGg5FKMCWQIHg8dB+gnoJnDAaDETWYEkgQKivpPgDz/TMYjFjClEACwPN0U5i5ghgMRqxhSiAB8Mf+63TxloTBYKQaTAnEGUKAigqmABgMRnxgSiDO+HzUEjAa4y0Jg8FIRZgSiDNuNz0bEIt0DwwGg1EbpgTiiCwDDgezAhgMRvxgSiCOeL3UHcTOBjAYjHjBlEAccbvpuQCWCI7BYMQLNv3ECVGkecLZ2QAGgxFPmBKIE14vSxPBYDDiD1MCccLhoHUBWJoIBoMRT5gSiAP+gjEsKojBYMQbpgTigL9imFYbb0kYDEaqw5RAjCGEuoLYXgCDwUgEmBKIMV4vtQRYVBCDwUgEmBKIMR4PtQZYmggGg5EIMCUQQ/xpIpgVwGAwEgWmBGKIx0PTRDAlwGAwEgWmBGKIy0XPBbCzAQwGI1FgSiBGCAJVAuxsAIPBSCSYEogRLE0Eg8FIRJgSiAGE0GRx7HAYg8FINJgSiAE8T9NGM1cQg8FINJgSiAFuNyBJNGEcI8WRZboqYDASBDYtRRlFYXUDGFUQApw+TT8QVit9GI0sXIwRV5gSiDJeL32kpcVbEkbcqawEystpdIDDQR9mM2CzASYTKzHHiAtMCUQZt5v+y77fKY7PB5w6Beh01Cw0GKhryOulysFkosrAbGZ+Q0ZMYZ+2KCJJ1PJnYaEpjixTN5AkBZuEajVgsVCfoc8HFBXRD4tfGbAPDiMGMCUQRbxe+t1OT4+3JIy44nDQ1X5DPkGViloCRiPdNC4tpfHEaWlUSRgMbN+AETWYEogiTidd7LHvbwrjdlMroDE+f46rdhUJAlBWBlRUUEWQlkaVBPMrMiIMUwJRQhDY2YCURxCAkyfp5K7ThXatTkcfskzzjdTcNzCZWC5yRsRgSiBKeDx0DrBY4i0JIy4QQlfyPN+00DC1moaS+jeRXS5qKdjtdN+AHUNnNBGmBKKAP01EqIs/RjPC4aCuHKs1Mv5A/yYyIXSjqbiYbhynpVHLoPY96rtnKG0sQillYH/pJuKuigE1mUzgqr5QTqeAkydFpKVpAFRHeHg8tK/BYISqyrcriiJEUYBKpYahxomyUPp6vR4QQqDXG6CuchNIkgSPh4fbrYJeb6wxrgeKQqDT0b6EAKIoQRB4cBztS4h/rvGCEAUajR5qNf2oyLIMSfJBpVLBaDRCpaJuakGgffV6PXQ6DVQqgBAZougDx3Ewm02BvjzvAyEydDodtFUrWVGU4fX6QAgHvd4ERaFBM16vD6IoQ63WQa3WQlEASVLg9XqhKIBebw709fl4SJIErVYLrVZXJYMCUfSC4wCTyRxwqUsSD0WRoNNpodPpqtoJvF5P1T6tCWo1B44DeF6AIIhQq7XQaHQgBFAUAo/HA0IAnc4EQujfnucFiJWV4E6dhkZnAXHoQBTA43NDUTjotPTvqVIRSJIAWRGh0Whg0OvBgf5+fLwbHEdg0BuhVqugUgGSJEKRBag1ahj0JnCcCZxLAF90HACg11X/7UVJhChLUKvU0OuqP39enxcEBHqtvvpzosgQRAEqTgWD3kD/9hwHryRBVmugs9ig0uuhcGqIhAMvigCCPyderxeyrECtpp8TQgBJkuHz+QBwMBjo357j6n5OOK76c6JScVXfI/9nin5OtFod9Hpt4O/J896q75w58N543gdJon39nylFUeDz0b5GY3VfQaCfE42G/u3p+6B/+9rj1u4LBPc1Gqu/94IgQJJEqNUa6HR1v/f19dVo6u9b83t/+rQIozGKFh9hEIfDQQAQh8MR8rUACABSWloaaPvHPx4nAMj1199Kjh8ngYfRaCIAyMaNBwNtDz/8HAFArrzy+qC+GRlZBAD57rvdgbYFC14lAMjIkVcE9W3Zsg0BQFas2Ey+/56QRx4hpFu3d6pkG171VfU/ula1f1+jbXlV26BafftVta+o0fa/qraetfoOqWr/qEbbuqq2DrX6jqpqX0pUKn/bjqq2glp9r6pqf7FG2/6qNlutvhOq2hfUaDtW1aap1ffOqvaHarSVB/6egFCjfWZV28wabUKNvuU12h+qaruz1v00Ve3HarQtqGqbUKuvrap9f422F6varqrVt6CqfUeNtqVVbaNq9e1Q1b6uRttHVW1DavXtWdX+vxptK6ra+tXqO6iqfXmNtu+r2rrW6ju8qv2dGm2bq9ra1Op7RVX7qzXadle1ZdXqe31V+3M12g5WtZlq9b21qv3xGm2lNf6eNftOrWr7R402V42+rhrt/6hqm1prDH/f0hptj1e13Vqrr6mq/WCNtufIDTdEb15jlkCE8Yd8x4rSUrr/AAATJtBAlFDwewMIof9qNNVtklQ1fWmqvQOyDIgi/VmtRtWqOHz5Q73Wv0qUZfrcaKRyqFTUZS6K1GWu16NqRVr9+/EbT9SaaJrcNVGr/VFgBJJEIMuARkNg0CtQVRURclTSvvY0BSqVQlfQPgIfX91XUVDVTvtqNQSEEBAAskwiI2wT4DgCQgAOgFargAPAqejel0IAg55ApyPgOECSCdxu+nux2WpMny76u9frqz8/skzHAOjnzP97IPF/yykBV7WaTWmcTidsNhscDgfSQtzEq+0OcruBAwcEmEwitFoN9PrIuoMcDgHbtqmxcaMBP/4I7N0LAB7QxYYBBoMa558PXHCBhP79eXTqpILJZAxUNKNmbF3Xkd8dZKwRzuT1UtNdp9NDo6l2B/G8r05fj4e6BPyuIzrRVrt4DAZTDbcNNd39Lh46GVCXgFpNXUf+0FpRpC4Bg0EHnU4Ljgs284NdAjxkOdh0b2xfqsyqXTx6vQkAV/U+BMiyCJ1OC70+2HUE1DDzKyshHDkMSauF2mAIcsV4qvoa9IbA31MQBUiSBLVafda+oihClETqstMZAsrX7a2SV6cPdgeJItQqFfT6Gm5DX92+kiRVuYM4GAz078lxgMBT15FBr4dGo6YTe1VfjuNgNBipALIMr9sFIgjQqdXQaLSAWg2Z48CDgDMYYLTZ6SY2x53lM8XBaDTV+Pz5IMsyNBr6OaFKnboC/e696r8n7VvbHVS/66ju56Tm3/NMfan76uzuIIOhce4gtfrsc8TJkyIKCrTIz0dINHZeY0oATVMCtSkpodkBInVATJaBXbuAtWvpY9u26lWTnx49gAsvBC64ADj/fJasLi7wPD3xqyh0ozaV8ZuLklRtsmVlARkZ7NBMGFRU0F9fTk5o1zV2XovryZO1a9di9OjRKCgoAMdx+Oyzz4JeJ4Rg3rx5yM/Ph9FoxPDhw3HgwIGgPmVlZRg/fjzS0tJgt9txyy23wOVyxfBdVCNJNJy7qWcDDh8G3n4buO024NxzgcsuA55+GtiwgSqAFi2A664DXn6ZKoiVK4F//IMqAqYA4oCiVIeDproCAKifx2Cg0Uz+cw0nT9LkeWzNmXDEdU/A7XajZ8+emDRpEv7+97/XeX3BggVYvHgx3nzzTbRt2xZz587FyJEjsWfPnoDrZPz48SgqKsKqVasgiiJuvvlm3HbbbXjvvfdi/Xbg8YSfJkJR6ET/5ZdUCdTEagUGDaKT/F/+ArRrxxZUCUVFBQ0JtVrjLUliotVWKwKAWgSMhCFh3EEcx2H58uUYM2YMAGoFFBQU4N5778XMmTMBAA6HA7m5uVi2bBmuvfZa7N27F127dsWWLVvQr18/AMDKlSsxatQoHDt2DAUFBY26d6TcQcePnzlFzJlYtw645hr6s0YD9OlTPen36sXCthMWjwc4cYJOdCzh25nhefrIzaUWAqNRNGt30Jk4ePAgiouLMXz48ECbzWZD//79sWHDBgDAhg0bYLfbAwoAAIYPHw6VSoVNmzY1ODbP83A6nUGPpsLz1Yc5w+G//6X/XnYZ8OuvwPLlwPTpQL9+TAE0Cn9YltNJJ+ZYrG0kqTociymAs6PX0xOUpaX078RICBJ2eikuLgYA5ObmBrXn5uYGXisuLkZOLfWo0WiQkZER6FMf8+fPxyOPPBJReT0euhcWjkdAUahfH6DWAEs10QgIob9wQaCTv9tdvRmpUtE/hM0WvcpdhNB9ALebrWpDwb9KKi2lfxfmQos7CWsJRJPZs2fD4XAEHkePHm3SeIRQl3C4aSK2baNRRVYrjfBhNIAo0km3rAw4dgw4coT+W1ZG/whGY3VOHZeLvlZcHB3LwF8lLFJpIVIJg4Eq6tLS6qpLjLiRsJZAXl4eAKCkpAT5NQJkS0pK0KtXr0Cf0tLSoOskSUJZWVng+vrQ6/VBsblNxeejB5XCXcF//TX9d/hw5lUIwn+KiOfpZMHz1SfVtFo6mdTnK6uZdK2ykj4iaRnUrBLGsnmGh8lE/6YlJUBeHouqiiMJawm0bdsWeXl5WL16daDN6XRi06ZNGDhwIABg4MCBqKiowLZt2wJ9vvvuOyiKgv79+8dMVrebzjfhzAeEVCuBUaMiK1fSUdOvX1ICHD1KH8XFVMtqNHTX3R92eLbNErWa9jebqSKIhGVQs0oYyxPeNMxm+jcvKaF/X0ZciKsl4HK58PvvvweeHzx4EDt37kRGRgZat26NadOm4fHHH0fHjh0DIaIFBQWBCKIuXbrgkksuweTJk/HKK69AFEXcdddduPbaaxsdGdRUZJm6gsKdD379lXo1DAbgr3+NqGjJgSjSid+/sy6KqMq7QFf7FkvTC6n4lYHfMnC5qmPYQ7UMzlYljBEaFgv9ffotAnbQJebEVQls3boVF110UeD5jBkzAAATJkzAsmXLcP/998PtduO2225DRUUFLrjgAqxcuTIovcK7776Lu+66C8OGDYNKpcLYsWOxePHimL0HfwlJuz286/1RQUOHpqBFTAj1C1dW0olep4tuwZSGlIHd3rgSjqFUCWM0HquVWn+lpTR8lPlEg+BlH7wiARAdyzNhzgnEk6acEygqovuD4SqBv/4VOHAAeOEFoJ7zcs0bQag2g+IRByvLdGL3R6nYbA0rA0Gg5wEkiYVvRQNCqo/b5+ayYhxVCLKA344Wo1WBEd0Ks0O6NunPCSQDokgXk+G6gg4coA+tlm4KpxyCUO36iQd+y8BkoivRY8eq/dM110b+cFCep35sRuTxK2KPh1oE/gCAFEZWZJx0n4RbrIzqfRI2OigZ4Hk6j4U7L/hdQX/5S4q6mHk+3hJQ/MpAkqp9/jWjiSJdJYxRPxxH/w4OB/05NzdlT0oSQlDmLUMlXwmNOrolRFPzN5wgpHRUECHUFZNIZr9GQyf+2srA7aZuIhYOGn38isDppD/n5KSkIqjwVaDMWwaLzgI3F90CJcwdFCeOHKEZQFUqYMSIeEsTBwSBPhKxULpfGRiNVBEQwqJWYolKVa0ITp2qTkedIrgFN057TsOgMUCjjr4CTD0VmyD4rYABA4DMzPjKEhf8KR4SeZWn0bC0BvHCn/qjooJaBNnZKRGR5ZN8OOk+CY7joNfEJkqq+f9WExT/fkBKuoIAGlebAl9qRhPwn/wuL6ehuZGqB5qgSIqEU+5TEGQBZl3sAhDYtzAOlJQAW7fSny+5JL6yxAX/fkAiuoIYiYVfEZw+XZ0jqhmiEAWnPafhEl2w6mtYn5IMldtbXXg6CiSwLd588buC+vRByHVDmwWCQN1BLO0CozGo1fRsxqlT1HpMT292UVoV3gqU+8ph1VnBEQLOx4PzeMBVuqEqdoLLiN6CiSmBOJDSUUFA9fkAFm3DaCwaDY3F9lcns9ubjTuxkq/EKc8pmBQNNJVucJUucD4fQAiIXgeoo6vwmBKIMWVltFYwAFx6aXxliRs+X7NbyTFigFZLrceTJ2nkkNlMo7Z0OvpIws+U1+vC6VOHoHd7YBQIOFkG0WpBzCZA5V8kscNizYpVq+giuGtXoLAw3tLEAUVJvPMBjOTBn75bFOmGsaJUJxs0m6urlyWyUqjKlitWOlBefBDEWwmzKQ1ErwfRxH6fjCmBGJPyUUGCwNIwM5qGWk0f/rMb/roTZWXVSsGfjNBgoAoiEZSCzxdIla54PSj3nIYHPKxZeSBc/FxbTAnEkMpKYO1a+nPKKwG2H8CIFGo1XVT4FxZ+pXD6NI0m0miohWAyVVsKWm1slIIk0VxUlZU0L5IkgWi1qNDIcOoJzLoMcHFUAABTAjHlu+/oZ7N9e+Ccc+ItTZzw+ZgCYESX2kpBkqj76NSp6tf1euo+8ruO/J/JSIWgCgJ1e7pcNEeWRhPIluviXSj3OGDUmqDm4v9dYEoghnz1Ff330kvjb5nGBUWhqyF2PoARSzQa+qipFASBKgVFoZ9Hf6RRpJSAKNIvuV5PU5BUfeG9ohenvaehU2mhVSXG94ApgRjh9VJLAAAuuyy+ssQN//mAlKuew0go/ErBjySd/TRyqKs2g6FOCKsgizjlOQUCAr0mcXJRMSUQI9asoYqgZUugR494SxMn2PkARiISg/xVsiLjtPc0eFmAVZdY+aiax2mLJMAfFZSyriCAakGmABgpBiEEFb4KuHkXLDoLuASbAJgSiAGCQM8HACkcFaQoVAmw/QBGiuHgnSj3VcCkM0MV50ig+kg8iZoh69fTA445OUC/fvGWJk7wPN0PYIfEGCmER/Sg3FsOvVoHjSoxve9MCcQAvyto5Mhmk+4kdPz7ASn7C2CkGj6Zx2nP6araAImzEVwb9o2MMrIMrFxJf05ZVxBAXUGJXECGwYggkiKjzFMGQRZh0iZ2NBxTAlFmyxZ6cNFuBwYOjLc0cUKW6SExth/ASAEUoqDcWwaP6IZFZ4m3OGeFLc2ijN8VdPHFKTwH+usJWxL/C8FILXhZgIuvhKCIDfbhEFo0DwGBW3DBnICRQPXBlEAUIYQljANAFYCisP0ARsIgKTIq+Uo4eSdERYAmwqd3TVpzQqSEaAxMCUSRnTuBoiKaouTCC+MtTRzxeNh+ACMhUIgCt+BBha8cvMTDoDXAqLXFW6y4wr6ZUcRfQWzYsOqstykH2w9gJACEEHglLxy8E27BDa1KC6s+LSncNdGGKYEoQUhwwriUxb8fYE2so/KM1MEn83B6nagUKqFSqWDRWRLy0Fa8YEogSvz2G3DoEE0iOGxYvKWJIzxPNSLbDwhCVX4aSnpmvMVo1oiyhEqhEg7eAUVRYNIlRurmRIMpgSjh3xAeMoTuCcQLSQp+Hmqm3Ib612dF12tZOz3gOA0gn71v7XvVd+86fRp4/UxW/pkcAI31DjT0e2no11uzP+f1oNVVF4Lv3BMVE+6B0HtAw/I01B5DL0ao7/VM1wANyx6p9yorMlyCG07BAV7mYdQYoVfrQGRAOvvliUeEsls3BFMCUcK/HxCPqCBJom54QaD7sY2esBsBx51hUqjdLkkgThHg9ICXa7jfGWSq3dTU93KmySlUBdngZFZPe9oXb4MD4Bw9HqZtG6FyVsC08XuYNn4P77n9UXbjPfD0Hxq4OJJyRoIzKtVGvhbOewrlvRIQeAQvnLwTXskDnVoPvdYOWaqzBkkq1OroGtJMCUSBP/8E9u6lE/DFF8fmnrUnfpOJ5iqqJ615EKEqg4aUQL1tLh4Q3FVFNeQz9m1InsZM+pFQAuFMrI1WAieOQ/3SQ+A8HmS2twHXXAal/4/g/rUE3McfwfjLJrS4bzxIt+5QptwNcsmlIKr63RbxUABAfJTAmf5Wte/pFb2o8DkgiZXI4VQwa8xVm74Nx/8nCy69B+n26KkxpgSigN8KGDyYnhSOFv7Am5rV67KzaQElvT4BUlZzAqAjgCHegsSZxx6mYbL9+kF92aX0nH77QmDB08CM6cCrrwJvvw3u191QT/k/4IcfgA4d4itzksBLPJy8Ew6fAwoUpJvMUAcUaJw0ZoTRaQmzBJKNmrUDIo2i0InfX6rXaAQyMujK32BIgIm/Ji4XCw399lv6gVCrgfnz65pleXnAvHnAXXcBS5cCxcXBCmDNGuD886tLIzIAAJIiwck7UeGtgKjQ/DxadYp/1sKEKYEIc/w4PSTGccAll0RmTEWhq32fj45rNAL5+Qk68fsRRXY+wOsFHnyQ/jx5MtC1a8N9MzKAe+8Nbjt8GLjxRmpOTp4MTJgApKVFTdx4QAgBATnjv7X7yUSGw+eAR/LAoDbAZkjtw15NhSmBCOPPGHr++dQ1Ey41J36ATva5udUTf8JHXPrPB6RyPeFFi4CjR4GCAmDGjNCvLyqi1x49Cjz1FPDSS8DEiVQhZCZHeKlCFFT4KiDJEhQoUBQFBASKokCBUneSr/qP/l+tDGpCQKBVaWHT29hhrwjAlECEaaoryOsNnvizs2mIqdGYBBN/TXy++nfwUoXTp4HXXqM/P/ZYeHHCAwYAP/4IfPEF8OKLwP79wAsv0HHHjwemTaMWRAJTyVei1F0aKKjC+f/juEBiNo7joOJUgTb/xF77OSM6MCUQQU6eBDZtoj+HExrq9dIon+xsuoA2GpO4JK/bndpVxDIzgeXLgRUraDWhcNFqgbFjgSuvBL75Bli8GPjlF+Cjj4LdR599RhUux9HVQs1HVhbQp09133Xrqg/wqVTV13AcdTd16hS+vDUQZAFlnjIY1AYYtKmaNyXxYUoggnzzDf1u9ewJtGgR+vU8T+eOnJzIyxZT/PsBqawEAPpB6NkzMmOpVNS8vOQSah0UFdHQW4B+6KZMafjaoUOBt9+ufj5hQrW5WZuBA4H//Kf6+Zo1dKM6xA80IQTl3nLwCg+7wR7StYzYwpRABGnqATFCmokL3V9POJ5HpeNFeTl1BUUrxJPj6qakJQS44AL6r6JU/+t/1Jalc2f6N6rZx39dbm51P54HJk2iCqNNG+qeGjCAKopWrc4opkf0wOFzwKJlNSQSHaYEIkRFBbWygfD2AwSBWv56fUTFig88T/9NRV/u/PnAxx8DjzwC3HRTbO6pUgEfftj4/v7Mhmfj5EmqMH75hUYqHT5cfZ+WLYFbb6Wb1LWQq0orqlQqaNRsikl02F8oQnz7LfXnd+4MtG8f+vU8T62AZuFBcbmayRsJka1bgXffpT+fc058ZYkELVtShVFZSeukbtwIbNgA/PwzcOxYcGKqEyeoAhw4EJW9u8CdZYCNuYGSAqYEIkRTo4JEsZlUXxQEqtFSTQlIEjB7Nv356qup26S5YLXSfYWhQ+lzt5sqvJqrnQ0bgE8/BT79FHYAlpwsiOf3hXBeX/AX/QVKbrJvdDVfkinoMGFxu+n+GRDefoAsU4u+WRSe8Z8PSLVDYv/+N7BnDz3YNXduvKWJLmYzTY/bsmV1W/fuIFOngu/XB4pWC03pKRhXfAPbQ08i69KroNl3IH7yMs4IswQiwHff0b2zwkKgS5fQrxcEuhfQLPYD/MeaU2k/4MQJ4Jln6M//+EfSHOSKKJ06wTX1DpyYdAUsigbGXXuh27wN+pXfQuXx0JUOIyFhSiAC1IwKCmfu43kayp1Uh8Hqg5DU3A946CFqDvbtC1x3XbyliQuSIqHMWwaNSgO1wQyhfz8I/fvBfcM1ACEgGenxFjHpEBURj+x4DjdmXwwb2kXtPsk+7cQdnqebwkB4+wH+aL5mkR9MFOO7H7BrF/Dmm4DTGbt7EgL06kX95k891Qw0eXg4fA54JS9M2uAYZ5JuD1IAqpOnYi1a0vLUzy/jjf0fYsL2uRCV6KXETs1PbATZsIEuAvPz6VwQKqJI58xmsR/A83SDNB77AW++SQ9S/eMftIjD5s2xuS/H0YNaW7acOUFcM8YrelHuK4dJYzpjigfjfz5H9rDLof9+bQylS15u7XQdutnPwT/OuRVaVfS+U0wJNJGaVkA4i0Cep1ZAs9hHbegUajQ4dIgWcvZz8cVUk+bk0PDFsWOBp5+mWjZa1Kx6YrVG7z4JjEIUlHnLoBAFOs2ZLUDNnt/A8Txs982F+tCRGEmYvOSbcvDNyHcwNPu8qN6HKYEmIIrA99/Tn8M9JSyKzWT+IISaRNHc3ZYkmptj/Hhasefxx6tfKyig8etr19IQTUWheXbGjKGl3iLN77/TnEAbNkR+7CTCxbtQKVTCrD376fDKB2ZA6NMTqkoX7HfNBOf2xEDC5OLjg19h88mdgefqBirMRRKmBJrATz8BDgcNBjn//NCvl2WaIK5ZRAX5zwdEw6QpLaVpmQcOpGkMfvihOulZzZW+xUI16qJFwCuv0Nw6O3fSFMyRhBDggQeAX3+l90lRBFlAmbcMerW+cZOVTouK55+GnJ0F7YE/kDbn0fjVy0xAfirZgmmbHsY139+JfY4oLFwagEUHNYHPP6f/XnJJeNk+BYHuBzQbJRANs+af/6RplP2nU9PTaQTODTfQfDYNMXo0zZz51FM0eieSfPoptQAMBpomOkVx+Bzwyb6QEsQpOdmoWLwAGTfdBuPXqyD26ArPLTFKr5Hg9M7sjovyByJbn4lz0trG7L5MCYSJogBffkl/DveUsM9H00Y3i4ASrzcyb6SykhZM9odLtWxJFUDfvjT75WWXNX4XvUULmn/fDyHArFl0jCFDwpOvogJ49FH687RpQOvW4Y2T5HgEDyp8FY1yA9VG7NMTztn3wvbo07A+8wKEAedB6hbGAZtmhkljxNK/PBvzGgoJPf3Isoy5c+eibdu2MBqNaN++PR577LFANSKApqydN28e8vPzYTQaMXz4cBw4EP3TiRs2ACUldOE7eHDo1xPSjLKGRuJ8wK+/0gm6T5/gZGhjxtB9gC++oJu9TQmj+vRTmtvn+uupdRDORvZTTwGnTgEdOwK33x6+LEmMrMgo85aB47iw6/p6r78anr9fDteUyZC6RKZ+QTJy2HUMnx76OvBcq9IGCvDEioS2BJ5++mksWbIEb775Jrp164atW7fi5ptvhs1mwz333AMAWLBgARYvXow333wTbdu2xdy5czFy5Ejs2bMHhijGXZaW0rDQvn3Dm/ua1SnhcEtJCgItuvLWWzTE0s+6dbSMIkAtgu7dIyPnqFHAtm00nPT11+mmzgsvNP6Y9/btwDvv0J/nz0+9Q3FVOHknXKILafom1DvmODifnJdaJ8trUSE4ce33U3DIdQySImFcu9FxkSOhLYH169fjiiuuwGWXXYbCwkJcddVVGDFiBDZXxYATQrBo0SI8+OCDuOKKK3DuuefirbfewokTJ/DZZ59FVbYrr6QRiv/4R3jXC0IzCg31nw/QhLimmDYNuPtuqgA0GurH//jj6rKMkcZoBJ58kiqBrCxg717qGnrtNerfOxuffkqtnquvppvUKQgv8Sj3lsOgNkDFNXH6qKkAfD4Y//N5Sm0U27RWXNryIrSxtMCQ/PglHExoJTBo0CCsXr0a+/fvBwD8/PPPWLduHS6tcsIfPHgQxcXFGD58eOAam82G/v37Y8MZQvd4nofT6Qx6hINKRavxhUOzyRoKhLcfIIrAqlX057vvpoe7XnkFGDQo+qvD4cPpAY+hQ6kCe/hh6oo6G489RiONmnuCuAYghKDCVwFRESNbLlKSkHHDbbDNeRTG9/9z9v7NBI7jMK/3NHwz8l3kGrPjJkdCK4EHHngA1157LTp37gytVovevXtj2rRpGD9+PACguLgYAJBbsxpS1XP/a/Uxf/582Gy2wKPVWaokRRp/aGizOCWsKHQ/IFS/1smTtJZtVhZw//3BFa1iQXY2dUM98QTNijlu3Nmv4Ti6R5GKCeIAuAU3HLwDFl2EVy8aDXyX0IVc2pPPQLv958iOn0BIioSP/vwyaF/TpovvQaGEVgIfffQR3n33Xbz33nvYvn073nzzTTzzzDN48803mzTu7Nmz4XA4Ao+jR49GSOLGwfN0zmwWLmV/aGiofq2CArofsGVL/MKjOI7uPWzeDJxX41TmTz/Rg29+3n+fHghJYfwJ4lScKioHmDy33AjvJcPBiRLs99wPVenJiN8jEXhg61OYuulhzNo6P96iBEhoJXDfffcFrIEePXrgxhtvxPTp0zF/Pv0F5uXlAQBKSkqCrispKQm8Vh96vR5paWlBj1giCDSqqFmEhgpCePsBfhJBE9rt1T8fOEDLQo4YAezYQY+Ez5xJXUee1D3h6vQ54RE9YYWENgqOg/PJhyB2bA/1yVOwT50FCFFM+REnzss6FzqVFn/NS5yiQwk9DXk8HqhqzZRqtRpK1SZe27ZtkZeXh9WrVwdedzqd2LRpEwYm6MadPzS0WWQNBejEGOpJOVmmLqRExOmkB9IOHQKuuAKoikLD5Zc3k3je0PEniDNqjVGNXydmEypefAaKxQzd9p9hXfBc1O4VL65pdzk2jP4co1oNjbcoARJaCYwePRpPPPEEvvrqKxw6dAjLly/HwoULceWVVwKgGyvTpk3D448/ji+++AK7du3CTTfdhIKCAowZMya+wjdAszolrCjUbRLqav7nn2nGzRtvjI5cTaFvX7ppPHo0VVZlZUBeHnDvvfGWLC74N4NlRYZeE/0PrVzYGo5naE4o44pvoDp1Our3jDZrizfBI3kDzwtMMd7/OgsJfU7ghRdewNy5c3HnnXeitLQUBQUFuP322zFv3rxAn/vvvx9utxu33XYbKioqcMEFF2DlypVRPSPQFHieuoKaRWio/3xAqGFOGzbQCTZRNaHdDixZAgwbBnzwAXUHNZtQrtCo5Cvh5J2R3ww+A/xFF8Lx+Fzwg/tDyUruTfhvj/+IiT/eiz6Z3fH+X1+EWZt41mRCKwGr1YpFixZh0aJFDfbhOA6PPvooHvUf5U9wJIkGozQLeJ5aA6G6g/zhu4lcjJ3j6HmAq6+OtyRxQ5RFlHnLoFVpY5LNsibeq8fE9H7RwqqzwKIxobWlBYyaxFyYJrQSaG74908T1EgJHY8nvPMB/oIvCbpvw6A4eJogzqa3xVUO/eo1MH7+FSqemx9epsY40j+7N/478i20Muc3/XBdlEhMqZopzSpVhKJQJRDqfsCuXXQfwW5vfLoGRszxil5UeGmCuFgmM6sNV1YO28wHYfhmNSzPRTgleJQ46i5Ckac08LydtXXYlcHUv/8JdZT3RZgSiCH+/YBmkS6F56t3uUOhpiuoWcTINj/81cIISNgJ4iIFyUiH8/EHAQCW196E/pvVZ7kiNMr54PMfx9xFeOv3T/DV0e+C2lce+wFv//4pjrpOBNqKPSfxzu+f4ssjqwJtHsmL23+ahdGrbsYfzsNhy6X+4yDSb70b2ZddDfv7n4U9TmNg38IY4T8g2GxcQTxfffQ5FPxKgLmCEpZKvpJWC9MlxuaV77KRcE+6AQBge+BhqH9vesEVQgj+uesV9P18FH4p2xto/63id8za8iRe2LM0qP/Le9/C/VuewO6KfYG2PysP474tT+CZXa8G2kwaI6waCziOg0nThDhwvR66nzaCcBzUzsrwx2kEbE8gRvhdQc1GCXg84R0QGz2aRtpccEHkZUoxFKLAxbvgEho+c0EQekI2n+SLTIK4CFJ5793Q/Pob9Ju2Iv22qfDcPB7eK/8GEmbUFsdxOFh5BF7Zh+WHV+LcDOqazDZm4dKWF6HQ0jKof//s3sg0pCPXkBVoy9DbMbLFELQwBR9MHd7iAlxc8Bfkm3IaJ0tZOUwffAJ1UTGcj1GrR25ZAOdjD0I4vw8cRhWiuSvDEZJCafsawOl0wmazweFwhHR62OUCDh+mZ4vOfg/qCmrRogmCJgqyTA9TAc3o1FvyQAiBR/Sg3FsOl+iCRqUBh/B8jPVdp+JUkU0QFyG4snJkjb0B6hPFIBo1Std/C2Kj31euwkF/PoOv9Q/nYeQYM2HVUsVR7DmJrad+xmWthsVl30P9x0GY33wfxs9WgON5EI7DqZWfQi4MLlTkOnkcthbtkF3YLaTxGzuvMUsgRkhSnELNBYGu2CPpf/fXE45xug1G1YatrwKVQiVUnApp+rSEWrFHE5KRjtP/eRuGL/4LdUlpQAEAQPrtU6E6eRr8sCHghw2B0K93kKX6+r738ejORbjlnGvxUO/pAIA8Uzb+1np4nftE900Q6DZugWnZuzD8sC7QLHbrDPfNN0BukR9becCUQEzwh4bGPCrI6wVOnKCrI6uVrtoNhvDz/Pjxnw8IVbF8+y1QWAi0b99Mdsdjh0/yweFzoJKvhAIFZq055rH7iYCSmQHPzTcEtXGVldDs/x0qjxeat96H+a33odht4IdcAN/wv0IYPABtra0hKhIOVh6FQpS4KU7Dl1/Dfh9NRU44DvzQC+G++QaI/XrH7TsR8mxQUlKCmTNnYvXq1SgtLUVtb5IsyxETrrkQt9DQigpaQtFgoKXQABrNYzJRs8QvVKgfPrc7dEUiScCdd9Jr//c/oFtopm2qIsgCKvlKVPgqICkSzFozNGq2dqsJsVpRun4V9D9tgn71Ghi+WwNVhQPGz7+C8fOv4L1kOIY9/zRWXLwMfTK7x9T1w5VXQF1cEiihyQ+/CHJeLnzDhsBz03V1XD/xIORP08SJE3HkyBHMnTsX+fn5cY0hThZ4nqbLj+mvyuul6Y/NZjrxG400REkQaDH38vLqk2sWC31drz/75C7LdGx2PiCqSIpEJ39vBXiFh0ljSphonYTEaAQ//K/gh/8Vnx1dg9fenYnr/jDi5oM28BddCADom9UD6oOHYZv9CPhhQ+Ab/lfIbdtERRz1n4dgfvM9GD9bAal1S5z+4kOA40BMRpz89vOEyhsTshJYt24dfvzxR/Tq1SsK4jQ//IZSzPdPKyqoy6bmZM1xwSaJJFGlUFJCX9PpqNIwm6sLHtTWXP7zAaHuB/hDQ/v3Z+cDzoA/4qfMWwafTKN07AZ7vMVKKgbkn4/7umTj00GdMei8OcjS2QOv6VevgW7Hz9Dt+BnWZxZDalcIflD/wPfEPekGKNk0Akj343rof9pIL/Q7PGp4PtwTr4dSQH34ug2boV+9BgCgPnIUhjU/VQuk0UB1uqw6D1ICKQAgDCXQqlWrOi4gRsP4C8jE1BXk8VRbAWdCo6EPk4l+uHmeKo/Tp+kHtbaVoFaHvx/AzgecEUII3IIbFb4KuEU3tCotbHobs7Qbwa/l+7Hy2A+4t8dtAGis/jcj30WWoW7Ynm/0pSBGAwzfrYFu4xZo/jwEzZ+HAq97/z66Wgns3AXz0ncbvK/3b5cElIDm199gfvuD4HtddCE8N4+HcH7fhN4DC1kJLFq0CA888AD+9a9/obCwMAoiNS/8i+am7sWGhMNBJ+pQVhwcRyd9/0EGUaQTvstVbUFYLPQNhbqSkSSWL+gMeAQPKvgKuAQXVJwKVr01ZSJ+mspJ72mM+t9NEBQR/bLODRRsr08BAICSmw3v+HHwjh8HzlkJ/dqfoNmzD1zVwlaxV0fkC316wnXrhOqL/fN41YSuZFdnOBV79YDr/yYBAIheD9+lF0fN1RRpQp6arrnmGng8HrRv3x4mkwnaWhNCWVlZxIRrDsQ8a6jHQ1fzTb2pVls92SsKVQinT9OfQ4113b2bKhO7ndYRYACojvhx8k4QkJSN+GkK2cZM3NThKpR4T6JDWmFI15I0K3x/uwT42yX1vi4MHgBhcOMy3Yr9etMInyQkLEuA0ThinjWUkOq9gEj6HVUq6hIKd2OD7QcEIcgCnLwzUKyFRfzURZBFHHQdBS/zgdO8ALCmaCNe/u0tvDDgUeQYqdvm4d7TmfJsAiF98kRRxJo1azB37ly0bds2WjI1G3ieKoCYldH1RwRZrTG6YSMZPx7o2DFlC7P4kRWZTv41In50iVBjOUYcdxfj4R0LoVVp8fKgJwLts7Y8ic8Of4NZ596JSedcAwAo9pbir/+9Gka1AX+Oq95k/eLoKqwt3oQFu17BM+fTFAtMATSNkJSAVqvFJ598grlz50ZLnmaFINCUEjHZEyKEhn0SEuMNiEaQlgYMj/HJzATDI3hQ5i2DS3SlTMTPr+X7YdYYUWhtBQAQFBErjq6GWRNcXUtUJDhFFypFd6DNpDEhXWeDWWsKOtzVJ7M7TGojpnSZAEZkCDl30IQJE9CrVy9Mnz49WjLFnGjkDiKELsrbtInRnoDbTYWxWBJPCaQwoizCwTtQ4a2gfn+dOSU2fb2SDxevvB5F3lIs/cuzuDCvP1yiGx8dXAGLxoRx7UYH+p7wlMAr+ZBtyECaLsGs2AQg4XIHdezYEY8++ih++ukn9O3bF+ZaM9w999wT6pDNEp6nbqCY7Af49wKAxFMAn38O7NkD/O1vQI8e8ZYmZhBC4BJcOO05DZ/sg1lrjntu/ljiEt3INWbBJbnRI70zAMCiNQfcPTVJtMLrqUbIM8Ybb7wBu92Obdu2Ydu2bUGvcRzHlEAVPE+DYWJSDc9/LiARfe7/+Q/w3XdAZmbKKAFe4lHuLYdTcELNqYPj/X0+cD4exGQCdM1XKWQbM/Hx0FdwzF2E9DiXp2ScmZCVwMGDB6MhR7NDlukZrKjj3wvguMSzAmqeDxg0KL6yxACFKHD4HCj3lkNURFh0luBNS0Kg8vFQ0qzg3B7A42nWykDFqdDa0hxypzdvEmzWaB6IYvWB26jjdtNiBYloBfz6K904sdmafb4gj+BBua8clUIlDGoDbIa6q1/O7YFiMkHJywUEASqHE5yjslkpg7nbnkGBKQe3dRrPonaShJCVwKRJk874+r///e+whWku+LOGRj36L5GtAKD6fMD558fILxZ7JEVCha8isPHbYH5/RQEnSVByc+jvwmiEYjQCtrRmowx2nv4Vr+9/HwCtxNUnKzXcf8lOyDNHeXl50HNRFLF7925UVFRg6NChERMsmeF5ICMjBqGhfisgUYu7rF9P/22GqSJC3fjl3G4oFguIpVaoWH3KwOsFMRkTLtHY2eiZ0RXPnj8Xh1zHmAJIIkJWAsuXL6/TpigK7rjjDrRv3z4iQiUzilKdhieq+K0AtToxV9nNeD/gjBu/9SFJgKyApNsbPjFdr2WQXMqA4zhc335MvMVghEhEApZVKhVmzJiB5557LhLDJTWCEJyHLWq43bQuQEwTE4XA8ePURZWW1mzyBSlEQYW3Asedx+HgHTBrzTDrzGfN9Mm5PSC2NBBzIyIFjEYoebmQW7cAsdvAeXlwDifdaEpQ9jsOQlQSVz7GmYmYI/mPP/6AJEmRGi5piUloqN8KUKkS0woA6Cm5X34BiooSV8YQ8IpelHnLzrjxWy+iCIADsdtD8w8miWVwyleOsatvQ4E5F0v/8iyL+U9CQlYCM2bMCHpOCEFRURG++uorTJjAjnLLcgwW5y4XtQISLUdQbVQqoEVyhwjWPvEbamF3zu0BSbfTyTscElwZ/O48BJnI4GUBmfr60zczEpuQlcCOHTuCnqtUKmRnZ+PZZ589a+RQc8cfGhrVAjKKkvhWgD8TSQIX0jgb/vKODp8DPtlHk71pQgz3EgRArQ7KUR82VcqAS0sD50wcZTAgpze+H/UhHIILenXqJMNrToSsBL7//vtoyJGciCLg8gF2C8Bxgf2AqIaG+vcCbAl8CnPXLuDmm4ERI4D58+MtTUhIigQX70KFrwI+yQeDxhB2hS/O7QHJzoroBhExGUFMNZRBhRNQ3HUVbkPy1tdeu4njQHQ6qlwa8b5zjdnINWY37g0wQkKURShEgSaKR7pCHnno0KH49NNPYbfbg9qdTifGjBmD7777LlKyJT48D5w6Bei9QGYmeF4d3dBQRQHKyuiGayLn5V+/HiguBk6ciLckjSZo8pd90Kv1sBmaUN7R5wN0Oihp0XHZ1VQG8PlqvVg3JyRHSL3tUJS6bZIEzusD5/WBACCG4HrThBDM2/4srmp7GXpmNO9DgPFEkAR4JS9yjFmwG+1Ru0/ISuCHH36AIAh12n0+H3788ceICJVUiCJQVgaFF8Hps2E0RtEMcLurT+AmMklUT7jeyb+ptX2r0kPIOdlRLy5NTEagEfsNIVcFF0VwPh6cxwPO5QYcTjqOXoePT3yL1/e/j/f+/AxbL/+K5QaKArzEg5d5ZJuzYfdxUa013Wgl8MsvvwR+3rNnD4qLiwPPZVnGypUr0SLJNwHDguOAtDTwJythMBHolUwAYW4CnolksQJkGdi0if6cwOcDJEWCW3Cj3Fseucm/Cs7rA9HrQaJkBcQErRZEqwWxWoBMiSoErxdwuTA8rRcuz/sremZ0QbomAdOVJDk+0QdBFqgCMNjB8Y6o3q/RSqBXr17gOKqR6jsZbDQa8cILL0RUuKRBpYJgtCEdDqiLfUBubuQjd/wRQbXccAnHr79WRy51Cy3/eSyQFRkuwRWVyR8AQAg4QYCcl5sQ0TsRQaMBsWjoaefMDKT7cvFKixeAShc4lwdQZKo09Lrm857jhE/0QVRE5FpyGx+G3EQarQQOHjwIQgjatWuHzZs3Izu7eiNIp9MhJycH6kSNVokBMlHBlJ8GKJX0oFRuLp2wIzG5+K0ArTaxrQAguJ5wAn0eoj75V8F5vFBMxuS2AhrgtOc0Mk2Z1S6ojHTIPE/3D1xucF4v3QxXq0BikjyreeEVvZAVGbnmXKQZYpcKptFKoE2bNgBoighGMKIIaNUEBj0BdBZa6/fECZouIDOz6RO3y0UfiW4FAAm3H+Cf/B0+BzySBzqVLiqTP4DqJHHZWQmlACPBn+V/4pJ3LsGN596IBy54gOZJqsqPQgwGELsNEATqNnK7wXm81C3GcdRCiPLeSLLjETxQoCDXkgurPrYLiLBmp7fffhuDBw9GQUEBDh8+DAB47rnn8Pnnn0dUuGSBFzgYDKTaEjYa6YmxkhIaJdOUk9SyTK0AnS7xrQAA6N6dPuK8HyArMpw+J447j6PIVQRJkWDT22DSmaK2ycZ5PFCMRupHb2Z88/s3cItu7D65u/4U0RwH6PUgtjQoBfmQ27SC3KIAxJYGThDBebyxFzpJcAkuEBDkmfNirgCAMJTAkiVLMGPGDIwaNQoVFRWQZRkAkJ6ejkWLFkVavqRAlIA0Cwn2/Oh0NIqnrIxaBTwf3uAuF40KikmFmggwcybwzTfAuefG5fY1J/8TrhMxmfzpjWVAkkEy0pNDWYfIHefdgaVXLMWzI55t3Inpqk1lJS+Xps8WhKYthpopLsEFFVTIs+TBoo/P4iHkT+sLL7yA1157DXPmzAnaA+jXrx927doVUeGSAUICi6C6qNXUhVNZCRw7RstAhoLfCkiGvYAEgBd9KCn5E0Vlh2I3+VfBuT0g1npSRTcjRrQfgZZpLUO+jqRZQdLt4CpdUZAqeankK6GGGnmWPJh18fvchDyzHDx4EL17967Trtfr4Xa7IyJUcsFBrwP0ugYisVUqID2dbhwcO0bz/zcWvxWQqJlCa7NnD90PiQNuxymUHNsHl+JDmqKDSY5ubHUQkgQQApKentSpMmojyiJe2PwC3EITv9ccByUjHTAYaFnNFIcQAifvhFalRZ41DyZdfK38kJVA27ZtsXPnzjrtK1euRJdmXkKwISxmcvZ9wLQ0qhCOHQNOn67/9GZN/FaAXp8cE4ssA2PH0rTRv/8es9sSnw+OE4dQ5CwCn5MBS4cuQEEBOFGO2YQTSBUdbpK4BOWFzS/gqXVPYdzH40DO9nk9GzodlKxMqjATOC12tPErAJ1Kh1xLLoza+H9mwsoiOmXKFPh8PhBCsHnzZrz//vuYP38+Xn/99WjImNBoNQRGQyO/ICYTVQRFRfTLkHWGKBK/FZCeJJkZ9+yhVo7VChQWRv9+kgTZ6UC54MQpowxtZmuYzDSumtjSoKhUUJWUgqusBIlmtlVBBDgVTRKXDMo6BM5rcR4KrAWY3HdyRKwqYrVQt9DpMhpN1Mx+X2fDrwCMGiNyLDkwaGJRhPzshKwEbr31VhiNRjz44IPweDy4/vrrUVBQgOeffx7XXnttNGRMaHS6qtDQxmIw0In/5EmqCHJy6h6wSTYrAKguJXn++dGtd6wogMsFQeJxSiui3KKHxZ5dp7QjsVqqFYHDSeP2o/C75NxukKxMGhHWzPhL679gzcQ1MGoi9N6q3EIqj4daT814/6Q2hBBaiEhjRo4lB3pN4oTMhvRtlSQJ7733HkaOHInx48fD4/HA5XIhJycnWvIlNCoVYNQpoZ+J0WrphnFFBTWN8/KCM01WViaXFQAAGzfSf6MVGkoI3VgXBHgMapSa1XDpVLAb0xuMViFmE+T8XKhKT1JFYEuLrCLgeUCrhWJL0BrPYSLIAnRVaaFN2gj7q7VakKwscMdPVOdeb+YoREElXwmL1oIcS07gd5sohLQnoNFo8H//93/wVWUtNJlMKasAAOrdyc2Sw7tYpaKKwOOhJ4xdVZETkkT3DJLJCqiZL2jAgMiP7/UC5eUgKhUcWVYcswI+owbpxoyzhytWlWskZjMt0xjBw46cx0vrBjejg1BbTmzBX5b+BT8c+iFq9yAWM0hGOlQu99n3xpIchShw8k5YdBbkWnITTgEAYWwMn3/++XUKy6QyTZqnOY4qAlmmiqCigloBXm/ynAsAgL17AYcDsFjoQbFIIQjULaYoUPJycSrLhOMqF9QaLdL0aY33U+v1UPJzQawWqBxO+vtuKj4fHTeO6SG2HN+C2768Lajt0TWP4sKlF+LjPR8H2o5XHscF/74AF799cVDf+T/Ox6A3BmHZzmWBtsUbF+OY8xg+++2z6Ale5RZSzCZwzTiiUFZkOHkn0vRpyLXk1nFZJgohO2/vvPNO3HvvvTh27Bj69u0Lc63wxXPjdEgoqbHUSDWh11PXULJYAUDk9wNkmSpDlQrIzoZoNaNUcqDcVw6LzhLeakqrhZKXC6jVUFU4oFjM4ctKCFReH+TcnLjlx/m5+GeM+XAMMo2ZQe0lrhL8Uf4HHL7qzJOyIuNgxcE6rp3T3tM47DgMJ18dtvyv0f/C4k2LcUe/O6L7BjQakMxMcMdO0M11XWJOkOEiKzIq+UrYDDZkm7OhUUVxn6yJcCTE2C9VPYeWOI4DIQQcxwVOECcTTqcTNpsNDocDaWkh+HddLuDw4cj57gWBKoO0CPuuo82ffwKrVwOtWgGXXBL+OIpC90JkmZ62Tk+HV8uhxFUCl+iCTW+rP2VBKMgyVKfLoDpdBsVkCmvy4TxeEBUHpVXL6G6CnwFCCG5afhPyrfl4YugTgVXmgdMHcNp7Gq1trVFgLQAA+CQfdpXsAsdx6FfQLzDGoYpDOO05jQJrAfKt+XF5H6pTp6EqPRnX6CpZkaGQxrkJuTpl2OqiEAUeyYN0QzqyTFlN/8xWVNBIwhBd742d10JWAv5cQQ3hTzSXTCSMEkhVCKHKz+ejVlFmJmCxwClUosRVAlERI5v0TVGgKisHd/IUje0PZTVPCFQOJ+T8PBrmGCMOVxzGkq1L8OhFjwYsIVEWE9bF0GgkCarjReB4X3RDeeuBEAK36IasyBFfqVt0lsgoACDqSiDkd56MkzwjgeF5uvo3GoGWLQGrFYqKQ5nnNErdpdCqtbAb7JG9p0oFJTMDHMdBVXoSRFEaXQeYc3ugmEwxTRInKRKu+c81OOo8inxrPqb2nwoAya8AAOoWyvK7hYSYudcESYBH9MCkNSHDnBHxkE2NShO7E+tNJHEdVYzkYM0amin1wguB/BBdCoRQBeCvvaDVQlIknHSV4rT3NMxac/TiqTkOJDMDiloNVclJQPGe/cSvP1V0bk5MU0VrVBrMGjwL7+9+H387528xu2+sIGYTlKwMqEtPQoly5TyFKDRpG6dCjjkHaYa0hPbXxwKWlYzRNN58E5gxA/jss9CvFUW68qtSAD7Jh+PO4zjtPY00fVpMDtQQuw1Kfi6d4F1nTnDGud1QLNFPEufwOfDAtw9gy/EtgbYxncfgw6s+RPv09lG9d7wgdhsUs5nWM44SPtEHJ++EWWtGC2sLZJgyUl4BAEwJMJqCojTtfIAo0mgorRaVfCWOO4/DLbiRbkiP6ZeTpFmh5OcBKjU4Z2X9nSQJkBV6LiDKGV0XblyIt395G7O+nQVZoYEW/tKuzRa1mp685rjw0643gKzIqPBVgBCCfEs+8q35CZGzJ1FgapBxVhSiQJRFSIoESaE54TlwUO/9DeaKChCzGd5O7cGJNIMox3HgwJ3xXwCAIICkp6PcU4ZSdylUKhXsRntc3iOxmCGr86AqLqk3zUQgSZw5+uc3pvWfhl0luzBz0MzIbCwmCcRkhJKZAXVJKZQIpE8nhMAjeiATGemGdNgN9oRK15AohKUEKioq8J///Ad//PEH7rvvPmRkZGD79u3Izc1FixYtIirg8ePHMWvWLHz99dfweDzo0KEDli5din79aKgbIQQPPfQQXnvtNVRUVGDw4MFYsmQJOnbsGFE5UgVJkQITviiL8EreagVApKAQufQ1q2AG4O7dDcc9RQAQNNnT/2v8VzWpchwHFVRQuyoh21VwKDJMWlP8E2oZjVDy86AqOQmuwkHTTKhUVVkvOZBI1YyugSiLeH376zjpOYl5Q+YBANKN6fj0mk8jep9kgaTboXi94Cpd9PcfJoIkwCN5YNKYkGvMhVlnbt6WVBMIWQn88ssvGD58OGw2Gw4dOoTJkycjIyMDn376KY4cOYK33norYsKVl5dj8ODBuOiii/D1118jOzsbBw4cQHqNkMwFCxZg8eLFePPNN9G2bVvMnTsXI0eOxJ49e2BoZMRHKkIIgaiIkGQJoiJCkAX4RB9tUyQoUMCBg0algZpTw6g11lmVpm3/lY41YABsBltgXAISSD1c83ngX0IgiV5IUCCqZKTpE2hzzmCAUjPfUJqVWgHp9qikit5ZshOP//g4OHAY23UsumV3i/g9koqqyC2110dDhkP8DitEgVtwgwOHbFM2bAZb4ny2EpSwUklPnDgRCxYsgLVGXO+oUaNw/fXXR1S4p59+Gq1atcLSpUsDbW3btg38TAjBokWL8OCDD+KKK64AALz11lvIzc3FZ5991mBWU57nwdfwOzpDKfSShMiKHDTh+yQfeImHRKSAz1nFqaBRaaBRaWDQGs6ek0dRoNtK04cI/fsGmmtaAWeCExQQsxF6cwIWYtHpoOTlQqVSQVVeAaLV0sNMEcJ/sBIAzis4D7f1vQ2dszqja1bXiN0jqTHWcgs1MhLLJ/rgk32w6qzIMGYwv38jCdnptmXLFtx+++112lu0aIHi4uKICOXniy++QL9+/XD11VcjJycHvXv3xmuvvRZ4/eDBgyguLsbw4cMDbTabDf3798eGDRsaHHf+/Pmw2WyBR6tWrSIqd6IgKRKKKotwpOIIjjqO4ljlMZS6S+EWaQSGQW1Amj4NNoMNVr0VRq0ROo2uUTVk1QcPQ+VwQjGZIHbtHLpwokijbBJNAfjRaKDk5kDJzKCbwRGwKgkhWLF/BS5777KgVA0PDXkI13S7hrkrakDsNihp1kZFC8mKDIfPAYUoyDPnsY3fEAlZCej1+npXzvv370d2dnZEhPLz559/Bvz733zzDe644w7cc889ePPNNwEgoHRyc3ODrsvNzT2jQpo9ezYcDkfgcfTo0YjKnSi4eBccvAMqlQomrQl2gx02gw0WnQUGrQEadfgHWuT2bVGy4VuUv7Y4rNQJHMfFLe9Oo1GrqSLIyjx730YgKiKe/ulp/FzyM17d9mpExmy2VLmFoNE0WLKUEAK34IZbdMNusKNFWgvYjfZGLWIY1YT87b388svx6KOP4qOPPgJAv8xHjhzBrFmzMHbs2IgKpygK+vXrhyeffBIA0Lt3b+zevRuvvPIKJkyYEPa4er0e+maU/rc+JEVCha8CerU+ahERJCMdYkYYKTNEEUSlAtEnuBIAmmypCLIArUoLjuOgU+vw5LAnsfHoRtx1/l0RErAZYzBAycqEuqgYik4X5BYSZRFu0U2rdBlzYNFZmCUVJiGrzGeffTZQSMbr9WLIkCHo0KEDrFYrnnjiiYgKl5+fj65dg/2kXbp0wZEjRwAAeXl5AICSkpKgPiUlJYHXUhW34IZP9sU/4qYeOFECqTof0JxZc2gNLlp2Eb7+/etA219a/wX3Db4vIf8u8YQQAoUodP9KpoEKvMTDY9TAY9KBLz8Ft+CGS3DB4XOAl3hkmbLQIq0FrHorUwBNIGRLwGazYdWqVVi3bh1++eUXuFwu9OnTJ8gvHykGDx6Mffv2BbXt378/kL+obdu2yMvLw+rVq9GrVy8AdJN306ZNuOOOKKfCTWBkRUa5txx6tT4qXw7N/t9hfWwB+CEXwHPrTaEPIAhApKt8JSCbjm/CIcchLNm6BJd2uDQlJ6qaSdpqRo0BAAEJCjlWcapAKHHgZ5UKSlYGdD4RWoFAZbFAzalh1pph0iVRzY0EJuzYqQsuuAAXXHBBJGWpw/Tp0zFo0CA8+eSTGDduHDZv3oxXX30Vr75K/akcx2HatGl4/PHH0bFjx0CIaEFBAcaMGRNV2RIZl+CCT/bBpo9Olkvdxi3Qb94G6PVhKQEOgNIM3XGiLMIluJBupC6yu8+/G1qVNmKF2pMRJ++EXqOH3WCHCioacVXfZH+Wn2EsAI4dAwzWuKXvjguERLQaXn2E/NtcvHhxve0cx8FgMKBDhw648MILoY5Agq3zzjsPy5cvx+zZs/Hoo4+ibdu2WLRoEcaPHx/oc//998PtduO2225DRUUFLrjgAqxcuTJlzwj4j8jrVLqoTTy6zdsAAML5fc/Ssx4kCUStBmlmRUR+KfkF01ZOQ1t7W7xxxRsAAKPWiOkDp8dZsvjhFb1Qc2rkWfKaXqs4LQ3IyKClVzMyIiNgIuOvqc3zgNlMU6xHiZDrCbRt2xYnT56Ex+MJHNoqLy+HyWSCxWJBaWkp2rVrh++//z5pQi+bUz0Bp8+JE64Tkc2/XxNFQc7A4VBVOHD6o2UQe/YI6XLO4wVRq6C0ad2s3EH7Tu3DiHdGIE2fhlU3rkKeJbX3pPxWUYu0FpFLBS4IwNGjdGVsjm4Sv7hRs7aG0UjrCFgsYWWtbey8FvLG8JNPPonzzjsPBw4cwOnTp3H69Gns378f/fv3x/PPP48jR44gLy8P06en7gooXsTCCtDs/4OWZ2zC+QCYE/h8QCORFRl7Tu4JPO+U1QlLLluCtRPXprwCIISgUqhElikrsi5JnQ7IzqafIY8n6m6SmOP1AuXl9LvRsiVQWEgr7EU5bXnI7qAHH3wQn3zyCdq3r05p26FDBzzzzDMYO3Ys/vzzTyxYsCDi4aKMs+MW3PBInqjtBQCAbvNWAIDYt1dY0T0cSf79gFOeU7hx+Y34o+wPrJm4JlCacVTHUXGWLDFw8A5YdVZkmjIjvxixWmn9iYoKwOGgE6bBEJHDfHHD66UPoxEoKKCurxjue4R8p6KiIkiSVKddkqTAAa2CggJUVjaQkpcRFWpbAaqSk1Bs1sCXw/jxchg//RJC/34QBpwHofe5NI1ziDR5P0ClSvr9gAxjBtScGmqVGvtP749bfd5ExCN6oFVpkWPOiU7OHo6j5UftdjpxVlZSt2xZGV2UGAzJE3rs81GLxmCgBZnS0uIie8h/pYsuugi33347Xn/9dfTu3RsAsGPHDtxxxx0YOnQoAGDXrl1BOX4Y0ae2FZA29zFo9+yD4+lHIAweAP2PG6Db/jN0238GlrwBotNB6NWDKoT+/SD26tGo1Qcxm6BYzEH5ghqNKILoNIl/UrgWhBB8/fvXGNl+JNQqNVScCosvXQyrzopsc2RPyScz/vj+FtYW0U/boFZTX7nFQt1DXi+1DDweqhj0ejq5xrACXKPxl1TV64G8PDr5x/E7EfLGcHFxMW688UasXr0a2iqtJUkShg0bhrfffhu5ubn4/vvvIYoiRowYERWhI02ybwzLiowTzhMQFAFmnRnan3chc9xEELUap77+BHKbVlAfPQ7dxi30sWkr1CdPBa4najVKN38fqJilKjkJJTO9YaUgy3RFFmK+d85ZSVMF5yTXxDnxs4lY9ecqPDnsSUzoGf5J9eaMQhRUeCuQY8lBtik7PiGxhNAJtqZCIKTaXRTvfShBoHOGv5qezRaWNd5YolZoPi8vD6tWrcJvv/2G/fv3AwA6deqETp06BfpcdNFFYYjMCBeP6IFH8sCqp1ldLS/SJHveK0ZBbkMjtORWLeBt1QLeq8cAhEB98HBAIXC8EFQy0X7XvdD8eQjCeX2o+6h/P0idz6me9MNcXXGKkpT7AUPaDMG6I+sQ4noppXD6nLAZbMg0RmEfoLHU3B+w2ai7xeWilkF5Of3cGgxRnXjrRRDoyl+joRvbdnvsZTgDIVsCzZFktgQUouC443i1FbBzFzKvqbICVn4CuXWIYbqCgJwLL4WqvCL4PnYbhJ7d4fjn4+EV+5BlcC4P5MJWCb2JV+GrwJItSzCi/Qj0LaAuL1mRUeIuQYG1IM7SJSb+/P0tbS0TMx2GJAXvHwgCXY0bDNHdgBVFOvmrVHTit9tj+tmPmiUAAMeOHcMXX3yBI0eOQBCEoNcWLlwYzpCMMHELbrglN9L09I9seZGepvZecVnoCgAAdDqU/vQ/aH7bH7AUdFt2QFXhgGHNT9CMvQHly16B3DLECVEUQfTahN8PWPDTArz585vYVrQNH1/9MTiOg1qlZgqgAQRZgCiLiasAADrRW630IQjUTeR00n9lma7KI7l3QAi1QjiOLhDtdhr5k6CErARWr16Nyy+/HO3atcNvv/2G7t2749ChQyCEoE+fPtGQkdEAfj+sRqWBilNBu+MX6H9cD6JWw33nLeEPrFZD6tYFUrcu8NxyEyCK0P76G3SbtkCxWCC3CD0ahhPE6nKNCYR/M9Oso+6wKedNwc7inbit721xlizxUYgCl+BCrjk3sAhJeHQ6+vC7izweun8gy5G9j81GFYAp8fMbhawEZs+ejZkzZ+KRRx6B1WrFJ598gpycHIwfPx6XXHJJNGRkNEBtK0C3YTMAwDvmb5BbtYzcjbRaiL160AiiMOEUBUqCrYZ+OPQD5qyeg+Hth+ORvz4CAGiR1gL/Hf/fOEuWHDh8DtgNdmQYkzCNA8fR1bnRSCfrSHvFVar4b0Q3kpCXZXv37sVNN9GkYRqNBl6vFxaLBY8++iiefvrpiAvIqB+FKHD4HAErAADcd96KU5++A9ddk+MsXS0UBYRLvPMBHDgcchzCfw/8F7zEn/0CRgCX4IJBY0C2KbtO7emkQ6Wi7qBIPpJEAQBhKAGz2RzYB8jPz8cff/wReO3UqVMNXcaIMB7BA7forpOYS+rWBUpBgh1eEgRAF//9gM3HN+P7g98Hnl/Y5kI8N/I5/DDhh6gV3mmO8BIPWZGRa8llv7dmQMjuoAEDBmDdunXo0qULRo0ahXvvvRe7du3Cp59+igEDBkRDRkYtFKKgwlcROLik/vMQiNEAJT8xc9Zwoghitcb14M6K/Stw+4rb0TKtJdZOXAu9htZaGNdtXNxkSkZkRYZbdCPfkg+LLnqZLRmxI2RLYOHChejfvz8A4JFHHsGwYcPw4YcforCwEG+88UbEBWTUpbYVkPb4AmRffAUMX359livjhKyAxGE/QCHVCcaGtR2GlmktMaTNEPgkX8xlaQ4QQuDgHUg3pAdqJjCSn5AsAVmWcezYMZx77rkAqGvolVdeiYpgjPpRiAIH7whYAdptO6H/aROIRg2xd894i1cXRaH+0Ri6gk55TmHhhoUochVh6RVLAdDc/j9M+CH66QyaMZVCJcxaM7LN2ayYezMipL+kWq3GiBEjUF5eHi15GGfBbwUYNXQyC5wL+PvlocfuxwJRBLTamBaVr+Qr8e6ud/G/P/6H3aW7A+1MAYSPT/IBBMgx50CnTuyzHozQCHlPoHv37vjzzz9Zgrg44DfHVZwKapUa2q07oF9PrQDX/zXhXEAU4UQRxGyO+n6AW3AHYv3bprfFgxc+iK5ZXdE9p3tU75sKSIoEj+hBgaUg8DtmNB9Ctukef/xxzJw5EytWrEBRURGcTmfQgxE9PGItK+CFfwEAvGOvgBLGAa6YIEpR3Q9QiIJXtr6Cwf8ejOPO44H2yX0mY3DrwVG7b6pACIGTdyLTmMn2AZopIVsCo0bRwhmXX355UKIoQgg4joMc6ZN3DAD091vhq6i2ArZsh37jFhCtBq7bJ8VbvPohhO4HRDFZlqRI+HLflzjpOYmPfv0opWv6RgMn74RFa0GWKSt+ieEYUSVkJfD999+fvVMqcepUTBLI+a0As5aa45rDR6GYjPCNHpXAVoAI6HRRPSSmU+vw8mUvY92Rdbi+x/VRu08q4i8Un2PJgVadWAf9GJGDZRFFE7KIvvwycO+9wNKlwIUXRk0+QghOVJ6AR6xOFw0AXFk5TceQlRm1ezcFzuUCMRqhtGwRsTElRcLCDQvRMq0lm/SjiKRIqOQrI1sonhFTolZoHgB+/PFH3HDDDRg0aBCOH6d+2Lfffhvr1q0LT9pkhBDgu+9oEqpbbgF+/jlqt/KIHriFuqeDSUZ6wioAAIAk003hCPL5b5/j+U3PY+53c1FUWRTRsRkU/z5AxAvFMxKSkJXAJ598gpEjR8JoNGL79u3geZpzxeFw4Mknn4y4gAkLxwH/+hcwYADNRHjDDUCNFBqRwh8R5E9prNm9F9ot2yN+n4jjNzAjfD7gyi5XYvQ5o7Fw5EJW27cKhShBD1mRm/SIaqF4RsIRsjuod+/emD59Om666SZYrVb8/PPPaNeuHXbs2IFLL700UGw+mWhSUZk9e4BbbwV27QJatgQ+/5zWDY0QHsGDY85jMOvMUKvUyLjxNug2b4Nz9r3wTExgd4gggBMlyIWtm1S4wyf58N6u9zCh54TkT1QWYRSiwOmjEXn1TdYcwpvAtWotCqwF7FxFkhO1ojL79u3DhfX4v202GyoqKkIdLvkxm4F33gGuuAI4dIhaBJ98QvOJN5HaVoBu01boNm8D0WrhGzm06bJHEU4QQZpYuYkQgus/uR6bjm+Ck3di2oBpkRMwyeElHm7RDZvehkxT5hlP8IaqDFScim0EpxAhu4Py8vLw+++/12lft24d2rVrFxGhko6sLOD994GcHGDvXuDjjyMyrFf0opKvpHsBhATOBXiuHpOwyeICSBKIuWkFNTiOw/U9rkeGMQPn5p4bIcGSG7+/npd45JnzUGAtgElrgkFjaPCh1+hDejAFkFqEvEybPHkypk6din//+9/gOA4nTpzAhg0bMHPmTMydOzcaMiYHrVsD774LrF5NN4qbiN8KAAdqBWzcAt2W7SBaLdy33xwBgaOI38MYxvkAt+BGha8CLdJoRNFVXa/C8HbDWYQKaBU0f/6eHHMOO73LiAghK4EHHngAiqJg2LBh8Hg8uPDCC6HX6zFz5kzcfffd0ZAxeejalT78SFLYBSZ8kg8uwUXPBdS0AsZdCSUvN1ISRwdRAjSakM8H7D+9H5M+nwST1oQvrvsiULOWKQCqHEVFRLYpG5mmTGhUUSyQzkgpQnYHcRyHOXPmoKysDLt378bGjRtx8uRJPPbYY9GQL3nxeICbbwaeeirkS/2ngwlItRWwdUdyWAEAOEEA0esBbWhKIE2fBifvRJm3DMecx6IkXXIhKzLKveXgwKFlWkvkmHOYAmBElJA/Te+88w7+/ve/w2QyoWvNVS8jmDVr6DmC774DsrNpBFEjCbICAECUILVqAX7IBVByc6Ik8BnkkXnsKd+PdL0dba2tAu2fHvoaRo0Bf80bCGPVqt0teiB5y2Gyt0ZjYnlEWQz4oPMseXhzzJsotBeyPDWgnwOP6EG6IR1ZpixWxYsRFUK2BKZPn46cnBxcf/31+O9//8tyBTXEpZcCs2bRnx96CFi+vFGXEULg8DmgECUQEilcOAinvv4Erul3RkvaAKIiYk/5AdSMHH5o+7O4bNVEvPtH9XuQFRlTNjyIST/OhFvyBtpf2/ceOn97BR7YEnxm5G/v/Q1XfHAFSlwlgba1h9di0L8HYXtR9bmH3vm9U14B+D8Doiwi35KPfGs+UwCMqBGyEigqKsIHH3xAS/ONG4f8/HxMmTIF69evj4Z8yc3dd1dvEk+bBjQi75JP8gU2/4LQakEskS3npxAFDqEy8FyQRXT5ZCiGrbwWxd6TgfaeGV2RobcHhSEKiogL8/qjX9a5sNQ4yewVPAAAU430FgpRsKN4B7ae2BoU67+zeCdOVJ7AM+ufiej7SmYEWUC5rxxGjRGtbK3OGv7JYDSVJuUO8ng8WL58Od577z18++23aNmyZVDh+WShSYfFDh8+cwI5RaHK4LPPAKMR+OgjoE+fBruXuErg8DmQZkiDdusOaP44CO+Vo2mh9iZACAEBCUwonx3+BvdtfgIX5vXHG3/5Z6Df0P9eg+OeYrw1ZBH6Z/cGQPPIqDl1o06Pcm4PeDUgtsyDsUo5KETBuiPr4BJcuLjdxQH3z5pDa7Dn5B7c2PNGVq8WgEtwQZIlZJmzkGHMYL5/RpOI2mGxmphMJowcORLl5eU4fPgw9u7d25ThmicqFfDcc0B5Od0nuPNO4McfIak5SIoESZYgKRIEWQAv8fBInkBEkPWZF6Db8TPURcVwTQvfFTR140NYdfxHvP6Xf2JQTl8AQI4xCy7JjX2OYKX90dBXkKG3Ba0+Q5qMRBHa9GxoalgHKk6FC9vUPWA4pHAIhhQOCfHdND/8qRqMGiPybbSAO0vXwIgVYSkBvwXw7rvvYvXq1WjVqhWuu+46/Oc//4m0fEkNIYRO9JwM+cXnoL9nOpx3TYbbXQRJkSATGQpRaBQQp4ZGpYFZa4ZGrYHup43Q7fgZRK+H5/qrmyRHpehGueDAz6d/DSiB3hnd8P2lH6JDWmFQ3yxD0/zxXIzrCSc7XtELr+RFhiEDWeYsVrqREXNCdgdde+21WLFiBUwmE8aNG4fx48dj4MCB0ZIvJjTVHUTsdjrZ13jwMg9e4iErMiRFggIFHGj6BzVHHxq1pn5/LyHIuG4SdDt+gXvC9aj8x71Nen+7yn6DTBR0tXeELpqnQUURnI+n+YKYIjgjClHg5J3QcBpkm7NhM9iY758RUaLmDlKr1fjoo48wcuRIqGvVjd29eze6d0+dmq6CLKDcXQovHJAJzcBIQHWqWqWGhtNAq9bCoDXU+YJrf94F8ytLUbHwCbpXUANqBfwCotfDfetNIct1xHUc921+AtO7T8aAnN7okdE5/DcZApwohXU+INXw5/1J06ch25TNErUx4krISuDdd98Nel5ZWYn3338fr7/+OrZt25ZSIaOCLMDhc0BnzoFOrYNGq2mcL1cQYL/7fqhLSmGf9gAqXnymeuIkBJYXXwUAeK67CkpOdshyPb/n31hbsgnggA9zXg75+rARBMCWFtYJ6VRAIQpcgguEEOSac5FhzGCZURlxJ2z7c+3atZgwYQLy8/PxzDPPYOjQodi4cWMkZUsKOHAwaA3QqrWN38zT6VDx3HwQvR6GH9bBNvfxQL4d3boN1AowhGcFAMD0brfixg5jcX+PO8K6Plw4gFoCjDp4RS8qfBUwqA1oZWuFbHM2UwCMhCAkS6C4uBjLli3DG2+8AafTiXHjxoHneXz22Wfs9HCIiH17oeL5p2CfMhPG5SsgZ2TAdf9UKJkZ4C8YAKljeyjZWWGN3dKcjwXn/SPCEp8FSQJRq6NaTzgZEWQBbsENvVqPFtYWSNOnscmfkVA02hIYPXo0OnXqhF9++QWLFi3CiRMn8MILL0RTtmYPf9GFcD7+IADA8sZbMP37bUhdO6P8jZdQOfOekMeLZ7loThCpAmAbwgCqDuL5HPCKXmSZstDa3hrpxnSmABgJR6Mtga+//hr33HMP7rjjDnTs2DGaMqUU3r9fDtXpclifWYy0pxdBblcI/q9/CasYy72bH4NJY8TUrpOQbYxx7WFRBNIy2H4AaE1oXuKRpk9DhjGDpXxmJDSNtgTWrVuHyspK9O3bF/3798eLL76IU6dORVO2lMF9601w3zwevosuBD/gvLDGOOw6hg/+/AJv7P8gKOVDrOAI2w/gJR5l3rJAxs8WaS2YAmAkPCGfE3C73fjwww/x73//G5s3b4Ysy1i4cCEmTZoEq9V69gESkHDPCbjKilG0exOsea3O3vlsKAp9hFmOkRCCdSVbsPHkdtzX4/+aLk8oSBI4jw9ym5aAwRDbeycAsiKjUqiECipkmDKQbkhn1bkYcaex81qTcgft27cPb7zxBt5++21UVFTg4osvxhdffBHucHEjIZRAMuP1AhwHpU1rmiYjRSCEwC26IcoibAYbMowZtBQog5EANHZea9I3tlOnTliwYAGOHTuG999/vylDMZqArMT3bAYnSoDZnFIKwCf5UO4th5pTo5WtFVpYWzAFwEhKIvKtVavVGDNmTFJaAcnOrrLfMHDFGLz/x+dxk4FTlJTZD5AUCRXeCkiyhDxLHtrY2yBNn8YSvjGSltRZujVTXt//Po66T+DHks3xEUCWQTgViL55h4YSQlDJV8IluGAz2NDa3hpZ5iyW7pmR9LBPcJLzVL/Z6GbvhKEFg+IjgCiC6Jv3+QB/pk+L1oJMUyZL9cxoVjAlkOQYNQbc1vn6uN2fE0QQW1qz3A8QZAEuwQW9Wo98Sz7sBjs77MVodjAlkKR4JV+guHs84RQFirF5ZcGsGfKZZcpCuiGd1fhlNFuYEkhSbv6R1hh4rO9MdExrGx8hFIXuBzSTfEEKUeAW3JCJHDjtyyJ+GM0dpgSSkMOuY/ipZAsAQK+Koy9eEGjt4yTfDyCEwCt54ZN8zO/PSDmSypH71FNPgeM4TJs2LdDm8/kwZcoUZGZmwmKxYOzYsSgpKYmfkDGgjaUlfvrbcjzbfx5aW1rETQ5OFEGMBkCdvH5yn+RDua88kOqhla0VrHorUwCMlCFplMCWLVvwr3/9C+eee25Q+/Tp0/Hll1/i448/xpo1a3DixAn8/e9/j5OUsaO1pQXGtf1bfIWQFZAk3Q8QZAFl3jIa72/OQ2tba7bxy0hJkkIJuFwujB8/Hq+99hrS06sLoTscDrzxxhtYuHAhhg4dir59+2Lp0qVYv379GQvc8DwPp9MZ9EgGCCE4zZfHWwyKogAcl3TnA2RFRoWvAj7Rh2xTNtrY2yDLnMVy/TBSlqRQAlOmTMFll12G4cOHB7Vv27YNoigGtXfu3BmtW7fGhg0bGhxv/vz5sNlsgUerVsmR+2ddyRb0+/wyPPVLDEtGNoQo0pKYSbIfoBAFlXwlKoVK2PT0sFeuJZdF/TBSnoTfGP7ggw+wfft2bNmypc5rxcXF0Ol0sNvtQe25ubkoLi5ucMzZs2djxowZgedOpzMpFMFXR1fDJ/Nwia54i0L3A8zmhN8PYJu+DMaZSWglcPToUUydOhWrVq2CIYIpivV6PfRJmOtmfr8HcHGLC9E9/Zx4iwKIUsLvB/gkHzyiB0aNES3TWsKqszKfP4NRi4R2B23btg2lpaXo06cPNBoNNBoN1qxZg8WLF0Oj0SA3NxeCIKCioiLoupKSEuTl5cVH6CjCcRyGFQxGrjE7voIQQiuIJagiZZu+DEbjSWhLYNiwYdi1a1dQ280334zOnTtj1qxZaNWqFbRaLVavXo2xY8cCoDUOjhw5goEDB8ZD5KhwzF2EXGMWtKoE2bwURUCnS7hDYgpR4OSdUEGFbFM27AY78/kzGGchoZWA1WpF9+7dg9rMZjMyMzMD7bfccgtmzJiBjIwMpKWl4e6778bAgQMxYMCAeIgccQghmPTjTFSKLiwZNB+9MrvGWyRwgkBdQWFWQYsG/iRvNr0NmaZMdtKXwWgkifMtDpPnnnsOKpUKY8eOBc/zGDlyJF5+OQGiZyLEUfcJnPCUwCfzaGXJD30ASaKbt5HcCJVkuimcAChEgcPngFalRb4lH+nGdKi4hPZyMhgJRZPKSzYXEr28pEfy4ueyvRiY0ye0C3kenNcHAOAAEJUKRKuhoZ3hruIJAedwQmndCsQc39V2zdV/likLRm1ib1QzGLGksfNa0lsCqYBJYwxdAQDgPF6QnGwQkwkQBHA8D3g84Lw+cLIMAgAaDfXtazSNSwdddT4gnofEZEVGJV8JjUqDAksB7EY7W/0zGGHClECCohAFv5bvR4+MzuENwPOAVgvFaqFRPCYjnfQJoQpBEOm/Xi84XgC8PnCEgKhUVDH4LYZacIIIYjDEbT/Av/q3G+zIMmXBkADptBmMZIYpgQTlyyPf4v/Wz8bYwlF4ceBjIV/Pebwg2Vl1wzirQjv9NYEJQPcNRBGcIFZbCz4ecHuo0tBqQLRaqhREESQzo+lvMERkRYaTd0Kn1qGFtQVsBhtb/TMYEYApgQTlgPMgVJwK7axh7DcIArUC0qyN66/R0NW/sYa1UKUUAtaCjwd8bkCljvn5AP/qP92QjkxTJlv9MxgRhCmBBGVmj9txZZtLkGPMDPlazuMFycoMf7LmuKpzADoAZhCkA7JMFYIk0/TRMcC/+ter9WiZ1hJp+jS2+mcwIgxTAglM+7Q2oV8kCIBG03groLGo1YDfUogBHtEDn+RDuiEdWaYsduiLwYgSbFmVYOw8vQenfOGni+bcHhC7LWFTOpwNSZFQ5i0DCNAqrRUKrAVMATAYUYQpgQRCVETcsX42+n85GutLt4U+QKh7AQmGW3Cjkq9EhiEDrWytYDPYWLZPBiPKMHdQAlHqPQ2bLg2Vohu9MkJPD8G5PfVHBDURQghERYQoixAVEQCg4lTgwNF/OQ4cOHAcF9Tuf+1sSIoEJ+8MZPtM06exyZ/BiBFMCcQRhSg4WHkUvzsPYWTLIWhhzsPXI97CUfcJmDQhnn71WwFWS5PlkhUZgixAVERIigQOHHRqHbQqLewGOzhwkBQJMpHpQ5GhEAWKokAG/ZkQQv+t2kUgIODABf71Kwz//TKNmcg0ZUKnTo4iNQxGc4EpgRhRKbrwa/kBmDXGwAGwStGNC76i9ZD3jv0edh1dAYdTPJ5ze2hEUIh1F/yrfEmRIMgCFKJAzamhU+tg1Vlh0pqgU+uoEjhDCcaak35NJVBfGwGBpEhUaRAZVr0VVh0r7s5gxAOmBCIMIQSHXMfwa8V+DC+4AAY1dc28tu99/HPXK0GHv2w6KzqmtYVZY0IZXwG7rvF5i4IIISJIIQpd5Ve5djhw0Kq00Kl1SDOlwaAxBCb9UMIxOY6DmmP5+hmMZIMpgSbgFj3Y7tgLHVeKgbl9A+2jV92M03w5vh7xdiD1czf7OSgw5cKuC56o14z6uMkr4DNZAf4VviALIIQEXDtmrRlmnRlatTbg6mErcQYj9WBKoAn899h3mLR9Nvqld8OXl7wFgK6Ie2d2Q4n3FHwyH+g7osWFGNlySJ0xmjzxnsEK4CWellfUGpFpzAxa5bMqWwwGA2BKoEn0yO+FHGM28g3ZgM8XWIm/deGiOpN7tFbZnMdLc/nUYwW4RTdyzDnINmWzVT6DwagXpgSaQOeszvhu4vfI8AKqklIoVTl4YjbhCiKgVtdrBfgkH/RqPWx6FmvPYDAahh0WiwAk3Q4l3Q6u0kWTr8UIzlN1OrgeK8AjepBuSGenbRkMxhlhSiAScByUzAwQswmcyxWbe57BCvCKXhjUBqQZwow2YjAYKQNTApFCqwXJygIIRwu6RJmGrABCCLyiF+nGdHbwisFgnBWmBCIIMZugZGVA5fXR1MvR4kxWgOSFUWtEmp5ZAQwG4+wwJRBhSLodii0NnLMyavsDZ7ICfJIPGcaMM57uZTAYDD9MCUQalQpKViaI0QDO7Y78+H4roJ4cQR7RA6PGCKs+ObOIMhiM2MOUQDTQ6UCyswEF9DBXBOE8HhBbGmAMTjBHCAEv88gyZUGjYpG/DAajcTAlECWIxUz3BzzeyO0PCCKgqn8vwC26YdaaYdE1PYsog8FIHZgSiCKB/YHKyISNUivAWscK8CeFyzRlsnQQDAYjJJjfIJpU7Q+ofTw4lxvEYg5/LL8VYKsb9eMW3LDqrMwKYEQFRVEgRNityWg6Wq0WanXTF31MCUQbnQ5KdhZUx4vo/oAuvNh9zuulEUG1rABZkSEpEjKMGSGlfmYwGoMgCDh48CAURYm3KIx6sNvtyMvLa1JqGKYEYgCxWqBkZ0JdepLmF1KFOFmLIsCp6rcCRDfS9GnMCmBEHEIIioqKoFar0apVK6hC/dwyogYhBB6PB6WlpQCA/Pz8sMdiSiBGELsNitcLzllJV/QhwHnqtwIkRYKsyEg3prMkcYyII0kSPB4PCgoKYDKZ4i0OoxbGqvmgtLQUOTk5YbuGmGqPFWo1lOwsQKcD5/Y0/rozWAEuwQW7wQ6ztgl7DQxGA8hVUW26MF2YjOjjV86iKIY9BlMCsUSvp4pAkuhGbyPgPN56I4JEmZaGtBvszApgRBX2+UpcIvG3YUogxhCrBSQzAyq3GzjbZtuZ9gIEN2x6G0xaZqYzGIzwYUog1nAclIx0KGnWs54faMgKEGQBKk7F9gIYjBApLi7GxRdfDLPZDLvdHm9xEgKmBOKBWg0lKxPQasF5vPX3qbICSFoDVoDBBqPWWM+FDEZqM3HiRIwZM6be15577jkUFRVh586d2L9/PwCgsLAQHMeB4ziYTCb06NEDr7/+egwlji9MCcQLg4EqAlGkj1r4rQBiCp7oeYmHRqVBujE9VpIyGM2GP/74A3379kXHjh2Rk5MTaH/00UdRVFSE3bt344YbbsDkyZPx9ddfx1HS2MGUQBwhaVaQjHSo3J7gtNOiCICr3woQ3bAb7DBo6paUZDAYDVNYWIhPPvkEb731FjiOw8SJEwOvWa1W5OXloV27dpg1axYyMjKwatWq+AkbQ9g5gXhStT8AH0/PD1RtAHNuWi+gthXgk3zQqXSwG+xxEJbBoLgFmiLdpDUF9qQEWYAoi9CoNEF1rf19jVpj4ES7KIsQZAFqlTpoMdNQ30jVxtiyZQtuuukmpKWl4fnnnw/E2ddEURQsX74c5eXlKRMayyyBeKPRQMnOBDQawOut3guw1T1Q5hE9yDBmsOLxjLhimW+BZb4FpzynAm3//OmfsMy34K7/3hXUN+eZHFjmW3DEcSTQ9tKWl2CZb8EtX9wS1Lfw+UJY5luw9+TeQNuyncsiJnd2djb0ej2MRiPy8vJgq/EdmzVrFiwWC/R6Pa666iqkp6fj1ltvjdi9ExmmBBIBoxFKdhY4XgBX6aJ7AcZgdw8rHs9gRI/77rsPO3fuxHfffYf+/fvjueeeQ4cOHeItVkxg7qAEgaRZQXgenMNJrYAaoZ+EEHglL/It+ax4PCPuuGbT0OaaZ1TuG3wfpg2YVqegUelMmtumZiTblPOmYHKfyXXSnh+aeqhO34m9JkZS9AbJyspChw4d0KFDB3z88cfo0aMH+vXrh65du8bk/vGEWQKJgv/8QG5OHSvAJ/lg1LDi8YzEwKwzw6wzB51R0al1MOvMdVyV/r41M9xq1VqYdeY6wQ0N9Y01rVq1wjXXXIPZs2fH/N7xgFkCiYRGA1KrapjfCiiwFLDi8QxGI3E4HNi5c2dQW2ZmZqOvnzp1Krp3746tW7eiX79+EZYusWBKIMHxF49newEMRuP54Ycf0Lt376C2W265pYHedenatStGjBiBefPm4b///W+kxUsomBJIYPzF41umtWTF4xmMRrJs2TIsW7as0f0PHTpUb/vKlSsjI1CCw/YEEhh/8Xirrm5heQaDwYgETAkkKP7i8RnGDFY8nsFgRA2mBBIUt+CGRWuBVc+sAAaDET2YEkhA/MXjM02ZrHg8g8GIKmyGSUD8xePNOlY2ksFgRBemBBIMWZEDxeOZFcBgMKINm2USDJfggs1gY8XjGQxGTGDB5wmArMgQZAGSIoEQgnQDKxvJYDBiA1MCMUZWZIiKCFEWISkSAEDNqWk+Fa0ZFr2FFY9nMBgxI6HdQfPnz8d5550Hq9WKnJwcjBkzBvv27Qvq4/P5MGXKFGRmZsJisWDs2LEoKSmJk8TByIoMn+SDS3Ch3FuOMm8Z3IIbhBCYtWbkWfLQ2tYahemFaJveFi1tLWE32JkVwGBECVZovi4JrQTWrFmDKVOmYOPGjVi1ahVEUcSIESPgdrsDfaZPn44vv/wSH3/8MdasWYMTJ07g73//e8xlrW/CdwkuEEJg1BiRZ8lDG1sbFKYXotBeiJa2lsg0ZcKqt8KgMbBNYAYjQjSl0LzRaERhYSHGjRuH7777LoZSx4+EdgfVzt2xbNky5OTkYNu2bbjwwgvhcDjwxhtv4L333sPQoUMBAEuXLkWXLl2wceNGDBgwIOoyEhCUe8uh4lTQqrQwaozIMNDqX1q1FlqVlp34ZTAShJqF5mvy6KOPYvLkyRAEAYcOHcI777yD4cOH47HHHsOcOXPiJG1sSGglUBuHwwEAyMjIAABs27YNoihi+PDhgT6dO3dG69atsWHDhgaVAM/z4Hk+8NzpdIYlj06tQ4G1ABqVBjq1jk34jNSghiVeB7UaMBga11elAmrW+W2orzkykXKFhYU4fPgwAOCtt97ChAkTAonm/IXmAaB169a48MILkZ+fj3nz5uGqq65Cp06dIiJDIpI0PghFUTBt2jQMHjwY3bt3B0D9ezqdro5vLzc3F8XFxQ2ONX/+fNhstsCjVatWYcmkU+uQYcxAmj4NBo2BKQBGamCxNPwYOza4b05Ow30vvTS4b2Fh/f0ixJYtW3DJJZdg3LhxKCoqwvPPP3/G/lOnTgUhBJ9//nnEZEhEkkYJTJkyBbt378YHH3zQ5LFmz54Nh8MReBw9ejQCEjIYjETmTIXm6yMjIwM5OTkNpppuLiSFO+iuu+7CihUrsHbtWrRs2TLQnpeXB0EQUFFREWQNlJSUBEy7+tDr9dDr9Q2+zmAwzoDL1fBr6lrWcGlpw31VtdagCTjZEkKafbReQisBQgjuvvtuLF++HD/88APatm0b9Hrfvn2h1WqxevVqjK0yQ/ft24cjR45g4MCB8RCZwWj+hOKjj1bfGHD69GmcPHmyzrzT3EhoJTBlyhS89957+Pzzz2G1WgN+fpvNBqPRCJvNhltuuQUzZsxARkYG0tLScPfdd2PgwIExiQxiMBjNl+effx4qlarBcNPmQkIrgSVLlgAA/vrXvwa1L126FBMnTgRA435VKhXGjh0LnucxcuRIvPzyyzGWlMFgJBKhFpqvrKxEcXExRFHEwYMH8c477+D111/H/Pnz0aFDhyhLG18SWgkQQs7ax2Aw4KWXXsJLL70UA4kYDEYyEGqh+Xnz5mHevHnQ6XTIy8vDgAEDsHr1alx00UXRFjXuJLQSYDAYjFCJVKH5VCFpQkQZDAaDEXmYEmAwGIwUhikBBoPBSGGYEmAwGIwUhikBBoPBSGGYEmAwGIwUhikBBoPBSGGYEmAwGIwUhikBBoPBSGGYEmAwGCkDKzRfF6YEGAxGsyLUQvMAsGPHDlx99dXIzc2FwWBAx44dMXny5ECfQ4cOBYrR135s3LgxFm8rajAlwGAwUoaaheZzcnIAACtWrMCAAQPA8zzeffdd7N27F++88w5sNhvmzp0bdP23336LoqKioEffvn3j8VYiBksgx2AwGgUhgMcTn3ubTEBTC3zVV2j+5Zdfxs0334xRo0Zh+fLlgb5t27ZF//79UVFRETRGZmbmGasWJiNMCTAYjEbh8US07ntIuFxNLzy2ZcsW3HTTTUhLS8Pzzz8Po9GIb775BqdOncL9999f7zWpsG/AlACDwUgJaheaB4ADBw4AADp37tyoMQYNGgRVrdrIrjPVXE4CmBJgMBiNwmQ6c435aN87GjSmcFVNPvzwQ3Tp0iU6wsQJpgQYDEaj4LiEqwXfZM455xwAwG+//YaBAweetX+rVq2aXblJFh3EYDBSlhEjRiArKwsLFiyo9/XaG8PNEWYJMBiMZkdjC82bzWa8/vrruPrqq3H55ZfjnnvuQYcOHXDq1Cl89NFHOHLkCD744INA/9OnT6O4uDhoDLvdDoPBEJX3EQuYEmAwGM2OUArNX3HFFVi/fj3mz5+P66+/Hk6nE61atcLQoUPx+OOPB/UdPnx4nevff/99XHvttZETPsZwJNSdkWaI0+mEzWaDw+FAWlpavMVhMBICn8+HgwcPom3btkm90m3OnOlv1Nh5je0JMBgMRgrDlACDwWCkMEwJMBgMRgrDlACDwWCkMEwJMBgMRgrDlACDwWCkMEwJMBgMRgrDlACDwWCkMEwJMBgMRgrDlACDwWCEyMMPP4xevXo1i/swJcBgMJolR48exaRJk1BQUACdToc2bdpg6tSpOH36dEjjcByHzz77LKht5syZWL16dQSljR9MCTAYjGbHn3/+iX79+uHAgQN4//338fvvv+OVV17B6tWrMXDgQJSVlTVpfIvFUm9W0mSEKQEGgxESbrcbbrc7qCqXIAhwu93geb7evoqiBNpEUYTb7YbP52tU33CYMmUKdDod/ve//2HIkCFo3bo1Lr30Unz77bc4fvw45syZA4AWn3/sscdw3XXXwWw2o0WLFnjppZcC4xQWFgIArrzySnAcF3he200zceJEjBkzBk8++SRyc3Nht9vx6KOPQpIk3HfffcjIyEDLli2xdOnSIDlnzZqFc845ByaTCe3atcPcuXPDfs/hwpQAg8EICYvFAovFglOnTgXa/vnPf8JiseCuu+4K6puTkwOLxYIjR44E2l566SVYLJY6qZ0LCwthsViwd+/eQNuyZctClq+srAzffPMN7rzzThiNxqDX8vLyMH78eHz44YcBJfbPf/4TPXv2xI4dO/DAAw9g6tSpWLVqFQBanB4Ali5diqKiosDz+vjuu+9w4sQJrF27FgsXLsRDDz2Ev/3tb0hPT8emTZvwf//3f7j99ttx7NixwDVWqxXLli3Dnj178Pzzz+O1117Dc889F/J7bgpMCTAYjGbFgQMHQAhpsBZwly5dUF5ejpMnTwIABg8ejAceeADnnHMO7r77blx11VWBiTg7OxsALRyTl5cXeF4fGRkZWLx4MTp16oRJkyahU6dO8Hg8+Mc//oGOHTti9uzZ0Ol0WLduXeCaBx98EIMGDUJhYSFGjx6NmTNn4qOPPorUr6JRsKIyDAYjJFxV1eZNNaq/33fffZg2bRo0muAppbS0FACCVuRTpkzB5MmToVarg/oeOnSoTt+JEyeGLWdjS6XUri08cOBALFq0KOT7devWDSpV9bo6NzcX3bt3DzxXq9XIzMwM/E4AWrh+8eLF+OOPP+ByuSBJUsxrmjBLgMFghITZbIbZbAbHcYE2nU4Hs9kMvV5fb9+ak6NWq4XZbK5TBKWhvqHSoUMHcBwX5Faqyd69e5Genn7GVX041JaV47h62/x7Hhs2bMD48eMxatQorFixAjt27MCcOXMgCEJE5TobTAkwGIxmRWZmJi6++GK8/PLL8Hq9Qa8VFxfj3XffxTXXXBNQYhs3bgzqs3HjxiBXklarhSzLEZdz/fr1aNOmDebMmYN+/fqhY8eOOHz4cMTvczaYEmAwGM2OF198ETzPY+TIkVi7di2OHj2KlStX4uKLL0aLFi3wxBNPBPr+9NNPWLBgAfbv34+XXnoJH3/8MaZOnRp4vbCwEKtXr0ZxcTHKy8sjJmPHjh0Dhez/+OMPLF68GMuXL4/Y+I2FKQEGg9Hs6NixI7Zu3Yp27dph3LhxaN++PW677TZcdNFF2LBhAzIyMgJ97733XmzduhW9e/fG448/joULF2LkyJGB15999lmsWrUKrVq1qlO8vilcfvnlmD59Ou666y706tUL69evx9y5cyM2fmNhhebBCs0zGPWRCoXmCwsLMW3aNEybNi3eooQFKzTPYDAYjCbBlACDwWCkMOycAIPBSFn8ZxNSGWYJMBgMRgrDlACDwTgjLHYkcYnE34YpAQaDUS/+tA6xPsHKaDwejwdAeCer/bA9AQaDUS8ajQYmkwknT56EVqsNSufAiC+EEHg8HpSWlsJut9fJwxQKTAkwGIx64TgO+fn5OHjwYFzSGTDOjj+7aVNgSoDBYDSITqdDx44dmUsoAdFqtU2yAPwwJcBgMM6ISqVqtieGGc1oY/ill15CYWEhDAYD+vfvj82bN8dbJAaDwUh4moUS+PDDDzFjxgw89NBD2L59O3r27ImRI0cGFW9gMBgMRl2ahRJYuHAhJk+ejJtvvhldu3bFK6+8ApPJhH//+9/xFo3BYDASmqTfExAEAdu2bcPs2bMDbSqVCsOHD8eGDRvqvYbnefA8H3jucDgA0Kx7DAaD0Rzwz2dnO1CW9Erg1KlTkGUZubm5Qe25ubn47bff6r1m/vz5eOSRR+q0t2rVKioyMhgMRryorKyEzWZr8PWkVwLhMHv2bMyYMSPwXFEUlJWVITMzM6hu6tlwOp1o1aoVjh49GpE6BGy8xBovGmOy8dh4sRqTEILKykoUFBScsV/SK4GsrCyo1WqUlJQEtZeUlDR4iEKv19cpiG2328OWIS0tLaLFaNh4iTVeNMZk47HxYjHmmSwAP0m/MazT6dC3b1+sXr060KYoClavXo2BAwfGUTIGg8FIfJLeEgCAGTNmYMKECejXrx/OP/98LFq0CG63GzfffHO8RWMwGIyEplkogWuuuQYnT57EvHnzUFxcjF69emHlypV1NosjjV6vx0MPPVTHtcTGax7jRWNMNh4bLxHGrAkrNM9gMBgpTNLvCTAYDAYjfJgSYDAYjBSGKQEGg8FIYZgSYDAYjBSGKYEwiWTq6rVr12L06NEoKCgAx3H47LPPmiTb/Pnzcd5558FqtSInJwdjxozBvn37wh5vyZIlOPfccwOHVQYOHIivv/66STLW5KmnngLHcZg2bVpY1z/88MPgOC7o0blz5ybJdPz4cdxwww3IzMyE0WhEjx49sHXr1rDGKiwsrCMfx3GYMmVKWOPJsoy5c+eibdu2MBqNaN++PR577LEmFR2vrKzEtGnT0KZNGxiNRgwaNAhbtmxp9PVn+wwTQjBv3jzk5+fDaDRi+PDhOHDgQNjjffrppxgxYkTglP/OnTvDlk8URcyaNQs9evSA2WxGQUEBbrrpJpw4cSJs+R5++GF07twZZrMZ6enpGD58ODZt2hT2eDX5v//7P3Ach0WLFp3xPTcWpgTCINKpq91uN3r27ImXXnopIvKtWbMGU6ZMwcaNG7Fq1SqIoogRI0bA7XaHNV7Lli3x1FNPYdu2bdi6dSuGDh2KK664Ar/++muTZd2yZQv+9a9/4dxzz23SON26dUNRUVHgsW7durDHKi8vx+DBg6HVavH1119jz549ePbZZ5Genh7WeFu2bAmSbdWqVQCAq6++Oqzxnn76aSxZsgQvvvgi9u7di6effhoLFizACy+8ENZ4AHDrrbdi1apVePvtt7Fr1y6MGDECw4cPx/Hjxxt1/dk+wwsWLMDixYvxyiuvYNOmTTCbzRg5ciR8Pl9Y47ndblxwwQV4+umnmyyfx+PB9u3bMXfuXGzfvh2ffvop9u3bh8svvzys8QDgnHPOwYsvvohdu3Zh3bp1KCwsxIgRI3Dy5MmwxvOzfPlybNy48aypIEKCMELm/PPPJ1OmTAk8l2WZFBQUkPnz5zd5bABk+fLlTR6nJqWlpQQAWbNmTcTGTE9PJ6+//nqTxqisrCQdO3Ykq1atIkOGDCFTp04Na5yHHnqI9OzZs0my1GTWrFnkggsuiNh4tZk6dSpp3749URQlrOsvu+wyMmnSpKC2v//972T8+PFhjefxeIharSYrVqwIau/Tpw+ZM2dOyOPV/gwrikLy8vLIP//5z0BbRUUF0ev15P333w95vJocPHiQACA7duwIW7762Lx5MwFADh8+HJHxHA4HAUC+/fbbsMc7duwYadGiBdm9ezdp06YNee655846VmNglkCI+FNXDx8+PNB2ttTV8cafKjsjI6PJY8myjA8++ABut7vJaTmmTJmCyy67LOh3GS4HDhxAQUEB2rVrh/Hjx+PIkSNhj/XFF1+gX79+uPrqq5GTk4PevXvjtddea7KMAP38vPPOO5g0aVJIyQprMmjQIKxevRr79+8HAPz8889Yt24dLr300rDGkyQJsizXKSFpNBqbZFH5OXjwIIqLi4P+zjabDf3790/o7wzHcU3KKeZHEAS8+uqrsNls6NmzZ1hjKIqCG2+8Effddx+6devWZJlq0ixODMeScFJXxxNFUTBt2jQMHjwY3bt3D3ucXbt2YeDAgfD5fLBYLFi+fDm6du0a9ngffPABtm/fHpLfuSH69++PZcuWoVOnTigqKsIjjzyCv/zlL9i9ezesVmvI4/35559YsmQJZsyYgX/84x/YsmUL7rnnHuh0OkyYMKFJsn722WeoqKjAxIkTwx7jgQcegNPpROfOnaFWqyHLMp544gmMHz8+rPGsVisGDhyIxx57DF26dEFubi7ef/99bNiwAR06dAhbTj/FxcUAUO93xv9aIuHz+TBr1ixcd911TUoCt2LFClx77bXweDzIz8/HqlWrkJWVFdZYTz/9NDQaDe65556w5WkIpgSaOVOmTMHu3bubvKLr1KkTdu7cCYfDgf/85z+YMGEC1qxZE5YiOHr0KKZOnYpVq1ZFpIB5zRXwueeei/79+6NNmzb46KOPcMstt4Q8nqIo6NevH5588kkAQO/evbF792688sorTVYCb7zxBi699NIm+XQ/+ugjvPvuu3jvvffQrVs37Ny5E9OmTUNBQUHY8r399tuYNGkSWrRoAbVajT59+uC6667Dtm3bwpYzGRFFEePGjQMhBEuWLGnSWBdddBF27tyJU6dO4bXXXsO4ceOwadMm5OTkhDTOtm3b8Pzzz2P79u1hW49ngrmDQiSc1NXx4q677sKKFSvw/fffo2XLlk0aS6fToUOHDujbty/mz5+Pnj174vnnnw9rrG3btqG0tBR9+vSBRqOBRqPBmjVrsHjxYmg0Gsiy3CRZ7XY7zjnnHPz+++9hXZ+fn19HuXXp0qVJLiYAOHz4ML799lvceuutTRrnvvvuwwMPPIBrr70WPXr0wI033ojp06dj/vz5YY/Zvn17rFmzBi6XC0ePHsXmzZshiiLatWvXJFkBBL4Xif6d8SuAw4cPY9WqVU1OBW02m9GhQwcMGDAAb7zxBjQaDd54442Qx/nxxx9RWlqK1q1bB74vhw8fxr333ovCwsImyQgwJRAyyZC6mhCCu+66C8uXL8d3332Htm3bRvweiqIElegMhWHDhmHXrl3YuXNn4NGvXz+MHz8eO3fuhFqtbpJsLpcLf/zxB/Lz88O6fvDgwXVCavfv3482bdo0Sa6lS5ciJycHl112WZPG8Xg8UKmCv7pqtRqKojRpXIBOXPn5+SgvL8c333yDK664osljtm3bFnl5eUHfGafTiU2bNiXMd8avAA4cOIBvv/0WmZmZEb9HuN+ZG2+8Eb/88kvQ96WgoAD33XcfvvnmmybLxdxBYRDp1NUulyto1Xrw4EHs3LkTGRkZaN26dcjjTZkyBe+99x4+//xzWK3WgN/VZrPBaDSGPN7s2bNx6aWXonXr1qisrMR7772HH374IewPoNVqrbM/YTabkZmZGda+xcyZMzF69Gi0adMGJ06cwEMPPQS1Wo3rrrsuLPmmT5+OQYMG4cknn8S4ceOwefNmvPrqq3j11VfDGg+gE8DSpUsxYcIEaDRN+9qNHj0aTzzxBFq3bo1u3bphx44dWLhwISZNmhT2mN988w0IIejUqRN+//133HfffejcuXOjP9Nn+wxPmzYNjz/+ODp27Ii2bdti7ty5KCgowJgxY8Iar6ysDEeOHAnE8vuVdl5eXr3WxZnGy8/Px1VXXYXt27djxYoVkGU58J3JyMiATqcLabzMzEw88cQTuPzyy5Gfn49Tp07hpZdewvHjxxsMCz7b+62tlLRaLfLy8tCpU6d6xwuJiMQYpSAvvPACad26NdHpdOT8888nGzduDHus77//ngCo85gwYUJY49U3FgCydOnSsMabNGkSadOmDdHpdCQ7O5sMGzaM/O9//wtrrIZoSojoNddcQ/Lz84lOpyMtWrQg11xzDfn999+bJM+XX35JunfvTvR6PencuTN59dVXmzTeN998QwCQffv2NWkcQghxOp1k6tSppHXr1sRgMJB27dqROXPmEJ7nwx7zww8/JO3atSM6nY7k5eWRKVOmkIqKikZff7bPsKIoZO7cuSQ3N5fo9XoybNiwM/4uzjbe0qVL6339oYceCnk8f5hpfY/vv/8+5PG8Xi+58sorSUFBAdHpdCQ/P59cfvnlZPPmzWG/39pEMkSUpZJmMBiMFIbtCTAYDEYKw5QAg8FgpDBMCTAYDEYKw5QAg8FgpDBMCTAYDEYKw5QAg8FgpDBMCTAYDEYKw5QAg8FgpDBMCTAYNSgsLDxj2b5Dhw41qpxhNImnDJEof8pILJgSYCQU9dXirfl4+OGH4ypfq1atUFRU1KTaDM1BBkbzgSWQYyQURUVFgZ8//PBDzJs3Lyijp8ViiYdYAdRqddzTHyeCDIzmA7MEGAmFPwtkXl4ebDYbOI4Lavvggw/QpUsXGAwGdO7cGS+//HLQ9bNmzcI555wDk8mEdu3aYe7cuRBFMajPl19+ifPOOw8GgwFZWVm48sorg173eDyYNGkSrFYrWrduHZQ9tLYr5ocffgDHcVi9ejX69esHk8mEQYMG1UlF/fjjjyMnJwdWqxW33norHnjgAfTq1avB30N5eTnGjx+P7OxsGI1GdOzYEUuXLq1XBoCWxOzYsSMMBgMuuugivPnmm+A4DhUVFQCAZcuWwW6345tvvkGXLl1gsVhwySWXBCndLVu24OKLL0ZWVhZsNhuGDBmC7du3NyijIAi46667kJ+fD4PBgDZt2jSppgEjPjAlwEga3n33XcybNw9PPPEE9u7diyeffBJz587Fm2++GehjtVqxbNky7NmzB88//zxee+01PPfcc4HXv/rqK1x55ZUYNWoUduzYgdWrV+P8888Pus+zzz6Lfv36YceOHbjzzjtxxx131JnUazNnzhw8++yz2Lp1KzQaTVBa53fffRdPPPEEnn76aWzbtg2tW7c+a9WquXPnYs+ePfj666+xd+9eLFmypMHShAcPHsRVV12FMWPG4Oeff8btt9+OOXPm1Onn8XjwzDPP4O2338batWtx5MgRzJw5M/B6ZWUlJkyYgHXr1mHjxo3o2LEjRo0ahcrKynrvu3jxYnzxxRf46KOPsG/fPrz77rsRKXLCiDERyUXKYESBpUuXEpvNFnjevn178t577wX1eeyxx8jAgQMbHOOf//wn6du3b+D5wIEDyfjx4xvs36ZNG3LDDTcEniuKQnJycsiSJUsIISSQdnjHjh2EkOoUwN9++23gmq+++ooAIF6vlxBCSP/+/cmUKVOC7jN48GDSs2fPBuUYPXo0ufnmm+t9rbYMs2bNIt27dw/qM2fOHAKAlJeXE0KqUy/XTLH90ksvkdzc3AZlkGWZWK1W8uWXXwbaAJDly5cTQgi5++67ydChQ4miKA2OwUh8mCXASArcbjf++OMP3HLLLbBYLIHH448/jj/++CPQ78MPP8TgwYORl5cHi8WCBx98MKgs5M6dOzFs2LAz3uvcc88N/Ox3R5WWljb6Gn9FM/81+/btq2Nt1H5emzvuuAMffPABevXqhfvvvx/r169vsO++fftw3nnnnXV8k8mE9u3bB8lZ832VlJRg8uTJ6NixI2w2G9LS0uByuRosqzlx4kTs3LkTnTp1wj333IP//e9/Z3xPjMSEKQFGUuByuQAAr732WlCZvd27d2Pjxo0AgA0bNmD8+PEYNWoUVqxYgR07dmDOnDkQBCEwTmMqq2m12qDnHMedtXRjzWv8xcCbUu7x0ksvxeHDhzF9+nScOHECw4YNC3LdhEN974vUKCcyYcIE7Ny5E88//zzWr1+PnTt3IjMzM+j3V5M+ffrg4MGDeOyxx+D1ejFu3DhcddVVTZKREXuYEmAkBbm5uSgoKMCff/6JDh06BD38NZTXr1+PNm3aYM6cOejXrx86duyIw4cPB41z7rnnBtW6jQWdOnXCli1bgtpqP6+P7OxsTJgwAe+88w4WLVrUYHnLTp06YevWrSGPX5uffvoJ99xzD0aNGoVu3bpBr9fj1KlTZ7wmLS0N11xzDV577TV8+OGH+OSTT1BWVhbyvRnxg4WIMpKGRx55BPfccw9sNhsuueQS8DyPrVu3ory8HDNmzEDHjh1x5MgRfPDBBzjvvPPw1VdfYfny5UFjPPTQQxg2bBjat2+Pa6+9FpIk4b///S9mzZoVNbnvvvtuTJ48Gf369cOgQYPw4Ycf4pdffkG7du0avGbevHno27cvunXrBp7nsWLFCnTp0qXevrfffjsWLlyIWbNm4ZZbbsHOnTuxbNkyANVWSWPo2LEj3n77bfTr1w9OpxP33XffGS2nhQsXIj8/H71794ZKpcLHH3+MvLw82O32Rt+TEX+YJcBIGm699Va8/vrrWLp0KXr06IEhQ4Zg2bJlAUvg8ssvx/Tp03HXXXehV69eWL9+PebOnRs0xv+3d78sCgRxGMcfb0EtJkFMgsm2iLhVRBRfhAi+ALFome4W7brFZFDwTzBs0WrxLVgEuzYxXrt0BxdcjnO+nzzzm0nzwG8GplqtarVaabfbqVgsqlar6XQ6RbrvVqslY4wGg8FXC6XT6SiZTP44Jx6Pyxgj13VVqVTkOI6Wy+W3Y/P5vNbrtbbbrVzX1XQ6/XodlEgkfr3P2Wym+/2uUqmkdrutXq+nTCbz4/hUKqXRaKRyuSzP83S5XBSGoT4+OFb+E/4YBv5Ao9FQNpvVfD6PpL7v+wqCQNfrNZL6eB+0g4CIPR4PBUGgZrMpx3G0WCx0OBy03+9ftsZkMpHneUqn0zoejxqPx+p2uy+rj/dFCAARi8ViCsNQvu/r+XyqUChos9moXq+/bI3z+azhcKjb7aZcLqd+vy9jzMvq433RDgIAi3GDAwAWIwQAwGKEAABYjBAAAIsRAgBgMUIAACxGCACAxQgBALDYJ0XAX8FjhgBRAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "# Plot\n",
        "plt.figure(figsize=(4, 5))\n",
        "plt.fill_between(range(NSTEPS), rl_mean + rl_std, rl_mean - rl_std, color='g', alpha=0.1)\n",
        "plt.fill_between(range(NSTEPS), demo_mean + demo_std, demo_mean - demo_std, color='r', alpha=0.1)\n",
        "plt.fill_between(range(NSTEPS), expl_mean + expl_std, expl_mean - expl_std, color='b', alpha=0.1)\n",
        "plt.plot(range(NSTEPS), rl_mean, 'g:', label='LfR')\n",
        "plt.plot(range(NSTEPS), demo_mean, 'r--', label='LfD')\n",
        "plt.plot(range(NSTEPS), expl_mean, 'b-', label='LfCE')\n",
        "plt.plot(range(NSTEPS), ref_mean, 'k:', label='Optimal')\n",
        "#plt.xlabel('N. samples')\n",
        "#plt.ylabel(r'$||V_\\pi||_2$')\n",
        "plt.xlabel('Teaching signals')\n",
        "plt.ylabel('Average return')\n",
        "plt.ylim(0, 120)\n",
        "plt.xticks(range(0, 15))\n",
        "plt.legend(loc='lower right')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SW_VyigwJ76s",
        "outputId": "24669735-eb43-407a-f121-ed0202bbee24"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Optimal:\n",
            "[-2.27468901 -0.25581261  0.2578652 ]\n",
            "Contrastive Examples:\n",
            "[-0.56565847 -0.57387602  0.7746314 ]\n",
            "Demos:\n",
            "[-0.12874652 -0.00433172  0.06174637]\n",
            "RL:\n",
            "[-1.14462068  0.          0.        ]\n"
          ]
        }
      ],
      "source": [
        "print('Optimal:')\n",
        "print(wopt)\n",
        "print('Contrastive Examples:')\n",
        "print(w_expl)\n",
        "print('Demos:')\n",
        "print(w_demo)\n",
        "print('RL:')\n",
        "print(w_rl)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "wopt = [-2.27468901, -0.25581261,  0.2578652 ]\n",
        "w_expl = [-0.56565847, -0.57387602,  0.7746314 ]\n",
        "w_demos = [-0.2506747,   0.,          0.05851921]\n",
        "w_rl = [-1.14462068,  0.,          0.        ]"
      ],
      "metadata": {
        "id": "WUK9Tl3KzGxZ"
      },
      "execution_count": 26,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import scipy.stats\n",
        "from scipy.stats import f_oneway\n",
        "from scipy.stats import normaltest, ttest_rel, shapiro, mannwhitneyu, wilcoxon, ranksums, ttest_ind, kruskal, sem"
      ],
      "metadata": {
        "id": "LJcjetSjvpaJ"
      },
      "execution_count": 27,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.stats.mstats_basic import Ttest_indResult\n",
        "ttest_rel(perf_demo, perf_rl)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "t4eV-vhmh3tr",
        "outputId": "8041c752-b6d2-4a26-c7ab-279b1cd68129"
      },
      "execution_count": 28,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-28-5d2c048752d8>:1: DeprecationWarning: Please use `Ttest_indResult` from the `scipy.stats` namespace, the `scipy.stats.mstats_basic` namespace is deprecated.\n",
            "  from scipy.stats.mstats_basic import Ttest_indResult\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "TtestResult(statistic=array([ 0.9758447 ,  0.45076895,  0.60194464,  0.40784648,  2.00385341,\n",
              "        1.85868644,  1.49148567,  2.0435796 ,  1.48253231,  1.47923358,\n",
              "        0.86432665,  0.23613798, -0.18256843, -0.18166277, -0.58730324]), pvalue=array([0.33721191, 0.65550848, 0.55188853, 0.68638184, 0.05450633,\n",
              "       0.07324797, 0.14663286, 0.05017322, 0.14898519, 0.14985948,\n",
              "       0.39450171, 0.81498547, 0.85640697, 0.8571112 , 0.56154529]), df=array([29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29]))"
            ]
          },
          "metadata": {},
          "execution_count": 28
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import scipy.stats\n",
        "from scipy.stats import f_oneway\n",
        "from scipy.stats import normaltest, ttest_rel, shapiro, mannwhitneyu, wilcoxon, ranksums, ttest_ind, kruskal, sem"
      ],
      "metadata": {
        "id": "VQfsUdTyul_5"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ttest_rel(perf_demo[5], perf_rl[5])"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "a2xUe-m2e87z",
        "outputId": "7b6a32b4-e572-4966-94d7-ffddd35ed212"
      },
      "execution_count": 29,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "TtestResult(statistic=-0.34563342286309545, pvalue=0.734756805129781, df=14)"
            ]
          },
          "metadata": {},
          "execution_count": 29
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "mannwhitneyu(rl_mean,demo_mean)"
      ],
      "metadata": {
        "id": "KOxoPso6utXY",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "35dda6f9-0782-4c04-b948-aca9ece95058"
      },
      "execution_count": 30,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "MannwhitneyuResult(statistic=53.0, pvalue=0.014396852751179462)"
            ]
          },
          "metadata": {},
          "execution_count": 30
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "mannwhitneyu(rl_mean,expl_mean)"
      ],
      "metadata": {
        "id": "L2ltypYFu0f4",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "795ea5d1-50fc-4b25-d6c7-e1bf57b4aaee"
      },
      "execution_count": 31,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "MannwhitneyuResult(statistic=1.0, pvalue=4.092505485946723e-06)"
            ]
          },
          "metadata": {},
          "execution_count": 31
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "mannwhitneyu(demo_mean,expl_mean)"
      ],
      "metadata": {
        "id": "wXeM5X1gu6o4",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "7f118929-382d-491e-df76-2ea12d320056"
      },
      "execution_count": 32,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "MannwhitneyuResult(statistic=5.0, pvalue=8.969852672402501e-06)"
            ]
          },
          "metadata": {},
          "execution_count": 32
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "ttest_rel(rl_mean,demo_mean)"
      ],
      "metadata": {
        "id": "ueGaHcIQFivE",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "f9a66c63-b322-4ee1-cb50-5178b67f4a35"
      },
      "execution_count": 33,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "TtestResult(statistic=-4.242349397289592, pvalue=0.0008202274506022108, df=14)"
            ]
          },
          "metadata": {},
          "execution_count": 33
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "ttest_rel(rl_mean,expl_mean)"
      ],
      "metadata": {
        "id": "tYcAfUjWuFHJ",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "3a4be44c-f2ac-40b5-ba75-9f5598c356dd"
      },
      "execution_count": 34,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "TtestResult(statistic=-11.778996836073606, pvalue=1.1869481010475595e-08, df=14)"
            ]
          },
          "metadata": {},
          "execution_count": 34
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "ttest_rel(demo_mean,expl_mean)"
      ],
      "metadata": {
        "id": "a5EWFCCrvOhQ",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "86a9f702-641c-40db-f406-19f82783da5c"
      },
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "TtestResult(statistic=-8.574141999189148, pvalue=6.056556742978775e-07, df=14)"
            ]
          },
          "metadata": {},
          "execution_count": 35
        }
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "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.7.4"
    },
    "latex_envs": {
      "LaTeX_envs_menu_present": true,
      "autoclose": false,
      "autocomplete": true,
      "bibliofile": "biblio.bib",
      "cite_by": "apalike",
      "current_citInitial": 1,
      "eqLabelWithNumbers": true,
      "eqNumInitial": 1,
      "hotkeys": {
        "equation": "Ctrl-E",
        "itemize": "Ctrl-I"
      },
      "labels_anchors": false,
      "latex_user_defs": false,
      "report_style_numbering": false,
      "user_envs_cfg": false
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "41b50f7bfd3b480a952593c96da0980b": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "HBoxModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "HBoxModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "HBoxView",
            "box_style": "",
            "children": [
              "IPY_MODEL_1131e727887c4030ad451a44f11870e9",
              "IPY_MODEL_581f0a2815d44294a7084bfb2ed3ba95",
              "IPY_MODEL_1764f7a1af18499c8ec767b06914f439"
            ],
            "layout": "IPY_MODEL_021806461ab44c5abdd4629d6d268199"
          }
        },
        "1131e727887c4030ad451a44f11870e9": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "HTMLModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "HTMLModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "HTMLView",
            "description": "",
            "description_tooltip": null,
            "layout": "IPY_MODEL_e633f0851f984feea42c0f676462a47a",
            "placeholder": "​",
            "style": "IPY_MODEL_e79702d34d7b426eb7f6612532abb0e1",
            "value": "100%"
          }
        },
        "581f0a2815d44294a7084bfb2ed3ba95": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "FloatProgressModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "FloatProgressModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "ProgressView",
            "bar_style": "success",
            "description": "",
            "description_tooltip": null,
            "layout": "IPY_MODEL_86248af124514eb995033fcd604c8fab",
            "max": 30,
            "min": 0,
            "orientation": "horizontal",
            "style": "IPY_MODEL_39a5ec80a5234a7dbcdf7615df616518",
            "value": 30
          }
        },
        "1764f7a1af18499c8ec767b06914f439": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "HTMLModel",
          "model_module_version": "1.5.0",
          "state": {
            "_dom_classes": [],
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "HTMLModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/controls",
            "_view_module_version": "1.5.0",
            "_view_name": "HTMLView",
            "description": "",
            "description_tooltip": null,
            "layout": "IPY_MODEL_faade32c915f4120895baada845e1c90",
            "placeholder": "​",
            "style": "IPY_MODEL_3274dd63da9747959237c72b1ab8a738",
            "value": " 30/30 [06:42&lt;00:00, 12.81s/it]"
          }
        },
        "021806461ab44c5abdd4629d6d268199": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "e633f0851f984feea42c0f676462a47a": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "e79702d34d7b426eb7f6612532abb0e1": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "DescriptionStyleModel",
          "model_module_version": "1.5.0",
          "state": {
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "DescriptionStyleModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "StyleView",
            "description_width": ""
          }
        },
        "86248af124514eb995033fcd604c8fab": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "39a5ec80a5234a7dbcdf7615df616518": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "ProgressStyleModel",
          "model_module_version": "1.5.0",
          "state": {
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "ProgressStyleModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "StyleView",
            "bar_color": null,
            "description_width": ""
          }
        },
        "faade32c915f4120895baada845e1c90": {
          "model_module": "@jupyter-widgets/base",
          "model_name": "LayoutModel",
          "model_module_version": "1.2.0",
          "state": {
            "_model_module": "@jupyter-widgets/base",
            "_model_module_version": "1.2.0",
            "_model_name": "LayoutModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "LayoutView",
            "align_content": null,
            "align_items": null,
            "align_self": null,
            "border": null,
            "bottom": null,
            "display": null,
            "flex": null,
            "flex_flow": null,
            "grid_area": null,
            "grid_auto_columns": null,
            "grid_auto_flow": null,
            "grid_auto_rows": null,
            "grid_column": null,
            "grid_gap": null,
            "grid_row": null,
            "grid_template_areas": null,
            "grid_template_columns": null,
            "grid_template_rows": null,
            "height": null,
            "justify_content": null,
            "justify_items": null,
            "left": null,
            "margin": null,
            "max_height": null,
            "max_width": null,
            "min_height": null,
            "min_width": null,
            "object_fit": null,
            "object_position": null,
            "order": null,
            "overflow": null,
            "overflow_x": null,
            "overflow_y": null,
            "padding": null,
            "right": null,
            "top": null,
            "visibility": null,
            "width": null
          }
        },
        "3274dd63da9747959237c72b1ab8a738": {
          "model_module": "@jupyter-widgets/controls",
          "model_name": "DescriptionStyleModel",
          "model_module_version": "1.5.0",
          "state": {
            "_model_module": "@jupyter-widgets/controls",
            "_model_module_version": "1.5.0",
            "_model_name": "DescriptionStyleModel",
            "_view_count": null,
            "_view_module": "@jupyter-widgets/base",
            "_view_module_version": "1.2.0",
            "_view_name": "StyleView",
            "description_width": ""
          }
        }
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}