kernel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import vonmises"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptm_torsion = {\n",
    "    \"SEP\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"OG\"],\n",
    "        [\"CA\", \"CB\", \"OG\", \"P\"],\n",
    "        [\"CB\", \"OG\", \"P\", \"O1P\"]\n",
    "        ],\n",
    "    \"TPO\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"OG1\"],\n",
    "        [\"CA\", \"CB\", \"OG1\", \"P\"],\n",
    "        [\"CB\", \"OG1\", \"P\", \"O1P\"]\n",
    "        ],\n",
    "    \"PTR\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"CG\"],\n",
    "        [\"CA\", \"CB\", \"CG\", \"CD1\"],\n",
    "        [\"CE1\", \"CZ\", \"OH\", \"P\"],\n",
    "        [\"CZ\", \"OH\", \"P\", \"O1P\"]\n",
    "        ],\n",
    "    \"M3L\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"CG\"],\n",
    "        [\"CA\", \"CB\", \"CG\", \"CD\"],\n",
    "        [\"CB\", \"CG\", \"CD\", \"CE\"],\n",
    "        [\"CG\", \"CD\", \"CE\", \"NZ\"],\n",
    "        [\"CD\", \"CE\", \"NZ\", \"CM1\"]\n",
    "        ],\n",
    "    \"ALY\": [\n",
    "        [\"N\", \"CA\", \"CB\", \"CG\"],\n",
    "        [\"CA\", \"CB\", \"CG\", \"CD\"],\n",
    "        [\"CB\", \"CG\", \"CD\", \"CE\"],\n",
    "        [\"CG\", \"CD\", \"CE\", \"NZ\"],\n",
    "        [\"CD\", \"CE\", \"NZ\", \"CH\"]],\n",
    "    }\n",
    "\n",
    "free_rotamer = {\n",
    "    \"SEP\": \"chi3\",\n",
    "    \"TPO\": \"chi3\",\n",
    "    \"PTR\": \"chi4\",\n",
    "    \"M3L\": \"chi5\",\n",
    "}\n",
    "\n",
    "nchis = {k: len(v) for k, v in ptm_torsion.items()}\n",
    "# rotamer bin definitions for special cases; sp3 by default\n",
    "rot_bin_lib = {\n",
    "    (\"PTR\", 1): 4,\n",
    "    (\"ALY\", 4): 6,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read ptm dataframes\n",
    "ptm_data = []\n",
    "for ptm in ptm_torsion:\n",
    "    ptm_data.append(pd.read_csv(f\"data/{ptm}_df.csv\"))\n",
    "ptm_data = pd.concat(ptm_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "backbone independent library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "header = \"#AA     Prob   \"\n",
    "line_formatter = \"{: <3}  {:8.4f}\" #\"{:8.1f}\"\"{:8.2f}\"\n",
    "\n",
    "def format_bb_indp_table(lines, residue):\n",
    "    n_ang = (len(lines[0]) - 1) // 2\n",
    "    formatter = line_formatter + \"{:8.1f}\" * n_ang + \"{:8.2f}\" * n_ang\n",
    "    table = header + \" \".join([f\"chi{i+1}Val\" for i in range(n_ang)]) + \" \" + \" \".join([f\"chi{i+1}Sig\" for i in range(n_ang)]) + \"\\n\"\n",
    "    #table = \"\"\n",
    "    for l in lines:\n",
    "        table += formatter.format(residue, *l) + \"\\n\"\n",
    "    return table\n",
    "\n",
    "def rotameric_class(ang, n_splits=3):\n",
    "    if n_splits == 3:\n",
    "        if 0 < ang <= 120:\n",
    "            return 0\n",
    "        elif -120 < ang <= 0:\n",
    "            return 2\n",
    "        return 1\n",
    "    elif n_splits == 4:\n",
    "        if 45 < ang <= 135:\n",
    "            return 0 \n",
    "        elif -123 < ang <= -45:\n",
    "            return 2\n",
    "        elif ang > 135 or ang <= -135:\n",
    "            return 1\n",
    "        return 3\n",
    "    elif n_splits == 6:\n",
    "        if -30 < ang <= 30:\n",
    "            return 0 \n",
    "        elif 30 < ang <= 90:\n",
    "            return 1\n",
    "        elif 90 < ang <= 150:\n",
    "            return 2\n",
    "        elif ang > 150 or ang <= -150:\n",
    "            return 3\n",
    "        elif -150 < ang <= -90:\n",
    "            return 4\n",
    "        return 5\n",
    "\n",
    "def weighted_circmean(angles, weights):\n",
    "    \"\"\"\n",
    "    Compute the weighted circular mean of angular data.\n",
    "\n",
    "    Parameters:\n",
    "    angles (array-like): The angular data (in radians).\n",
    "    weights (array-like): The weights for each angle.\n",
    "\n",
    "    Returns:\n",
    "    float: The weighted circular mean of the angular data.\n",
    "    \"\"\"\n",
    "    # Convert angles to complex numbers\n",
    "    complex_numbers = np.exp(1j * np.array(angles))\n",
    "\n",
    "    # Compute weighted sum of complex representations\n",
    "    weighted_sum = np.sum(weights * complex_numbers)\n",
    "\n",
    "    # Compute the angle of the weighted sum\n",
    "    mean_angle = np.angle(weighted_sum)\n",
    "\n",
    "    return mean_angle\n",
    "\n",
    "def weighted_circstd(angles, weights, mean_angle):\n",
    "    \"\"\"\n",
    "    Compute the weighted circular standard deviation of angular data.\n",
    "\n",
    "    Parameters:\n",
    "    angles (array-like): The angular data (in radians).\n",
    "    weights (array-like): The weights for each angle.\n",
    "    mean_angle (float): The weighted circular mean of the angular data.\n",
    "\n",
    "    Returns:\n",
    "    float: The weighted circular standard deviation of the angular data.\n",
    "    \"\"\"\n",
    "    # Calculate angular differences from the mean, accounting for circularity\n",
    "    angle_diffs = np.angle(np.exp(1j * (np.array(angles) - mean_angle)))\n",
    "\n",
    "    # Calculate weighted variance\n",
    "    weighted_variance = np.sum(weights * angle_diffs**2) / np.sum(weights)\n",
    "\n",
    "    # Circular standard deviation is the square root of the variance\n",
    "    circ_std = np.sqrt(weighted_variance)\n",
    "\n",
    "    return circ_std\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For backbone independent rotameric torsion regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_rotameric_bb_indp(ptm, data, last_free_rot=False):\n",
    "    results = []\n",
    "    nchi = nchis[ptm]\n",
    "    if last_free_rot:\n",
    "        f_class = data[f\"chi{nchi}\"].apply(rotameric_class)\n",
    "        f_stats = np.zeros((3, 3))\n",
    "        for l in range(3):\n",
    "            ang_rad = data[f_class==l][f\"chi{nchi}\"].apply(np.radians)\n",
    "            kappa, loc, _ = vonmises.fit(ang_rad, fscale=1)\n",
    "            weights = vonmises(loc=loc, kappa=kappa).pdf(ang_rad)\n",
    "            wm = weighted_circmean(ang_rad, weights)\n",
    "            f_stats[l, 0] = len(ang_rad)/len(data)\n",
    "            f_stats[l, 1] = np.degrees(wm)\n",
    "            f_stats[l, 2] = np.degrees(weighted_circstd(ang_rad, weights, wm))\n",
    "        nchi -= 1\n",
    "\n",
    "    #bins = product(range(3), repeat=nchi)\n",
    "    chi_classes = []\n",
    "    for n in range(nchi):\n",
    "        chi_classes.append(data[f\"chi{n+1}\"].apply(lambda x: rotameric_class(x, 6 if ptm==\"ALY\" and n==4 else 3)).values)\n",
    "    chi_classes = np.array(chi_classes)\n",
    "\n",
    "    for c in np.unique(chi_classes.T, axis=0):\n",
    "        mask = np.ones(len(data), dtype=bool)\n",
    "        for i, l in enumerate(c):\n",
    "            mask &= chi_classes[i] == l\n",
    "        if np.sum(mask) == 0:\n",
    "            continue\n",
    "\n",
    "        angles = data[mask]\n",
    "        if np.sum(mask) == 1:\n",
    "            line = np.concatenate(([1/len(data)], angles[[f\"chi{n+1}\" for n in range(nchis[ptm])]].values[0], np.zeros(nchis[ptm])))\n",
    "        else:\n",
    "            line = np.zeros(2*nchis[ptm] + 1)\n",
    "            prob = np.sum(mask) / len(data)\n",
    "            for n in range(nchi):\n",
    "                ang_rad = angles[f\"chi{n+1}\"].apply(np.radians)\n",
    "                kappa, loc, _ = vonmises.fit(ang_rad, fscale=1)\n",
    "                weights = vonmises(loc=loc, kappa=kappa).pdf(ang_rad)\n",
    "                wm = weighted_circmean(ang_rad, weights)\n",
    "                line[1+n] = np.degrees(wm)\n",
    "                line[1+n+nchis[ptm]] = np.degrees(weighted_circstd(ang_rad, weights, wm))\n",
    "            if last_free_rot:\n",
    "                for j in range(3):\n",
    "                    line[0] = prob * f_stats[j, 0]\n",
    "                    line[-1-nchis[ptm]] = f_stats[j, 1]\n",
    "                    line[-1] = f_stats[j, 2]\n",
    "                    results.append(line.copy())\n",
    "                continue\n",
    "            line[0] = prob\n",
    "        results.append(line)\n",
    "        \n",
    "    return results\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = calc_rotameric_bb_indp(\"ALY\", ptm_data[ptm_data.residue==\"ALY\"], last_free_rot=False)\n",
    "t = format_bb_indp_table(r, \"ALY\")\n",
    "with open(\"aly_bi.lib\", \"w\") as f:\n",
    "    f.write(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For non-rotameric torsions, use DBSCAN to cluster before regression (SEP as an example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAFfCAYAAACr9rwWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcrklEQVR4nO29eXxTVfr4/05Smi50pZC0yCbigkX21VFkgIILjt9xHWcc3BiXYRQQZXFGwJHN0ao/cRs/juC+zOg4KgOlDuCCQmVRFhdQoCxdWNqmBdpCcn9/5CZkuUlu0iw3yXnzyuuS+5x7cpKnee7Jc57zPDpJkiQEAoFAkJDoYz0AgUAgEEQOYeQFAoEggRFGXiAQCBIYYeQFAoEggRFGXiAQCBIYYeQFAoEggRFGXiAQCBKYlFgPIBrYbDYOHjxIVlYWOp0u1sMRCASCNiNJEo2NjRQVFaHX+56vJ4WRP3jwIF26dIn1MAQCgSDs7Nu3jzPOOMOnPCmMfFZWFmD/MLKzs9vUl80mUf1TPcctrWRkp2LumYtef/rXwcGddXz0zLdteg2AK/54AYBiX+nH23HWzyZ2nVnDiYyTALTPNXLDQ0PdxhIKDStWUv3ggwHbmefPB1DVNv+22yi44w/oDAYkq5WfJ1zJqZqaoMZ1xgvPkzloUFDXxIRa4B3gOqBTdF7y+W+eZ+n2pWHt84oeV/DA4AdITUlVlJftLuOhLx8K2M/cYXMZf+b4sI4tZGKgm1Bw2JgDO+vYXLbPZ7vm1uP85fUbnPbNF0lh5B0umuzs7DYZ+Z821/LZ2zs5Vt/iPJeZa+Si63vRs7/9r6Z9/ywKOu13axMKepsRgPTUTC9ZQUN7xlX0p77LBg6nNgFgPQ7Hamx0PievTa9r6NaNJoMhYLv8bt0AVLVtXbqU2hUrMM2eRUpOLmmHD4OK61zJOH6izTfoqLALeBAYD5wVejdWm5VNtZs4dPwQHTM6MqDTAAx65c8sPTMdQ3pwn2cg/lv9X/774X/5RdEvuLXPrV6v383UTdVrLvl+Cfl5+YzpNias4wuJMOkmknjaGKXvvyeBXNBi4VUlP22uZcUL27yM97H6Fla8sI2fNtcCoNfruOj6Xm1+vcxsI5nZxqCuOWZp240FIGPQQFLMZvD1h6PTkWI2kzFoYOC2LpyqqeHAvVNo/N//QhpXSseOIV0Xj5TvLWfcv8Zx68pbmfHZDG5deSvj/jWO8r3liu2HFA6J2Fg+P/i54uvXNdepur6upY5pa6b5HLvgNL5sTFsRRl4FNpvEZ2/v9Nvm83d2YrPZc7317N+J8XcUk5nrbqTb5xkZN+l8r/OetM8zUtgrl8JeuQHbuhLsTUEJncGAafYs+YmH8Zafm2bPQmcw+G/riZwHz/Lhh0GPyXFTSQbK95Yzbc00ao67u7Nqj9f6NJaDTIPISc2J6Lhqjtc4X99qs/JoxaNBXb94w2KsNmuERhf/qLExoSKMvAqqdtYHvLs21bVQtbPe+bxn/078fsEIrpran7G39eaqqf25af4IzhpoCjjT/8V1vdDrdUH9KnDcGMJBdkkJnZ96khSTye18islE56eeJLukJGBbRSQJ69Gj6PPyVM3+AdDpnDeVRMdqs7JowyIkvBPDSvK/RRsWeRlLg97A3BFzIz4+CYnFGxZTUV3hdRMKdF318Wo21W6K4OjiGzU2JlSEkVeBWjeIZzu9Xkfnc/I4e7CZzufkORdF/c30x99R7PTvu7Y1Zp5ePmlNtbK72yFaU09/2R03hnCRXVLCWZ+U03XZMooee4yuy5Zx1iflbgbes22HO+9U1XfOlRPs/wlg6FPMZq+biubJASbIxyDZVLspoPGsOV7Di9++6HV+TLcxPHHJE+Sm5gb/wkFQfbyaipqKkK49dPxQmEcTJG3QTaQJh6vVF0mx8NpW1LpBgnGX9OzfiR59O9rv4JYWMrPtM3ElQ+1o+/V/d/PtJ/uxcILl47cC9hvDL67r5XZjCBc6g4HMoer8vTqDgczhwzny/PMB22b9cjQZAwdSs2Ahp6qrvfvKzqbDxN9TcOed8TeD7wn8J7RL1RrBZ755hrPyzvJazBzTbQwXd76Yi9++mGOnjoU2CDWEWIGiY0aM11XaoJtIEw5Xqy+EkVeBwzfu7+dUKO4Sx0xfbdshl5/JoEt7ULWjnuaDraQVpVLYW/nGEAscC7GnamqcPng3dDpSTCYyBg20h1PabBycMtWrmdTYyOElz2Ds1Su+ZvEAJ4F6IBdoF9ylwRjBxRsWM6rLKK+Im9SUVH7f+/c89+1zwb14EAwpHMIHP31A7fFaRdeSJzp0mDJMDOg0IGJjUkUbdBNp1NiYUBHuGhWo8Y2H213ibyydW/PoOc5E59Y8zRh4CG7RVrJaqV20WLkj+QZRs2AhkjXOFuu2Yo/B3hr8pQM6DcCUoWJtA/z6uAeaI7dInZOawyDTIGYOmQnYDbg/HPIZQ2b4DAGNGm3QTaQJV1SeYt8R6TUBCcaPnsz4Wog15OVR9Phj6LOyqH3yKQ7OnKXoqnEiSZyqrub41xsjPGLtYNAbnMZTDav2rKKiusJrIfbIiSPhHpqT4UXDMegNjOk2htJLSumU4f53r9e5mxRThonSS0q1ESevcXzZmLYi3DVBEIwfPV6RrFaOf72RU4cOkdKxo9O1Ekz77JISJJuN6nkPY6uzx1Nbjx7l4H3Tld04fjh1KMaLdVFmTLcx/LHfH3lmyzMB2775w5u8+cObmDJMzBwy02lII+n73lyzGavN6jT0o7qMctu01begL98c/sbvJq5gNnolG0o25rillbKXtofcpzDyQRKMHz3esJSVeS2GppjNmGbPUvSN+2qfffllHP3Hy94GPYSa8Ya8XI6t36D6ppMITOoziXd/eJfaE7Wq2jti6B0zZofbR63P3MEQ0xA21Gzw26bmRA2bajcx2DwYsP/6cPzfgedzV8r2lPHIV49Q13J6M5XnTSrZUbIx+hSd12779ipn/DpJCuGbF2dYLBZycnJoaGiIj63xgdgEDAQ2AmFay7KUlXHg3inehlj2pXuGMvptH8Y/KV1uLlJ9vfO5v5uOJgiTbhybooCgFjdXXL0Cg97gvF7NteYMMzOGzKDV2sqMz2YEbL/4osVcduZlgd+EB6Vfl/Ly9pd9jj/ibp0IfG+iic0muc3wM0168vJyA9o1MZP3wPOD1KQ7pi/QAAROa6EKyWqlZsFCZeMsSaDTUbNgIVmjRzsXTf22DyOuBh5Op0dAq/HzYdKNw+e9aMMiVRuPXDccOWbS2anZNLQ2uLXLNebyl2F/IS8tz8tdUlGtLv690lIZ9Psp21Pm08A7xu8rYihshPl7E208Z/gWi0XVdcLIu6AmAZkmMABh/EFy/OuNqhdBM4cOCdw+kijcdDRFGHXj6vNetWcVb/7wZsBrDh0/5HcWX99Sj16nV3SpONw8gW4q//zxn0y6YJJqY2y1WXnkq0cCtvO8SYWdMH9v4gURXSOjNgGZJtgJjJOPYUDt4qajXcwXQ7UceRNm3Th83mO7j1XVvkN6B5+pEcDuFvGXR2Z44fCAr+Hwy6tlU+0mNx+8PyK6KzbMuokXhJEn+ARkMacRKJOPQSBZrRxbv4GGjz7m2PoNzhh0tRkeHe20khEy5jcbJULUTSAcs2xfcek6dJgzzEiS5Hcm7iuPjCPz5b9/+req8QRjjINpG9FdsRHSjdYR7hqCS0AWr5E1/iJnskaPVr1TFVx2tkbIZaPPysLWGPibqJWbTTRwxNBPWzMNHTq3mbrrhqOjzUdV9edqeINZpHUQjDFW2zY/LT/2u2ITEDGTJ/QEZPGCIxLG0yg7c7x/8onqnargsbM1zBjy8znr07Wqc9onE742ILluOFJrUB3t/GW+9IUpPbgUBWp38j449EERLx8BhJEnMgnItIKaSBjHIqba9MJg39la8Kc/hX285of+giE9PaibTjIxptsYVl69kn+M+weLL1rMP8b9gxVXr3CGHgYyqA63jsNIq8l86UmzrZnV+1arbu/4FeIvBcIt599CSXcNRku1AatN4sufjvDBlgN8+dMRrDFy9wp3DZFLQBYxugBL5GMAgomcyS4pIWv0aJ87Xj13t3aYdDt177yD1Ve9Vp0OXXo60vHjqt9azaLFoNfbbypPPentYjKZtB0nH4RuQkVpA5KD1ftW03yqWVGmlEcmlIVOS4vFbfOVGnyFhOYZ8/jzsD9Hx8BHQTcOVmyrYt6HO6hqOK2Lwpw05kzozfjiwsgPwIWIbob69NNP+dvf/sbGjRupqqri/fff56qrrnLKb775ZpYtW+Z2zdChQ/nqq6+cz1taWpg+fTpvvvkmJ06cYPTo0Tz77LN+q5N7omYzlCO6xhfxmp+m4aOPOTh9esB2RY89Rs4Vl/uUB9zdCu6/FuRNUfr27bE1NakfsMfmq2DTLCQzgXzrucZc5gyf42aYK6oruHXlrUG/lufmK7UkQ0qDFduquOu1TV5acPyOee53A8Ji6NVu8oyou+bYsWP07duXJUuW+Gwzfvx4qqqqnI/ly5e7yadMmcL777/PW2+9xeeff05TUxNXXHEF1jBnJ4yrBGRHgdfkYwCCjZxRwp9P/+g/Xib/1lsU3TwFf5ocnIEHrwyUjpz2OVdcTubQIdo38EHoJpyo8a0b9UZGdRnldi5Q1I4vQq325PgVctmZlzHYPDi6Bj4KurHaJOZ9uENRC45z8z7cEVXXTUTdNZdeeimXXnqp3zZGoxGz2awoa2ho4KWXXuLVV19lzBj77OO1116jS5culJeXM27cuLCON24SkO0BbsK+PTvff9NgcrwroWY3rOXj5fQsW8mJzVvcZtyW/64I8o2d7td181VcsQfVugnnrFZVVSmPvDPgP2pHDTGv9hQMe1Ctm1DZsPuom4vGEwmoamhmw+6jDO/ZITKD8CDmC69r1qyhU6dOnH322UyaNIna2tObjjZu3MjJkycpcfG/FhUVUVxczLp163z22dLSgsVicXuoxVfJvnglmBzvSqj26W/0ntG1NcRRk3HwYcIRl37ryluZ8dkMbl15K+P+NU6xULca1BpbpXa+onbUEPNqTxqjttG3gQ+lXTiIqZG/9NJLef311/nf//7H448/TkVFBb/85S9pabEvgFZXV5OamkpenntsuslkotqP4Vm4cCE5OTnOR5cuUVhp0TDBFOb2RK2hPTBlCpUTJ3Jw+nQqJ05k1+gxnKo76j8UMgCJGgfv8J17zrwd2SRDMfTBhk564hq184cL/qCqLxHX7k2nrLSwtgsHMY2uuf76653/Ly4uZtCgQXTr1o2PP/6YX//61z6vkyQJnR/DMWvWLKZNm+Z8brFYhKEPEDnjC7WG1tbgngjrVE0NB6dOI//WW+wLs8FkpwzgQopn/PnOJSRn2oFgE3UFSi+spgSfw1+u9lfB5T0uT7hF07YypEc+hTlpVDc0Kzq+dIA5J40hPSLkL1Ig5u4aVwoLC+nWrRs7d9pTDJjNZlpbW6mrc897UVtbi8nkOxbYaDSSnZ3t9kgoMoFhBJ1NL5RFTIdPP+jZuGzQLR8vp/MTpRhyc4O6Nm7j4APoJpDvvC0Lmr5K8gVbgk/tr4JRXUcFbqQlQvzeBINBr2POhN4AXkvZjudzJvTGEEU3sKaM/JEjR9i3bx+FhfbwooEDB9KuXTtWrVrlbFNVVcW2bdsYMWJErIYZe84BvpSPEcavTz8Qsr/ekJdPp1nqd8jqg7khaI0AummL7zwQanbEqkFNxI3rhqq4IUrfm/HFhTz3uwGYc9xdMuactLCFTwZDRN01TU1N7Nq1y/l89+7dbNmyhfz8fPLz85k7dy5XX301hYWF7Nmzh9mzZ1NQUMD/+3//D4CcnBxuu+027rvvPjp06EB+fj7Tp0+nT58+zmgbQeTxtTFJn5Pj5aZR4tShQ7Tz88vLE1tDg7ZzxrcBtbPknxp+oqK6IuiIG6WSfMH2oTZPjnDV+GZ8cSFje5vZsPsotY3NdMqyu2iiOYN3ENHNUGvWrGHUKO+fdBMnTuS5557jqquuYvPmzdTX11NYWMioUaP461//6uY/b25u5v777+eNN95w2wwVjI9dVIYKD54bkySblX23BN5I0+XlfwBwYMpUVTcFwOmXP+uT8vhy2wTQjdVmZdy/xqkuzRfL0njle8u9dqg6qkjFZam+OK8M5YlauybK/8UjGvhjlaxWjm3Y4N9w63QYcnLAaPSd+iAAXZcti69YeRW6Cba0H8Af+/2RSX3UF+oIhNoY/YTaoaqB7004UWvXRO4aQdAopTjwQo6msXqU7wuWRIyVD7a0H8AzW57h3R/eZdbQWW2eRSvN0H39YvCXJ0cQH2hq4VWgfXylOPAkxWQKywJqosbKu8Wl91EXl157IrQ4eqvNSkV1Bct/Xs7zW55n6pqpYY3RF2gbMZMXqMZvigNX0tNJHzaUxn9/EPqLJXCsvINg49IdBBNHrzRrV6ItMfoCbSNm8vFIb+x1KntH92VVF/A+cSI4A59IOeND0E0wqQGCiaP3tbM2HH3HJTH63sQaYeTjkTTgLPkYBRy1YS1lZWHvu+BPk32mW8gaPVqxJq2asQZzTVgJQTehZIIMNPsPpeKT2r7jlih/b7SCcNfEI7uBvwB/BXpE9qVULbKGguyOKbjzTgruvNMZmmno0AGQaFqzhqo5c7G57HZ21KT1FTvvr45t1OLtQ9CNa1y6WgLN/kOp+KS277glit8bLSFm8hpCdbmwOuB1+RhB1C6yBo2HO8aRbkGX2o6qWbPYd8ut1C17xc3Aw+matEq/KALVsY3ErxBFQtSNI+ImUC1Uz/J9vghlNq6277glSt8brSFm8hpBS+XCIIhF1hBwlPBzuGNOHTpEy+7dHHnmmQCDsuevd9SkdS1LGCjnvec1WsSxW/XFb1/kmW+8P4tgdpsGOxsXO1kTFzGT1wCOcmGexQaqG5q567VNrNhWFfUxqV5kDRJDfj49y1aCzcbOX1zkTE8c0MA7kPPhHH31Nae/PZg6tlrHoDdwZ787eeKSJ7xm9cHkoAnWzx9sfhtB/CBm8jEmULkwHfZyYWN7m6Oa9yJSm5CsR49ycOZMGpf/t0391C5axNGlSzHNnoXUelLVNfG0saqtOWjU5J+5u9/ddM3qGv87WQV+EUY+xoRULqwQmCMfI0QkNyG11cA7cPjbCyb/UVX7qGysCqNu2rrb1NfOWlOGKX7zz7SFKHxvtIgw8jEmpHJhhcDciAzHidrasJ1mzqB24SJ7u2gjj+vIslfQ5+V5LdQ6iebGqijoJhjCkZUyYdCYbqKFMPIxJqRyYRbsebGHAxHKt+bII3/g3ineVZ1comOyS0rIHjuW419vpLWqitr587E1NkZmUD6QLBbf0eDR3lgVBd0Ei8g/I6NB3UQDsfAaY4b0yMec7dvQ67BH2biVC9sFjJePEURtbVhHCGTeVb+icP4jdsMaYl3XSJB/6y3Ri5OPkm4EIRDvurFZYfdnsPWf9qNN3UY/MZOPMat2VNN8SllZsSoX5kqwtWF9FRiJGZLE0Zf+QXrfvglXgESQROz4D6yYAZaDp8+lqivEI4x8DHGETvpyNeRmtGPhr/vEJE7eFcdMXS2eN4aWXbs48vzzERxhYOIhTl4gUGTHf+Cd34OnpWhUtw4m3DUxwl/opANjip6xvc1RG1Nb8MwZA5wuGj58eIxHR9zEyQsEbtis9hm8zyDrwIiZfIwIFDoJUG1pcQ+ddGAEespHDRAoZ0zGoIHoc3OxtbGASFuJSpy8xnQjcCEedbN3nbuLJgSSy8hvAdq7PM/DnqioGdih0N6RwuMH4JiHrDuQDxwC9nnIsoBegBX4RqHfPvaQyK51GWS1uKugJquZw5mtZDen0KU+g5YNVnBU10sHzgPOB94FWrCXNHNwntxmL3DE4zVNQGegEXu6VVfa2ccEwFbAc29RL/k9HQA8fiE2bv+UAwumoLOmktbqksN1Lxy+4wV4wf40pSYfvVTkdm1ru/3YDBYMpzrQ7pS7f9Gqb+Rk6j6QUkhrOdtNpktP44RtC+hstGvthsGW6SY/mVKFNaUOvTWH1JOdnedTj3Sz6/Ic+YRSRt3e2LMU7sY7x0mh/LDgvXhnxK6X84H38NbN2dj/9vYDtR7XFgBdgePA9x4yPdBP/v8O7H+rrpwJ5ALVgKctyJXlrcA2vOkn9/8j0OQh6yqP6zBQ6SFrL78fG/bvlCfFQCrwM1DvISsCzPL5nz1kaZxOA7xF7t+Vc4EMeTyHPWSdgDPk9/GjhywFuAC7bj7AWzdnYY+2qZIfrsTQRtAO+K4Jqvp6XHsQ2h+CEznYv9ABkJKAhoYGCZAaaJAkpNOP38oNdkru5x0PB8MUZK/KsiUKshLHC/vot1aS1u06LJWdVe0le3jUdqnbjI+ku371tfd1/V3GlKrQ7zZZdpuCbKYsW60g6+zSb2cF+WpZNtNbVm/+SNpxzrnSru5XeMmsuhbph4tHSj9efLF0wrjdS76v6F5pxznnStUdF3rJLJmfSDvOOVf6oedwxc/w+16DpB3nnCs1ZnzmJavq9LC045xzpf3m+72vHebyXpV0s1OW/VZBNkeWrVCQ9XTpt0BBvk6WTVWQ3S3LNirIslz67a0g/0CWLVCQXSPL9vl4r82yfKSC7EVZ9qKCbKQsa/bR7z5Zfo2CbIEs+0BB1tvlvWYpyDfKsrsVZFNl2ToFWYFLvz0V5Ctk2RwFWQxthCRJkjT6sMK1syRpTrbU8Ksb7HatoUHyR3IV8l7bQHZ7lwDZGN6lrQaJ6x/4kubDVjfPmmMmn9OcQr9Tefzj5sGnI2scM/lvgZHA8/LrOIjyTP7E1m3sf/heTrU7iM5mxNja0+0yCYmWtO8ASG05E73kHioa6kzePG8eluoVHH35JdUz+Y73T6f98OGQSWRn8t8ClwDP4a4bMZO3E8uZvK/vjZZn8jut8PSV0OIyY5dn8pb6bHKe2h+wkHdyGfkAH0a0cUTXgPsSiiNY8rnfDVCOrNFI1fmGjz7m4PTp0XtBeedqz7KV/FQyTlWIpiE3F/PD86IXPqkR3QgUCJNurDaJDbuPUtvYTKcs+x6WiIU47/gPvHOTosjSIpGzqDGgXUsun7zGGF9cyHO/G+CVYtgcwxTDwRDVItsuO1dPbN6iysB3nDGDDr+/SYRNCsJGVFOCOyNr2oYw8jFmfHEhY3ubozczCCMB89uEEV1mJoWP/JXskhIaPvpY1TXtOnYUBl4QNnzta3GkBPf5yztUwhBZA8LIawKDXucdJhkH+M1vE2akpiZqFy1Gp9er/gUR1V8agoQmJinBm8KT9E9shopHzgbWyccY4yu/ja59ex9XhM6p6moO3DsFa91RUsxm3/lxdDpSzOboZJ30REO6EXjQBt0EkxI8bLRXl7YgEMLIxyPtsWfSC78dDYnskhLO+qScrsuWkTdxIoa8PKQmz3CN8FGzaDGdZsq+Sk9DH+2sk55oTDcCF9qgm5BSgreVY0dA589Eq/vFIIx8PLIfmCYfNYLOYMDaUE/dK69g9ZXXPRxI9lJ+KXn5qjJkRh0N6kYg0wbdhJQSvC3s+A/882aQPONIg0f45OORWuAJ4HfYY4M1QCQLfytx6tAhcq64PKgMmVFBg7oRyLRBN0N65FOYk0Z1Q7OiX16HPSrOLSV4qPjNV+N4QQNctQQW/TZgd8LIC8JCpAp/+8KxqBpshkyBIBQMeh1zJvTmrtc2oUN5X0tIKcFtVnsUTWMVHDsEmR3tx0BRNZIV0tXdUISRF4SFqBXJjmYpP4HAhbDva1HKER8MTeomVcLIC8JCNMMV1SyqSlarttw4goQgbPtafOWID4byeaqaRXTh9dNPP2XChAkUFRWh0+n497//7SaXJIm5c+dSVFREeno6l1xyCdu3b3dr09LSwp/+9CcKCgrIzMzkyiuvZP/+JF/VKgDulo8awbExKqJl//R6VaX8LGVl7Bo9hsqJEzk4fTqVEyeya/QYLGVlkRubAw3qRiATJt049rX8ql9nhvfsEJqLJpDPXQ0nPBNUKRNRI3/s2DH69u3LkiVLFOWPPvoopaWlLFmyhIqKCsxmM2PHjqXRpRD0lClTeP/993nrrbf4/PPPaWpq4oorrsBqVVffMCHpCjwjHzWCY2OU/YlCWKNOhz43t20vYrNx9B8v+zXWlrIyDtw7xWt94FRNDQfunRJ5Q69B3QhktKKbMO1kVUvUEpTpdDref/99rrrqKsA+iy8qKmLKlCnMmGGPeW5pacFkMrF48WLuuOMOGhoa6NixI6+++irXX389AAcPHqRLly4sX76ccePGqXptrSYoCxlHtkJHVj4N4a+ACGDfHduWPznZJ3/WJ+Ve7hfJamXX6DG+F4D9XBs2NKybpEcrutn6T/jXbW3uRm2CspjFye/evZvq6mpKXH56G41GRo4cybp16wDYuHEjJ0+edGtTVFREcXGxs40SLS0tWCwWt0dC8T32bHqeaWk1gOvGqKLHHqPrsmWc9Uk52SUlp3fHmt1LGurz8si75Rav84rIcfJKpfwCRvj4uTZsaFg3SY9WdHPYM9d3ZInZwmu1/GU0eWxmMZlM7N2719kmNTWVvLw8rzbVfr7MCxcuZN48dYsSgvDjL6zRs8i366Koafp9HHp6iaqi30rRPGojfKIWCSQQeGKzwoYXo/qSMd/xqvPw30qS5HXOk0BtZs2aRUNDg/Oxb59nxn5BLHHcBHKuuJzMoUOcrhOdwaC66LdSNI9IXCbQPHvXqV4wDRcxM/Jm+ae554y8trbWObs3m820trZS57FN3rWNEkajkezsbLeHID4IGKXjJ/lYW64VCKJCmDJLBkPMjHyPHj0wm82sWrXKea61tZW1a9cyYsQIAAYOHEi7du3c2lRVVbFt2zZnm5CwWWH3Z/YFkN2f2Z/HE3rs5cNi/jss/ASM0sF3nHxbrg0bCaybuCdWunG1NzEw8hH1yTc1NbFr1+limLt372bLli3k5+fTtWtXpkyZwoIFC+jVqxe9evViwYIFZGRkcOONNwKQk5PDbbfdxn333UeHDh3Iz89n+vTp9OnThzFjxgQ/oG3vQe16+OG/cMLl10F2EYxfDL2vbOtbjg79sNcZTVCyS0rgqSe9o3RMJkyzZ/mNk2/LtWGhHwmtm7imH9HXjdKuVp0+LInH1BLREMo1a9YwatQor/MTJ05k6dKlSJLEvHnzeOGFF6irq2Po0KE888wzFBcXO9s2Nzdz//3388Ybb3DixAlGjx7Ns88+S5cuXVSPwxlCOTOLbKPST3n53HWvxI+hTwLasmtV7HgVxJxw7Gr1g9oQyuQq5O3TyAPo7DP6KVtBr3FjsAO4FniX09XtBdpA6Ea7RFo3jmRjTTWQUQAf3OV/01MbZ/SikHfQSGA5YFdSj4tiPRj/NGP/gw1jfQJBmBC60S6R1E0oycYkG3QdAZW+9/yEA2HkPYnBwohAIIhj2uKWibCBBxED4E2Y6ioKBIIkIFzJxiKImMm7klEA3XyEZrr629qb7O207rsXCASRJcrJxkJBGHlXLrhO2XAr+dtiGXZ5JvCBfBRoC6Eb7RIJ3cSBe1cYeVfOucz7nC9/m6XKfj4WYZe5gIj0jHqYpNUmBS4WkYvQjVbJJfy6iQP3rjDyDrI7e7tq/PrbJEAHK2bCuZdH13VTDbwM3AKoSNyYiPhLaRyJDU8rtlV5lX0rVCr7JnSjXSKhm24j7L/qLVUo2wkdZBXC/3vevuv1s7+F6YXVIxZeHRRf7W2oA/rbXMIuo8lBYLZ8TEKiXRhkxbYq7nptk5uBB6huaOau1zaxYlvV6ZNJrhtNEwnd6A12ty1wuqQ37s8vXQxnjoRRs+w3BK92kUUYeQfrnra7ZlxR62+LA79coiBZrdQsWKhceEQ+V7NgIVKYKodZbRLzPtzh87ccwLwPd2C1aTe6QhBhel9pd9tmexTyzi5yd+f6vSFEDuGucSJ5u17U+tviwC+XKARTGMRXTvtg2LD7qNcM3u3lgKqGZjbsPsrwnh3a/HqCOKX3lXbbESgCz3FDCHbjVBsQM3lXPF0vDn+bz7uuTtmXL4gY0S4MUtuobnuk2naCBEZvsO+W73ON/ehrna73lTBlG4xbEJ1hReVV4olGF/+qGn/b+EXRj5fPBa6Rj0lGtAuDdMpKC65dLkmrG82Ti3Z0ozdEzQMgjLwnxzxmgGr9bdHkTOxJlpIwFjvahUGG9MinMCfN3285CnPs4ZRAUutG82hNN1Ey8sIn70mmwgxQrb8tWrQCtUAnIDU2Q4gVjsIgB+6dYjf0rguwESgMYtDrmDOhN3e9tgkd7kFyDsM/Z0Lv0/HySawbzaM13TjDLyPrmxczeU+O/KR8Xq2/LRpsA7rIxyQku6SEzk89SYpHCcgUk4nOTz0Z9jj58cWFPPe7AZhz3F035pw0nvvdAPc4+STXjabRmm6c7uDIRtqImbwnaxfBwc0w4k8iP42GyS4pIWv06KjteB1fXMjY3ubAO14FgmCIQrSNMPJK7Fxpf8RbWcAkQ2cwhCVMUi0GvU6ESQrCj6s7+PuPYP3zYe1euGv84chP47lJSiAQCMKJwx186WIYcU94uw5rbwmHvMy2YqY9j41AIBBEmpK/wjXLICU9LN0JIx+QGOWn8Uc/7CXM+sV2GAIF+iF0o1X6ET+6Kb4KZlaGxdALI68WLeWn0QNGhPa0iNCNdok33aSkwlVt98/Hy9sNDxP+P3v1p1DQUn6aH4FL5KNAWwjdaJd41E3xVTB8cpu6SC4j394Mxw8HeZEG89M0AWvlo0BbCN1ol3jVzbj5MPxPIV+eXEb+WG2QF8QwP41AIBA4GPcIXLssJE9EcsXJZ3YKrn12kd3Aizh5gUAQa86/Cs6bcDq9itQeFl0a8LLkMvJdhwYo1YX9Tjl+ob1kl9jxKhAItIQjnh7AYlF3SQSHoz0Cpg7WwRVPwAXXxT4/jT+6Ai/KR4G2ELrRLkmqG50kKdVRSywsFgs5OTk0NDSQnZ1t38HqmSsiu7NwzQgEgrjBy675ILncNQ60ljo4WA4D/wauAkKMCBVECKEb7ZKkuklOIw/uvq14oxKYBAwgqf5Y4wKhG+2SpLpJLp+8QCAQJBnCyAsEAkECE3MjP3fuXHQ6ndvDbDY75ZIkMXfuXIqKikhPT+eSSy5h+/btMRyxQCAQxA8xN/IA559/PlVVVc7H1q1bnbJHH32U0tJSlixZQkVFBWazmbFjx9LY2BjDEceY9sBI+SjQFkI32iVJdaOJhdeUlBS32bsDSZJ48sknefDBB/n1r38NwLJlyzCZTLzxxhvccccd0R6qNjgbWBPrQQgUEbrRLkmqG03M5Hfu3ElRURE9evTghhtu4OeffwZg9+7dVFdXU+JSmNloNDJy5EjWrfOd372lpQWLxeL2SChsQIt8FGgLoRvtkqS6ibmRHzp0KK+88gorV67kxRdfpLq6mhEjRnDkyBGqq6sBMJnc0/yaTCanTImFCxeSk5PjfHTp0iWi7yHqbAHS5KNAW2xB6EarbCEpdRNzI3/ppZdy9dVX06dPH8aMGcPHH38M2N0yDnQ69xQEkiR5nXNl1qxZNDQ0OB/79u2LzOAFAoFA48TcyHuSmZlJnz592Llzp9NP7zlrr62t9Zrdu2I0GsnOznZ7CAQCQTKiOSPf0tLCd999R2FhIT169MBsNrNq1SqnvLW1lbVr1zJihIaKeAgEAoFGiXl0zfTp05kwYQJdu3altraWRx55BIvFwsSJE9HpdEyZMoUFCxbQq1cvevXqxYIFC8jIyODGG2+M9dAFAoFA88TcyO/fv5/f/OY3HD58mI4dOzJs2DC++uorunXrBsADDzzAiRMnuPvuu6mrq2Po0KGUlZWRlZUV45HHkGJgHxBkDRRBFBC60S5JqpvkTDUsEAgEcY5au6Y5n7xABT8D18pHgbYQutEuSaobYeTjkXrgn/JRoC3qEbrRKvUkpW6EkRcIBIIERhh5gUAgSGCEkRcIBIIERhj5eKQIWCAfBdpC6Ea7JKluRAilQCAQxCEihDKRqQf+Q9JFCcQF9QjdaJV6klI3wsjHIz8DvyLp4n3jAqEb7ZKkuhFGXiAQCBIYYeQFAoEggRFGXiAQCBIYYeTjkTSgt3wUaAuhG+2SpLqJeaphQQj0BrbHehACRYRutEuS6kbM5AUCgSCBEUY+HtkCZJN0Vefjgi0I3WiVLSSlboSRj0dsQKN8FGgLoRvtkqS6EUZeIBAIEhhh5AUCgSCBEUZeIBAIEhhh5OORc4GN8lGgLYRutEuS6kbEyccjGcCAWA9CoIjQjXZJUt2ImXw8Ugn8UT4KtIXQjXZJUt0IIx+PHAaelY8CbSF0o12SVDfCyAsEAkECI4y8QCAQJDDCyAsEAkECI4x8PNIJmCofBdpC6Ea7JKludJIkSbEeRKRRW9VcIBAI4gW1dk3M5OORJuBL+SjQFkI32iVJdSOMfDzyIzBCPgq0hdCNdklS3QgjLxAIBAlM3Bj5Z599lh49epCWlsbAgQP57LPPYj0kgUAg0DxxYeTffvttpkyZwoMPPsjmzZu56KKLuPTSS6msTLL9yQKBQBAkcZGgrLS0lNtuu43bb78dgCeffJKVK1fy3HPPsXDhQq/2LS0ttLS0OJ9bLBb7f7YA7V0a5gE9gGZgh8ILO5IZ/QAc85B1B/KBQ8A+D1kW0AuwAt8o9NsHaAf8BDR4yDoDJqAO2O0hSwfOw661XGCnh/w8uc1e4IiHzCT33ahwXTt5TABbgZMe8l7yezoA1HjIOgDdgBPAdx4yHdBf/v93chtXemDXQY3ctys5QE95LFvxpi9gkN9Lo4esC9AROArs8ZBlAufI/9+k0G9vIA37Z1/nISuUHxZgl4fMCJyPb92cjf1vbz9Q6yErALoCx4HvPWR6oJ/8/x3Y/1ZdOVN+vWrgoIcsV5a3Atvwpp/c/494L0Z2lcd1GO9cL+3l92NDuZReMZAK/AzUe8iKALN8/mcPWRr2zx+5X88KTudiTzJWiXdqgk7AGfL78PS5pwAX4Fs3Z2EvC1glP1zRuo1Qg6RxWlpaJIPBIL333ntu5++55x7p4osvVrxmzpw5EuD1aKBBkpBOP34rX7BTcj/veDgYpiB7VZYtUZCVyLIGH/3WyvIJCrLHZdk7CrL+LmNKVZBvk2W3KchmyrLVCrLOLv12VpCvlmUzFWS3ybJtCrJUl377K8jfkWWPK8gmyLJaH59hgywvUZAtkWWvKsiGuYxJqd+dsuy3CrI5smyFgqynS78FCvJ1smyqguxuWbZRQZbl0m9vBfkHsmyBguwaWbbPx3ttluUjFWQvyrIXFWQjZVmzj373yfJrFGQLZNkHCrLeLu81S0G+UZbdrSCbKsvWKcgKXPrtqSBfIcvmKMg0bCMaljbY7VqD48ugjObj5A8ePEjnzp354osvGDFihPP8ggULWLZsGT/88IPXNUoz+S5dutCwtoHs9i7xpFq/S/uayQNsxn7rckXM5O3EciYP8C1wykMuZvJ2YjmTB9gOtHjI43Qmb9lrIad74Dj5uDHy69atY/jw4c7z8+fP59VXX+X77z2/Fd4k3Gao7cCvgA84bVgE2kDoRrskmG4SZjNUQUEBBoOB6upqt/O1tbWYTKYYjSrGtGC/w3vOSASxR+hGuySpbjRv5FNTUxk4cCCrVq1yO79q1So3941AIBAIvImL6Jpp06Zx0003MWjQIIYPH87f//53KisrufPOO2M9NIFAINA0cWHkr7/+eo4cOcLDDz9MVVUVxcXFLF++nG7dusV6aAKBQKBpNL/wGg4SbuHVgj3R0nDsUQEC7SB0o10STDdq7VpczOQFHmQD42I9CIEiQjfaJUl1o/mFV4ECVcBcvGN6BbFH6Ea7JKluhJGPR6qAeSTdH2tcIHSjXZJUN8LICwQCQQIjjLxAIBAkMMLICwQCQQIjjHw8kgf8Vj4KtIXQjXZJUt2IOHmBQCCIQxImQZlAgWbsqW49084KYo/QjXZJUt0IIx+P7MCei1opv7UgtgjdaJck1Y0w8gKBQJDACCMvEAgECYww8gKBQJDACCMvEAgECYzIQhmPDMC7iLdAGwjdaJck1Y2YyQsEAkECI4x8PPID9sIHP8R6IAIvhG60S5LqRhj5eOQY8JV8FGgLoRvtkqS6EUZeIBAIEhhh5AUCgSCBEUZeIBAIEhhh5OOR7sCr8lGgLbojdKNVupOUuhFx8vFIPvC7WA9CoIjQjXZJUt2ImXw8cgh4Rj4KtIXQjXZJUt0IIx+P7AMmy0eBthC60S5Jqhth5AUCgSCBEUZeIBAIEhhh5AUCgSCBEUY+HskCSuSjQFsI3WiXJNWNCKGMR3oBK2M9CIEiQjfaJUl1I2by8YgVsMhHgbYQutEuSaqbmBr57t27o9Pp3B4zZ850a1NZWcmECRPIzMykoKCAe+65h9bW1hiNWCN8A+TIR4G2ELrRLkmqm5i7ax5++GEmTZrkfN6+fXvn/61WK5dffjkdO3bk888/58iRI0ycOBFJknj66adjMVyBQCCIK2Ju5LOysjCbzYqysrIyduzYwb59+ygqKgLg8ccf5+abb2b+/PlkZ2dHc6gCgUAQd8TcJ7948WI6dOhAv379mD9/vpsr5ssvv6S4uNhp4AHGjRtHS0sLGzdu9NlnS0sLFovF7SEQCATJSExn8vfeey8DBgwgLy+PDRs2MGvWLHbv3s3//d//AVBdXY3JZHK7Ji8vj9TUVKqrq332u3DhQubNmxfRsQsEAkE8EPaZ/Ny5c70WUz0fX3/9NQBTp05l5MiRXHDBBdx+++08//zzvPTSSxw5csTZn06n83oNSZIUzzuYNWsWDQ0Nzse+fQmWrKIPUCsfBdpC6Ea7JKluwj6Tnzx5MjfccIPfNt27d1c8P2zYMAB27dpFhw4dMJvNrF+/3q1NXV0dJ0+e9Jrhu2I0GjEajcENPJ5oB3SM9SAEigjdaJck1U3YjXxBQQEFBQUhXbt582YACgsLARg+fDjz58+nqqrKea6srAyj0cjAgQPDM+B45CdgKvAE0DNwc5vNyoHvttNUX0f73Dw6n3c+er0hwoNMUoLUjSCKJKluYuaT//LLL/nqq68YNWoUOTk5VFRUMHXqVK688kq6du0KQElJCb179+amm27ib3/7G0ePHmX69OlMmjQpuSNrGoAPgbmBm+5cv47/Lf07TUcPO8+1zy/glzf/gV5DR0RqhMlLELoRRJkk1U3MomuMRiNvv/02l1xyCb179+ahhx5i0qRJvPnmm842BoOBjz/+mLS0NC688EKuu+46rrrqKh577LFYDTuu2Ll+Hf8pXeBm4AGajh7mP6UL2Ll+XYxGJhAIokXMZvIDBgzgq6++Ctiua9eufPTRR1EYUWJhs1n539K/+22zetnf6Tl4qHDdCAQJTMzj5AWR4cB3271m8J40HjnMge+2R2lEAoEgFggjH490Bh6Xjz5oqq9T1ZXadgKVqNCNIEYkqW5intZAEAImYJr/Ju1z81R1pbadQCUqdCOIEUmqGzGTj0fqgHflow86n3c+7fP9h7JmdSig83nnh3VoSY8K3QhiRJLqRhj5eGQ3cJ189IFeb+CXN//BbzejJv5BLLqGGxW6EcSIJNWNMPIJTK+hI7hy2myvGX1WhwKunDZbxMkLBEmA8MknOL2GjqDn4KFix6tAkKQII59g+Eph0OX8C2I9NIFAEAOEkY9H0oH+8tEFkcJAA/jQjUADJKludJIkSbEeRKSxWCzk5OTQ0NCQsDlvHCkMfCF88AJBYqHWromF1wRAbQoDmy3JytQLBAJh5OOSzYBRPiJSGGgKD90INESS6kYY+XhEAlrlIyKFgabw0I1AQySpbsTCq4ZRW+xDpDAQCAS+EEZeo/iNlGnnvoDqSGHgz2UjUhgIBMmJcNdokEDFPiq3fuN2XqQwEAgEvhBGXmOoiZRZ9cUSbN9a4bzT50QKA41wHrANN90INEKS6ka4azSGmkiZ+oYqDui30yXdfRerSGGgAdIB4RXTJkmqGzGT1xhqImCyjnUka2ZH2Ostc6QwOO/CkXQ5/wJh4KPNXuB2FHUjiDFJqhth5DWGmgiY9NZscj8qhCNRGJAgOI4ALyF0o0WSVDfCyGsMNcU+MnJyozMYgUAQ9wgjrzHURMoMvuLXqvqy2azs2/4t332xln3bvxVpDQSCJEQsvGoQR6SMZ5x8VocCRk38A13b9Q3Yh8hIKRAIQBh5zeI3UuYAMBN7YWIFfGWkdMTZRyqkUu0O3YTGhF/dCGJIkupGpBpOMGw2Ky/+8baAu19vX/JSWA2w+OUgEEQXkWo4kWkE1shHD2KRkTLQDt2d69eF7bU0jx/dCGJMkupGGPl4ZCcwSj56EGpGylAXaUUuew/86EYQY5JUN8Inn2CEkpGyLa6WYH45iDqzAkH0ETP5BENNnL1rRsq2ulpELnuBQNsII59gBJORMhyuFpHLXiDQNsLIxyPtgM7yUQG1GSnDsUgb7C+HhCeAbgQxJEl1I3zy8UgfYL//JmoyUobD1eL45aAUl+8gqXLZq9CNIEYkqW4iOpOfP38+I0aMICMjg9zcXMU2lZWVTJgwgczMTAoKCrjnnntobW11a7N161ZGjhxJeno6nTt35uGHHyYJwvv9oiYaJlBGynC5WkQue4FAu0R0Jt/a2sq1117L8OHDeemll7zkVquVyy+/nI4dO/L5559z5MgRJk6ciCRJPP3004A94H/s2LGMGjWKiooKfvzxR26++WYyMzO57777Ijl8zbL39S10uKMr/7vwBQ7n2vOmhrLxKJxlA0Uue5mtwKXAf7HPHAXaIUl1E1EjP2/ePACWLl2qKC8rK2PHjh3s27ePoqIiAB5//HFuvvlm5s+fT3Z2Nq+//jrNzc0sXboUo9FIcXExP/74I6WlpUybNg2dThfJt6A5dq5fx1fL3uamY0+hl06rL5SUBeF2tTh+OSQ1J7GnnTgZ64FEBqvNyqbaTRw6foiOGR0Z0GkAhni5kSe4bnwR04XXL7/8kuLiYqeBBxg3bhwtLS1s3LjR2WbkyJEYjUa3NgcPHmTPnj2K/ba0tGCxWNweiUAkNh4JV4tALeV7yyn5Zwm3rryVGZ/N4NaVt1LyzxLK95bHemgCP8R04bW6uhqTyT1bUF5eHqmpqVRXVzvbdO/e3a2N45rq6mp69Ojh1e/ChQudvyISCUc0TAY5PtuEsvFIuFrii1jMpsv3ljN1zVSv87Unapm6ZipPXPIEY7qNiegYBKER9Ex+7ty56HQ6v4+vv/5adX9K7hZJktzOe7ZxLLr6ctXMmjWLhoYG52Pfvn2qx6NlIrnxSJQNjA/K95Yz7l/j3GbT4/41LqKzaavNypwv5vhtM2/dPKzJkroizgh6Jj958mRuuOEGv208Z96+MJvNrF+/3u1cXV0dJ0+edM7WzWazc1bvoLa2FsDrV4ADo9Ho5t5JFBxRLnXtD/L2JbOoa3/QbztBDOgFrJaPYaZ8bznT1kxDwj2yrPZ4LdPWTKP0ktKIzKZf+PYFLCf9uzzrW+v5uuZrhhYODfvrh40I6kbLBG3kCwoKKCjwv/lFLcOHD2f+/PlUVVVRWFgI2BdjjUYjAwcOdLaZPXs2ra2tpKamOtsUFRWpvpkkCq7RMPs7bVVsk1Qbj7RIFnBJ+Lu12qws2rDIy8ADSEjo0LF4w2JGdRnl5bppi3unfG85z33znKq2pV+XckXPK7j+7OtJTUlVdU1UiZButE5EF14rKyvZsmULlZWVWK1WtmzZwpYtW2hqagKgpKSE3r17c9NNN7F582Y++eQTpk+fzqRJk5z5kW+88UaMRiM333wz27Zt4/3332fBggVJGVnjiIZpf7wDv/h2Iu2Pd/Bqk1Qbj7TIAWCWfAwjm2o3UXO8xqdcQqL6eDWbaje5nW+Le8dxY1HLjqM7eLTiUQa/MZjSr0tVXxc1IqQbrRNRI//QQw/Rv39/5syZQ1NTE/3796d///5On73BYODjjz8mLS2NCy+8kOuuu46rrrqKxx57zNlHTk4Oq1atYv/+/QwaNIi7776badOmMW3atEgOXbP0GjqCS6+ZxtDvryWjJdd5XkTDaIQaYJF8lLHaJL786QgfbDnAlz8dwWoLfiPfoeOHgm7ncO943hwc7p1Ahj7QjcUXNsnGy9tf1p6hV9BNMhDR6JqlS5f6jJF30LVrVz766CO/bfr06cOnn34axpHFN1372Gu8jp00mbqigyIaRsOs2FbFvA93UNXQ7DxXmJPGnAm9GV9cqLqfjhkdg2rXFveOA7U3Fl8s27GMyf0ma9N1k0SIBGVxjLlnLxENo2FWbKvirtc2uRl4gOqGZu56bRMrtlWp7mtApwGYMkzoUHZR6tBhzjAzoNMAIHT3jitqbyy+sEk23v7x7Tb1IWg7wsgLBBHAapOY9+EOhXk0znPzPtyh2nVj0BuYOWQmgJehdzyfMWSGc1YeinvHPm4rFdUVLP95OVabldzUXFX9+GKvZW+brhe0HZGFMh7pANwmHwXaQtbNluP1XjN4VySgqqGZDbuPMrynOkWO6TaG0ktKWbRhkdss3ZRhYsaQGW7hk2pn4Ueaj2C1WTHoDZTvLffqOyMlQ1U/vtjfqC7tY1Q2eCXp90YnJUE6R7VVzQWCcPHBlgPc+9aWgO2euqEfv+rXOai+1RhEq83KuH+No/Z4raJf3hVThonLelzG0u1LA7YNloyUDNb9Zp1fg610czFlmJg5ZKbYResHtXZNuGvikRPAdvmYwEhWK8fWb6Dho485tn4DkjUOdlTKujG3S1PVvFOWunauGPQGBpsHc9mZlzHYPFjRgPpz73hSe7yWl7e/HHYDD3D81HEqqit8ytsaARQUSfK98UQY+XjkO6BYPiYolrIydo0eQ+XEiRycPp3KiRPZNXoMlrKyWA/NP7JuBp3IpzAnzad51WGPshnSIz9iQ3G4dzpldPLbLhLG3ZWKGmUjHygCCGDxhsXhS5eQBN8bJYSRF2gOS1kZB+6dwimPdBanamo4cO8U7Rt6wKDXMWdCbwAvQ+94PmdCbwz6yG7oG9NtDCuvXskDgx9oc1+p+hBDIX3cQ8IRASQIjDDyAk0hWa3ULFgISktF8rmaBQvjwnUzvriQ5343AHOOu0vGnJPGc78bEFScvFpco2Mqqiuci6od0tq+2vjM6Ge464K7gl6MbWscflvj9ZMdEV0j0BTHv97oNYN3Q5I4VV3N8a83kjl0SPQGFiLjiwsZ29vMht1HqW1splOW3UUTiRm8vwXMtsS869BhyjAx2DyYYUXDuKPvHby49UVe3f5qwMRlAO/tfI87LrjDy9gHu8FLEBrJZeS3AO1dnucBPYBmYIdC+wHy8QfgmIesO5APHAI8MxlnYc90ZwW+Uei3D/aK8T8BDR6yzoAJqAN2e8jSgfOw/95PAb73kJ8nt9kLHPGQmeS+G4GdHrJ2nC6HthXvyjm95Pd0AO8t4R2AbtgXszx9nTqgv/z/7/Be8OqBXQc1OPOJ2L5sIa25N1Z9IydT94GUQlrL2Xhyqlqe3e2U35MrXYCOwFFgj4csEzhH/r+SF6A3kIb9s/fM2FwoPyzALg+ZETgfRd0Y0DH87A7QE3sh6S0e1xYAXYHjeOtUD/ST/78D+9+qK2cCufBpxae8uPxF8uV/AI0ZjRzodIAHyh/gGdMzXFh9IXXNdU5/9/ddv0fSS3Sr7kZGs/vsvKpDFfVZ9eQ15lF4pJBpg6Zh2GI30usb1vNs5bNgg/Mqz8OTnWfs5FTKKc6oPYOs41n218r5nvMLzociwAzUw4D9A9zG1NKuhZ87/wzAOXvPwSAZyE/LZ8D+AXAQOBfIACoBz4qVnYAzgCbgRw9ZCnABvr83ZwHZQJX8cEXrNkINUhLQ0NAgAVIDDZKEdPrxW7nBTsn9vOPhYJiC7FVZtkRBVuJ4YR/91sryCQqyx2XZOwqy/i5jSlWQb5NltynIZsqy1Qqyzi79dlaQr5ZlMxVkt8mybQqyVJd++yvI35Flj3vLLJmfSDvOOVf6oedwxc+w6ZMK+7UlCv0ukft9VUE2zGVMSrrZKct+qyCbI8tWKMh6uvRboCBfJ8umKsjulmUbFWRZLv32VpB/IEmnrKek//vN/3nJVg5aKRUvLZbGlI5RfK/9X+wvFS8tljacs8FL9tAtD0nFS4ulx+943Eu2pfcWqXhpsdT/xf6K/Y4uHS0VLy2WVg5a6S1fIL+XD7yv21m0UypeWiwVLy2WGtMava/dKF97t8LnMFWWrVOQFbh8hj0V5Ctk2RwFmYZtRMNS2a41NEj+SKo4+S/e2MvQXl1O/1TW+l3a10weYDPeC1oJMJOXrFb23XkXrfV7Odmu0nsmrwNDfge6fPkculSD9mbyAN8CpzzkZ2P/FbkfqPWQtXEmX9Fcwf1v30+nevcomsaMRvZ32k/KqRR67e/FnOFzaDrZxMvbXuZo81HnTH5g/UDuOusuhhYOxWqz8v3R7zmYf5DcM3IZoB+AYf9pN8v2w9uZsXkGe8170dl0nFt5Lp4ozeTnDJ/jNZPHPmlnfdV6Xt72MlXWKudM/sLaC7ml9y3u+enbOpMHewhli4c8Tmfylr0WcroHjpNPKiPfZco7dO6UH3RyKM3xHfBb4HVOG/0oI1mtdv/5oUOkdOxIxqCB6Azh2aHoiK6xv5DLn6ecWrrzU0+SXVISltcKOzHQzfKflzPjsxkB2y2+aDGXnXlZm3aXqn0tBw5//oqrV/h9jajseNXA9yacqN0MlVw+eU4nh4pUdENUOIF9Jh+jTR2WsjJqFix0WyBNMZsxzZ4VFuObXVICTz3p/RomU9heI2LEQDfBLmA6NlNF8rVccc2p44u2jEk1Mf7exIqkM/ISdi/CvA93MLa3OeJxyvFAMLNy5yzb4wegI4adMM2ys0tKyBo9OmK/FhIJR4ZKXykMdOjolNEJm2Rj+c/L2zRTdryW2jzzj170qEhNEGOSzshDaMmhEpVgZuUBY9h1OmoWLCRr9OiwGGOdwRAXYZKxxpHCYNqaaejQuRl6x/PmU83cXna783youWEMegPXnH0Nz2x5RlX7Q80ixj3WJPVmqNpG31kCk4Fgd5YGE8MuiC6+UhjkGHMAaGh1X71rS26YrlldVbfdZ/FccRREm6ScyTsIJTmUJugBvCMfQySUWfmpQ+pmZWrbJSRh0E2ojOk2hlFdRjkXMPPT8nnw8wcV20pydahF6xfRvl17jjYfVe3GCcYv3yW7S1DvIaLEUDexJCmNvA771vJIJoeKKHnAtW3rIpSdpSkd1X251baLB6w2KbjdqmHQTVtwXcCsqK6g9oRnvOZpJCRqTtQwadUk5zk1bpwBnQbQMa1jQFeMXqfn+rOvD/IdRJAY6yZWJJ27JprJoSJGDVBKmwoShzIrzxg0kBSz2RnK6IVOR4rZTMaggaEPTEOs2FbFLxb/j9+8+BX3vrWF37z4Fb9Y/D//ZfvCoJtwEUrOFzVuHIPewOxhswP2NbH3RG3Vd9WQbqJJ0hn5SCaHihoHgPtwbiAKhVBm5TqDAdPsWfITD0MvPzfNnqV60VXL+eJDrs8aBt2Ei1DCHSWVKX7HdBvDE5c84TNZ2fju47l3wL1Bv35E0ZBuoklSuWv+MXEwoy7oFr8z+DDimJWfqqlR9svrdKSYTF6z8nDFsEc61r4tBKrPGi8huIFCK30huaT4HWwe7MxsWVFTARIMKRzCINMg5xrAhuoNvPDNC2w7so0Wq3076Yo9K9hcu1lUd9IASWXkh5wZmex/8YhjVn7g3in2WbjCzlJfs/K2xrBHK9Y+VDbsPhr2+qyxwF9opRoOHT9E+d5y5n45l4aW09E5f9/6d3JSc5g7Yi5juo3h2MljbKrd5NW/w/VTekmpMPQxJOncNYLTZJeU0PmpJ0kxmdzOp5hMAVMHOGLYc664nMyhQ4Jy0Wg9X7za0Np4CMFVWx1KicrGSqaumepm4B00tDYwdc1UyvaURbe6kyBokmomnzDkABPkYxuJ9s7SeMgXrza0VrFdGHUTLjxDKzukd+DBzx/0v0M2vRPv/vBuwL7/+uVfqW+t9yn3dP3EFA3qJhoIIx+P9AT+E77uormzNB5i7Yf0sNdnrW5oVnRw+A3BDbNuwoVnbhh/O2QB+pv6s2LPioD9+jPwrmiiupNGdRNphLsmHjmJPX2pZ0rgOKCtsfbRiMhpU33WONGNLzeOKcPEzeffrMrAB4MmqjvFiW7CjZjJxyNbgYHARk7ns44TQo3qgehG5Djqs877cIfbIqw5J81/quo40o2nG6djRkf6FvTlsvcvC6qfPGMe9S31Pl0/pgwTAzpp4MOII92EE2HkBVEl1KieWETkRLM+a6zwdON8dfAr1RkmAUzpJh4Y8gDT10736fpRk2pYEDmEkRdEnWBj7aOd/dIVg16n6TDJtuJarKPSUskrO14J6vqZQ+1x8KW6UsUi4jOGzBDhkzFGGHlBTAgmqiceInJiTSiVlcr3lnsZZrVkpmTyyC8ecRpwJddPRKo7CYJGGHlBRPFXkERtVE88ROTEEiVjHSjRWPnecqatmRb0Bimw++DLryknNSXV6+bSt6Av3xz+hkPHD7GpdpMw9BogokZ+/vz5fPzxx2zZsoXU1FTq6+u92ugUkl0999xz3Hnnnc7nW7duZfLkyWzYsIH8/HzuuOMO/vKXvyhemxT0xV7cNzPWA/FPuBZK4yr7ZZR148tY1xyvYeqaqTw+8nFKup/+rK02K1/XfM3cdXNDMvAAfx72Z1JTUhVvLnqdHptkcz4PtThJRIiT7024iWgIZWtrK9deey133XWX33Yvv/wyVVVVzsfEiROdMovFwtixYykqKqKiooKnn36axx57jNLS0kgOXdsYsFeX1/AEKZiCJIHCIuMq+2UUdWO1WX3uNnUwfe10VvxsD4cs31vOuH+N4/ay272KiARDXlqe8+bi6epxNfDgndXSkQdn+c/LqaiuiO5O2Dj43kSCiM7k582bB8DSpUv9tsvNzcVsNivKXn/9dZqbm1m6dClGo5Hi4mJ+/PFHSktLmTZtWnLO5ncCk4ElQK8Yj0WBYBZKGz/5JOBsP9SInGBq14aNKOpmU+2mgP50CYn7P7uf5XuWs2bfmpBn767UHKvhyU1PqurLUZxk8YbF2CQbj1Y8GpRbqa24upPOqD6DPnP7oHtGp8nvTaTQxGaoyZMnU1BQwODBg3n++eex2U7PBr788ktGjhyJ0Wh0nhs3bhwHDx5kz549iv21tLRgsVjcHglFI1AmHzWI2oXSw8+/oHq2H2yeHUtZGbtGj6Fy4kQOTp9O5cSJ7Bo9xqukYdiJom6C2UW6et/qsBh4gLqWuqAWax2pDe5be5/XdW0pQxgIxy+XW1feyozPZvDIqkfQrdKx/sf1YX8tLRNzI//Xv/6Vd999l/Lycm644Qbuu+8+FixY4JRXV1dj8vhiO55X+zAkCxcuJCcnx/no0kVDJciSALULoEdfeSWoRGXZJSWc9Uk5XZcto+ixx+i6bBlnfVKuaOCDqV0br0R7F6kOHeYMM3nGvLD1GakkZr7cSQClX5dG5KaiVYI28nPnzkWn0/l9fP3116r7+/Of/8zw4cPp168f9913Hw8//DB/+9vf3Np4umQk2Qj4ctXMmjWLhoYG52PfvsQuJmyzSRz4oY4fK6o58EMdNlt4ZmyhonYB1Nbgxy/soyh4oOyX8ZDlMlw48sVHA9eNTabM8L6maxIzaLvfXs1aRTJlxgzaJz958mRuuOEGv226d+8e6ngYNmwYFouFmpoaTCYTZrPZa8ZeW2uvW+k5w3dgNBrd3DuJzE+ba/ns7Z0cq29xnsvMNXLR9b3o2T/49LLhIGDqAvD2rfsg2LDIZIqpd+SLn7pmatj61KNnbLexfFX1ldvirOvGJqvNiinDFFJ8vT8c+euDDQf1JNBahaYyY0aBoI18QUEBBQUFkRgLAJs3byYtLY3c3FwAhg8fzuzZs2ltbSU11V4vsqysjKKiojbdTOKaLsAS2F1/iBXvbPMSH6tvYcUL2xh/R3FMDL3bQqkvVBh4CD4sMuYx9bJuCLOH0NdmpzHdxvDbc3/L69+/HpbXsWHj+nOvZ/HFi31ubDLoDTww+AHuW3tfWF7TQaWlkme/ebbNxUd8rVVU51cz/3fzqc6v9tsu0YhodE1lZSVHjx6lsrISq9XKli1bADjrrLNo3749H374IdXV1QwfPpz09HRWr17Ngw8+yB/+8AfnTPzGG29k3rx53HzzzcyePZudO3eyYMECHnrooeSMrAHoCLa7JNbO/tFvs8/f2UmPvh3RxyDXSnZJCdITpRycdh/YbIEv8MRPojJ/xDymviPwx/B2GWh2O7rb6LAZebAbP8+cNg4ccfble9T5tHXo+H3v37Nizwr/+eszOvHuj+/6LD7iiNAZ1WVUwM1VvtYq6rLreGvMWwHbJRoRXXh96KGH6N+/P3PmzKGpqYn+/fvTv39/p8++Xbt2PPvsswwfPpwLLriAp556iocffpjHH3/c2UdOTg6rVq1i//79DBo0iLvvvptp06Yxbdq0SA5d2xyF+sePc6rav0+xqa6Fqp310RmTAil5+SEbeAiuKDjY/fGSzYY+x09ViEjH1B8FXpOPYaBsTxlT10z1G5USbt/8Zwc+4+lNT7O+ar2b37p8bzkj3x7J7WW389+9/1XVl4TEKzte4bIe9syWOo/kzY7n15x9DbUnav324+q394fj8/B8reymbK5YdwU5TTmYM8zayIwZBXSSpPJ3cxxjsVjIycmhoaGB7OzsWA+n7WwCBsLbv97A4Y5NfpuOva03Zw9W3oMQaRo++piD06cHbKfPyXFbhA1lV6zS7lov5JtHoNKGbULWTTjS2ZbtKeP+T+/32mDkwJHGd8XVK1i9b7XfNAWjuoxi9b7VQY/BqDdy9dlX09TaxH9+Dq3ihmOc9w++n8UbFrsZc1O6iZlDZ9JqbWXGZzMC9rX4osVcdmbgVMiO6Bo4HcFz3p7zeGfuO1w/93om3TxJG7tw24BauyZy1yQ4mdmxW4BW6xLp/OQT6PSGkDct+UpD7DUeH1kutUj53vKAPm/H7PaN79/gxnNvpPQS70yQecY8/jzsz5R0L6F8bzn3r72fU9Ip1eNosbXwxvdvhPw+XMe5q24XJ20eFTvkybZa14nado6iKEoJ2KYNmsbQbkNV9ZMICCMfx6RnpfqVt88zUtgrN+h+w7VTVFWUjV6PraGB7PHjg+7fMVafIZOOl8jJofOTT5A5RH3B8VjiCAFUy6MVj7Js+zJmDpnJyqtX+lww3VK7JSgDH26e+/Y5r3M1x2uYtmYav+/9e6+8N66EUnzEMzNmt8JuAAwtTB4DD8LIxzX9S7qy71Pfzt9fXNcr6EXXcFZfUhVlY7NxYOo00OtDmmEHDJnEHo+v0xuiauCtNiubqkNLu6smXYEnviJQHDHnVU1VLNuxLKg+o4WEpGpsoRQfcVtArg9hcAmAMPLxSCYwDLoMzGf8ecVecfLt84z84rrg4+T9Vl+6516OT5xI1i9/GdTMXm2UTahFP2IeMulJJtT3r2fKF1PYmHt6I1cwsd6hhPY5IlAWbVhEVmoWR04codJSybs/vut3QTMe0Ov0/O3iv7Xdhy5/b5ItC6Uw8vHIOcCX9v/2pBM9+nakamc9xywtZGbbXTTBzuDV7BStW7aMumXLMJhM5F13LandujvdOYBPF0/AKJs2bFCKecikB+Vp5Uy713sBNJhY71BD+yQkao7XcHvZ7SFdr1Vsko28tDCkUnD53iQTSWfkbTbJzSCaeuZQ81NDmwxkrNHrdXQ+p21fAjVuDwfWmhoOP73k9OvLG9dsLvUCXF08kZxtt6UweLjxt53eV6y30iYnRwigr7jyZCRZNi5FgqQy8j9/c4gtH291c2147q6PdUoAVYQxTM9BW9wZNoViMK4FtiM52w41DXEk2FS7ifwd+ZTPLee6udfxXffv3OSe2+n9bXIKd7qCeCcsG5ci8L2JB2KehTKalL+8w83Ag/fkz5ES4KfN8e3HDJawuzNckoGl9+8X0aIfwaYhjhRqZ5uOHC1KWRIdbp1vD30biSHGHY7Ml8mycSkSJNVMPhhimRIgFqgKdwwW2dd+YvOWiM+2gykM7o9QCmI7UDvbzE/L589f/NmnWwfQbCRMJBjVZRRr9q0BcPtMXDNfijqxoSOMvA8cKQHa6uuOF/y6PdrIqUOHyLnicnjqSe/wzDBuUFJbGNwXbc2AOKDTAPLT8n2PT4711ul0AUMkfcWLJxKeG7WUPntH5ktB6Agj74djlpbAjRKI7JISRUPcVhyuoHDNtiOBr4LYwUTFGPQGbim+BfCdo2XGkBkcOXEkLGPuW9CXG8+7kbX717J89/Kw9Bkt/nDBH7i7793OGbrnxqVgf0W1lbb8gtM6wsj7IZYpARx4RgMV9spF31tnryV6Rvhfz2GIj776GrWL1O+6VEQhssVztu0o4h1Lox9KVIwvho4eyhdrvsCy3wIum0tdZ6UV1RVhGfc3h78h9+dcBhcOjjsjP6xwmNdn6SvzpS+CNsy9UfzehCOHvZYRRt4HoaYECCexKgiiMxhIaWvNABW+9nDurm0LYS0ykQYXjryQj2wf+TRAakIk/W3xd2XtgbWsPbA2YDstkWPMwSbZsNqsIc+WQzLMacBZ3v209RdcuFCc0IVhTTCpomuCIZSUAOHkp821rHhhm1c00LH6Fr742y4aL2uG3e7XhLMMYFujbQJFtmipDmswUTEB2Q38Dgx77bPSy868jMHmwW7GzFHRCZTdOjp0TOw9UfX4442GlgZuL7udcf8aF1Kt1UCRST77lHXj+N4E+gUH0SsT+NPmWpbN+oJ/P7GZVS/t4N9PbGbZrC/CEuWXVEZ+zC29ycx1d8F4RvW1zzPGrKKSA5tN4rO3d/qUG1tSyPpvGrYjp/84f9pcyyuz17n9kbwye13IfySOaBufYY+e6HQY8vMpfPRRnwW2HWitDmtYMyDWAa/LRz84siR2ynD/OzNlmCi9pJRpg6bxi6JfqBpXvBLQKCvQJsPsoZtgfsFFEseE7nhDq9v54w2tYQnnTip3zZl9O3LBhWdqfsdr1c56rxm8EocrG+k0KNv5R+JJW8oABh1tI0lYjx6lnckUMMJFa3VYA7lPQsmAqIZAi40jOo/g84Ofh/U1o8m4buO45uxruG/tfVhaLV7yYNc7ILyutbD+ggsRm01i9avf+22z5rXv2xTOnVQzeTidAuDswWY6n5NHSore7XmsDTyoj+o5caw14Kwf7DH/obhufG0y8oeanbNaSyoWyH0CkYvVdiw2Krl1rj/7evS6+P2KflX1FYCigXcQ7Gw5nIY53DnsQ+HgD3W0HPef/rn52CkO/hDgp6EfkmIm7yh+ZbH4/mPTEjZ9Cydaj/mUHz8JFiyc0jXz4+ZKDtf6rzV3ouYYP26upKhXCDH/w4bR6f33qHvrbQ6VlgZsfjwjHV2Az/l4RjpNKlwxavoKF0PyhvDXgX/liY1PeFUumjJwCkPyhqj7+2lyOYZh6Nd3vZ7Xvn+t7R3FgKMnjvLpT59iPRFY13tr9nJOxjkB26Vb01X1l25N99aXh27OSjuLDnTg0IlDvmvPpnfirLSzImY7fvim0u933bVddmd3c+0YU6DifklR/m///v106dIl1sMQCASCsLNv3z7OOMN3PHVSGHmbzcbBgwfJyspCp3YhMcGxWCx06dKFffv2JUbdW40gPtfIID5XbyRJorGxkaKiIvR63269pHDX6PV6v3e6ZCY7O1t8aSKA+Fwjg/hc3cnJyQnYJn5XdQQCgUAQEGHkBQKBIIERRj5JMRqNzJkzB6Mx9vl5EgnxuUYG8bmGTlIsvAoEAkGyImbyAoFAkMAIIy8QCAQJjDDyAoFAkMAIIy8QCAQJjDDyAoFAkMAII5/gzJ8/nxEjRpCRkUFubq5im8rKSiZMmEBmZiYFBQXcc889tLa657beunUrI0eOJD09nc6dO/Pwww8HTIyUbDz77LP06NGDtLQ0Bg4cyGeffRbrIWmaTz/9lAkTJlBUVIROp+Pf//63m1ySJObOnUtRURHp6elccsklbN++3a1NS0sLf/rTnygoKCAzM5Mrr7yS/fv3R/FdaB9h5BOc1tZWrr32Wu666y5FudVq5fLLL+fYsWN8/vnnvPXWW/zrX//ivvvuc7axWCyMHTuWoqIiKioqePrpp3nssccoVZGVMll4++23mTJlCg8++CCbN2/moosu4tJLL6WysjLWQ9Msx44do2/fvixZskRR/uijj1JaWsqSJUuoqKjAbDYzduxYGhsbnW2mTJnC+++/z1tvvcXnn39OU1MTV1xxBdYoFZyJCyRBUvDyyy9LOTk5XueXL18u6fV66cCBA85zb775pmQ0GqWGhgZJkiTp2WeflXJycqTm5mZnm4ULF0pFRUWSzWaL+NjjgSFDhkh33nmn27lzzz1XmjlzZoxGFF8A0vvvv+98brPZJLPZLC1atMh5rrm5WcrJyZGef/55SZIkqb6+XmrXrp301ltvOdscOHBA0uv10ooVK6I2dq0jZvJJzpdffklxcTFFRUXOc+PGjaOlpYWNGzc624wcOdJtt+G4ceM4ePAge/bsifaQNUdraysbN26kxKPcYUlJCevWrYvRqOKb3bt3U11d7faZGo1GRo4c6fxMN27cyMmTJ93aFBUVUVxcLD53F4SRT3Kqq6sxeVR+ysvLIzU1lWq5RJ9SG8fzan9l/JKEw4cPY7VaFT8j8fmEhuNz8/eZVldXk5qaSl5ens82AmHk45K5c+ei0+n8Pr7++mvV/Snl2Jckye28ZxtJXnQV+flPo/QZic+nbYTymYrP3Z2kyCefaEyePJkbbrjBb5vu3bur6stsNrN+/Xq3c3V1dZw8edI5izKbzV4zo9pae4k8z5lWMlJQUIDBYFD8jMTnExpmsxmwz9YLCwud510/U7PZTGtrK3V1dW6z+draWkaMGBHdAWsYMZOPQwoKCjj33HP9PtLS0lT1NXz4cLZt20ZVVZXzXFlZGUajkYEDBzrbfPrpp25hlWVlZRQVFam+mSQyqampDBw4kFWrVrmdX7VqlTA2IdKjRw/MZrPbZ9ra2sratWudn+nAgQNp166dW5uqqiq2bdsmPndXYrrsK4g4e/fulTZv3izNmzdPat++vbR582Zp8+bNUmNjoyRJknTq1CmpuLhYGj16tLRp0yapvLxcOuOMM6TJkyc7+6ivr5dMJpP0m9/8Rtq6dav03nvvSdnZ2dJjjz0Wq7elOd566y2pXbt20ksvvSTt2LFDmjJlipSZmSnt2bMn1kPTLI2Njc6/R0AqLS2VNm/eLO3du1eSJElatGiRlJOTI7333nvS1q1bpd/85jdSYWGhZLFYnH3ceeed0hlnnCGVl5dLmzZtkn75y19Kffv2lU6dOhWrt6U5hJFPcCZOnCgBXo/Vq1c72+zdu1e6/PLLpfT0dCk/P1+aPHmyW7ikJEnSt99+K1100UWS0WiUzGazNHfuXBE+6cEzzzwjdevWTUpNTZUGDBggrV27NtZD0jSrV69W/NucOHGiJEn2MMo5c+ZIZrNZMhqN0sUXXyxt3brVrY8TJ05IkydPlvLz86X09HTpiiuukCorK2PwbrSLyCcvEAgECYzwyQsEAkECI4y8QCAQJDDCyAsEAkECI4y8QCAQJDDCyAsEAkECI4y8QCAQJDDCyAsEAkECI4y8QCAQJDDCyAsEAkECI4y8QCAQJDDCyAsEAkEC8/8DNWveWNrduIkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "data = ptm_data[ptm_data.residue==\"SEP\"][[\"chi2\", \"chi1\"]]\n",
    "fit_data = np.stack((np.cos(np.radians(data)), np.sin(np.radians(data))), axis=1)\n",
    "clustering = DBSCAN(eps=0.4, min_samples=3).fit(fit_data.transpose(0, -1, 1).reshape(-1, 4))\n",
    "\n",
    "plt.figure(figsize=(4, 4))\n",
    "for l in np.unique(clustering.labels_):\n",
    "    plt.scatter(data[clustering.labels_ == l].chi1, data[clustering.labels_ == l].chi2)\n",
    "\n",
    "plt.xlim(-180, 180)\n",
    "plt.ylim(-180, 180)\n",
    "\n",
    "# Highlight grid lines at x = [-120, 0, 120] and y = [-120, 0, 120] with magenta\n",
    "for val in [120, 0, -120]:\n",
    "    plt.axvline(x=val, color='magenta', linestyle='--', linewidth=1)\n",
    "    plt.axhline(y=val, color='magenta', linestyle='--', linewidth=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ptm = \"SEP\"\n",
    "#calculate terminal rotamer\n",
    "data = ptm_data[ptm_data.residue==ptm]\n",
    "\n",
    "r3_stats = np.zeros((3, 3))\n",
    "r3_class = data[free_rotamer[ptm]].apply(rotameric_class)\n",
    "for l in range(3):\n",
    "    ang_rad = data[r3_class==l][free_rotamer[ptm]].apply(np.radians)\n",
    "    kappa, loc, _ = vonmises.fit(ang_rad, fscale=1)\n",
    "    weights = vonmises(loc=loc, kappa=kappa).pdf(ang_rad)\n",
    "    wm = weighted_circmean(ang_rad, weights)\n",
    "    r3_stats[l, 0] = len(ang_rad)/len(data)\n",
    "    r3_stats[l, 1] = np.degrees(wm)\n",
    "    r3_stats[l, 2] = np.degrees(weighted_circstd(ang_rad, weights, wm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = []\n",
    "\n",
    "for l in np.unique(clustering.labels_):\n",
    "    angles = data[clustering.labels_ == l]\n",
    "    if l == -1:\n",
    "        # for outliers just append individual point    \n",
    "        for _, pt in angles.iterrows():\n",
    "            line = np.concatenate(([1/len(data)], pt[[f\"chi{n+1}\" for n in range(nchis[ptm])]].values, np.zeros(nchis[ptm])))\n",
    "            results.append(line)\n",
    "    else:\n",
    "        line = np.zeros(1+nchis[ptm]*2)\n",
    "        for n in range(3):\n",
    "            ang_rad = angles[f\"chi{n+1}\"].apply(np.radians)\n",
    "            kappa, loc, _ = vonmises.fit(ang_rad, fscale=1)\n",
    "            weights = vonmises(loc=loc, kappa=kappa).pdf(ang_rad)\n",
    "            wm = weighted_circmean(ang_rad, weights)\n",
    "            line[1+n] = np.degrees(wm)\n",
    "            line[1+n+nchis[ptm]] = np.degrees(weighted_circstd(ang_rad, weights, wm))\n",
    "        \n",
    "        for r3 in r3_stats:\n",
    "            #factor in free rotamer\n",
    "            line[0] = len(angles) / len(data) * r3[0]\n",
    "            line[-1-nchis[ptm]] = r3[1]\n",
    "            line[-1] = r3[2]  \n",
    "            results.append(line.copy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = format_bb_indp_table(results, ptm)\n",
    "#print(table)\n",
    "with open(\"sep_bi.lib\", \"w\") as f:\n",
    "    f.write(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "backbone dependent rotamer library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "header = \"#AA   PHI   PSI    Prob   \"\n",
    "line_formatter = \"{: <3}  {:4d}  {:4d}  {:8.4f}\" #\"{:8.1f}\"\"{:8.2f}\"\n",
    "\n",
    "def get_closest_angle(input_angle, degsep=10):\n",
    "    \"\"\"\n",
    "    Find closest angle in [-180, +170] with n degree separations\n",
    "    \"\"\"\n",
    "    if input_angle > 180 - degsep/2.:\n",
    "        # the closest would be 180, namely -180\n",
    "        return -180\n",
    "    return round(input_angle/degsep)*degsep \n",
    "\n",
    "def format_bbdp_table(lines, residue):\n",
    "    n_ang = (len(lines[0]) - 3) // 2\n",
    "    formatter = line_formatter + \"{:8.1f}\" * n_ang + \"{:8.2f}\" * n_ang\n",
    "    table = header + \" \".join([f\"chi{i+1}Val\" for i in range(n_ang)]) + \" \" + \" \".join([f\"chi{i+1}Sig\" for i in range(n_ang)]) + \"\\n\"\n",
    "    #table = \"\"\n",
    "    for l in lines:\n",
    "        table += formatter.format(residue, int(l[0]), int(l[1]), *l[2:]) + \"\\n\"\n",
    "    return table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_bb_dependent(ptm, data, bb_sep=30, last_free_rot=False):\n",
    "    results = []\n",
    "    nchi = nchis[ptm]\n",
    "\n",
    "    #scbins = product(range(3), repeat=nchi)\n",
    "    chi_classes = []\n",
    "    for n in range(nchis[ptm]):\n",
    "        chi_classes.append(data[f\"chi{n+1}\"].apply(lambda x: rotameric_class(x, rot_bin_lib.get((ptm, n), 3))).values)\n",
    "    effect_chi_classes = np.array(chi_classes)\n",
    "    if last_free_rot:\n",
    "        nchi -= 1\n",
    "        effect_chi_classes = effect_chi_classes[:-1]\n",
    "\n",
    "    #bbins = product(range(360//bb_sep), repeat=2)\n",
    "    bb_classes = []\n",
    "    for ang in [\"phi\", \"psi\"]:\n",
    "        bb_classes.append(data[ang].apply(lambda x: get_closest_angle(x, bb_sep)))\n",
    "    bb_classes = np.array(bb_classes)\n",
    "\n",
    "    for b in np.unique(bb_classes.T, axis=0):\n",
    "        bb_mask = np.ones(len(data), dtype=bool)\n",
    "        for i, l in enumerate(b):\n",
    "            bb_mask &= bb_classes[i] == l\n",
    "        if np.sum(bb_mask) == 0:\n",
    "            continue\n",
    "        if np.sum(bb_mask) == 1:\n",
    "            line = np.concatenate(([int(b[0]), int(b[1]), 1/np.sum(bb_mask)], \n",
    "                                    data[bb_mask][[f\"chi{n+1}\" for n in range(nchis[ptm])]].values[0], np.zeros(nchis[ptm])))\n",
    "            results.append(line)\n",
    "            continue\n",
    "\n",
    "        for c in np.unique(effect_chi_classes.T, axis=0):  \n",
    "            mask = bb_mask.copy()\n",
    "            for i, l in enumerate(c):\n",
    "                mask &= effect_chi_classes[i] == l\n",
    "            if np.sum(mask) == 0:\n",
    "                continue\n",
    "\n",
    "            angles = data[mask]\n",
    "            if np.sum(mask) == 1:\n",
    "                line = np.concatenate(([int(b[0]), int(b[1]), 1/np.sum(bb_mask)], \n",
    "                                       angles[[f\"chi{n+1}\" for n in range(nchis[ptm])]].values[0], np.zeros(nchis[ptm])))\n",
    "            else:\n",
    "                line = np.zeros(2*nchis[ptm] + 3)\n",
    "                prob = np.sum(mask) / np.sum(bb_mask)\n",
    "                line[0] = int(b[0])\n",
    "                line[1] = int(b[1])\n",
    "                for n in range(nchi):\n",
    "                    ang_rad = angles[f\"chi{n+1}\"].apply(np.radians)\n",
    "                    kappa, loc, _ = vonmises.fit(ang_rad, fscale=1)\n",
    "                    weights = vonmises(loc=loc, kappa=kappa).pdf(ang_rad)\n",
    "                    wm = weighted_circmean(ang_rad, weights)\n",
    "                    line[3+n] = np.degrees(wm)\n",
    "                    line[3+n+nchis[ptm]] = np.degrees(weighted_circstd(ang_rad, weights, wm))\n",
    "                if last_free_rot:\n",
    "                    for j in range(3):\n",
    "                        sub_mask = bb_mask & (chi_classes[-1]==j)\n",
    "                        line[2] = prob * np.sum(sub_mask) / np.sum(bb_mask)\n",
    "                        ang_rad = data[sub_mask][f\"chi{nchis[ptm]}\"].apply(np.radians)\n",
    "        \n",
    "                        if len(ang_rad) == 0:\n",
    "                            continue\n",
    "                        elif len(ang_rad) == 1:\n",
    "                            line[-1-nchis[ptm]] = np.degrees(ang_rad.values)\n",
    "                            line[-1] = 0.\n",
    "                            results.append(line.copy())\n",
    "                            continue\n",
    "                        kappa, loc, _ = vonmises.fit(ang_rad, fscale=1)\n",
    "                        weights = vonmises(loc=loc, kappa=kappa).pdf(ang_rad)\n",
    "                        wm = weighted_circmean(ang_rad, weights)\n",
    "                        line[-1-nchis[ptm]] = np.degrees(wm)\n",
    "                        line[-1] = np.degrees(weighted_circstd(ang_rad, weights, wm))\n",
    "                        if line[-1] > 30:\n",
    "                            print(np.degrees(ang_rad))\n",
    "                        results.append(line.copy())\n",
    "                    continue\n",
    "                line[2] = prob\n",
    "            results.append(line)\n",
    "        \n",
    "    return results\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = ptm_data[ptm_data.residue==\"ALY\"]\n",
    "results = calc_bb_dependent(\"ALY\", data, bb_sep=30, last_free_rot=False)\n",
    "table = format_bbdp_table(results, \"ALY\")\n",
    "with open(\"aly_bd.lib\", \"w\") as f:\n",
    "    f.write(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mda",
   "language": "python",
   "name": "mda"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}