HPLC-Code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1ae7b0c4",
   "metadata": {},
   "source": [
    "# Import Moduli"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8cd4ad2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import csv\n",
    "import openpyxl\n",
    "import re\n",
    "import glob\n",
    "import os\n",
    "from os.path import exists as file_exists\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abcece84",
   "metadata": {},
   "source": [
    "# Set Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2d501858",
   "metadata": {},
   "outputs": [],
   "source": [
    "uL_plasma = 5\n",
    "\n",
    "cys_factor = 1.0213\n",
    "hcy_factor = 1.6548\n",
    "\n",
    "x1_values = [12.5,25,50,100,200,400]                      #Standard cysteine concentrations\n",
    "x2_values = [18.75,37.5,75,150,300,600]                   #Standard homocysteine concentrations\n",
    "\n",
    "peak_1_min = 1.6                                          #Set time window to look for peaks\n",
    "peak_1_max = 1.8\n",
    "\n",
    "peak_2_min = 3.2\n",
    "peak_2_max = 3.4\n",
    "\n",
    "peak_3_min = 6.1\n",
    "peak_3_max = 6.2\n",
    "\n",
    "minimum_peak_area = 10                                    #Sensitivity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b0c357b",
   "metadata": {},
   "source": [
    "# Define Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "aec9daa4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_standards(file_name):\n",
    "    with open(file_name, encoding=\"utf-16\", errors='ignore') as f:\n",
    "        contents = f.read()\n",
    "\n",
    "    l = re.split('[^0-9.]', contents)\n",
    "    l1 = [l[i] for i in range(0, len(l), 2)]\n",
    "    l2 = [l[i] for i in range(1, len(l), 2)]\n",
    "\n",
    "    time = []\n",
    "    area = []\n",
    "    for i in range(len(l2)):\n",
    "        A = l1[i]\n",
    "        B = l2[i]\n",
    "        floatA = float(A)\n",
    "        floatB = float(B)\n",
    "        if peak_1_min <= floatA <= peak_1_max or peak_2_min <= floatA <= peak_2_max or peak_3_min <= floatA <= peak_3_max:\n",
    "            if minimum_peak_area <= floatB:\n",
    "                r1 = l1[i]\n",
    "                r2 = l2[i]\n",
    "                newtime = time.append(r1)\n",
    "                newarea = area.append(r2)\n",
    "                if len(area) > 3:\n",
    "                    print('check ' + file_name + ' manually!')\n",
    "                    return\n",
    "            else:\n",
    "                print('check ' + file_name + ' manually!')\n",
    "                return\n",
    "    \n",
    "    cys = area[0]\n",
    "    hcy = area[1]\n",
    "    std = area[2]\n",
    "    adj_cys = float(cys)/float(std)\n",
    "    adj_hcy = float(hcy)/float(std)\n",
    "    \n",
    "    cys_values = y1_values.append(adj_cys)\n",
    "    hcy_values = y2_values.append(adj_hcy)\n",
    "\n",
    "        \n",
    "        \n",
    "def calculate_HPLC_results(file_name):\n",
    "    with open(file_name, encoding=\"utf-16\", errors='ignore') as f:\n",
    "        contents = f.read()\n",
    "\n",
    "    l = re.split('[^0-9.]', contents)    \n",
    "    l1 = [l[i] for i in range(0, len(l), 2)]\n",
    "    l2 = [l[i] for i in range(1, len(l), 2)]\n",
    "\n",
    "    time = []\n",
    "    area = []\n",
    "    for i in range(len(l2)):\n",
    "        A = l1[i]\n",
    "        B = l2[i]\n",
    "        floatA = float(A)\n",
    "        floatB = float(B)\n",
    "        if peak_1_min <= floatA <= peak_1_max or peak_2_min <= floatA <= peak_2_max or peak_3_min <= floatA <= peak_3_max:\n",
    "            if minimum_peak_area <= floatB:\n",
    "                r1 = l1[i]\n",
    "                r2 = l2[i]\n",
    "                newtime = time.append(r1)\n",
    "                newarea = area.append(r2)\n",
    "                if len(area) > 3:\n",
    "                    print('check ' + file_name + ' manually or change bounds!')\n",
    "                    return\n",
    "            else:\n",
    "                print('check ' + file_name + ' manually or change bounds!')\n",
    "                return\n",
    "                \n",
    "    cys = area[0]\n",
    "    hcy = area[1]\n",
    "    std = area[2]\n",
    "    adj_cys = float(cys)/float(std)\n",
    "    adj_hcy = float(hcy)/float(std)\n",
    "\n",
    "    conc_cys = ((adj_cys-int_cys)/slope_cys)*cys_factor\n",
    "    conc_hcy = ((adj_hcy-int_hcy)/slope_hcy)*hcy_factor\n",
    "\n",
    "    final_cys = (10/uL_plasma)*conc_cys\n",
    "    final_hcy = (10/uL_plasma)*conc_hcy\n",
    "\n",
    "    data = pd.DataFrame({'Time':[time[0],time[1],time[2]],'Area':[area[0],area[1],area[2]],'Adj':[adj_cys, adj_hcy, '1'],'Conc':[conc_cys, conc_hcy, None],'Final':[final_cys, final_hcy, file_name]})\n",
    "    index_ = ['Cys', 'Hcy', 'Std']\n",
    "    data.index = index_\n",
    "    \n",
    "    new_file_name = file_name.removesuffix('.CSV')\n",
    "    running_ids = ids.append(new_file_name)\n",
    "    final_cys_values = final_cysteine.append(final_cys)\n",
    "    final_hcy_values = final_homocysteine.append(final_hcy)\n",
    "    \n",
    "    if file_exists('compiled.xlsx'):\n",
    "        with pd.ExcelWriter('compiled.xlsx', engine=\"openpyxl\", mode='a', if_sheet_exists='replace') as writer:  \n",
    "            data.to_excel(writer, sheet_name=file_name, index=True)\n",
    "    else:\n",
    "        data.to_excel('compiled.xlsx', sheet_name=file_name, index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3751892b",
   "metadata": {},
   "source": [
    "# Get Standard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4c3528d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvsElEQVR4nO3dd3hUZdrH8e8dCCV0CD2EhN5RiCCyKogFLIvYVoyKWLLu6urqqoCAiIqirr52WVbBhroqoCgoICooawEU0uglEEINEEoIpNzvH+fgRpyEIWRyJpP7c11zzZw6vxwgN6c8zyOqijHGGHO8MK8DGGOMCU5WIIwxxvhkBcIYY4xPViCMMcb4ZAXCGGOMT5W9DlCaIiMjNSYmxusYxhhTbixbtmy3qjb0tSykCkRMTAxLly71OoYxxpQbIpJW1DK7xGSMMcYnKxDGGGN8sgJhjDHGp5C6B+FLbm4u6enp5OTkeB3llFWrVo2oqCjCw8O9jmKMqQBCvkCkp6dTq1YtYmJiEBGv45SYqpKZmUl6ejqxsbFexzHGVAAhf4kpJyeHBg0alOviACAiNGjQICTOhIwxpWTaNIiJgbAw533atFLdfcifQQDlvjgcEyo/hzGmFEybBgkJkJ3tTKelOdMA8fGl8hUhfwZhjDEhafTo/xWHY7KznfmlxAqEMcaUR5s3n9z8ErACYYwx5VF0NAAHqlRnafOOv5tfGqxAHC8AN33eeustunXrRvfu3RkyZAixsbHk5uYCsH//fmJiYsjNzeWFF16gU6dOdOvWjWuvvfaUv9cYE8ImTGB+57O54NZXue2KMWSHV4WICJgwodS+okLcpPZbAG76pKSkMGHCBBYvXkxkZCR79uzhH//4B7Nnz+byyy/n/fff58orryQ8PJyJEyeyceNGqlatyr59+0rnZzLGhJxdB47wsHRk9qUj6LBvKxNnPk5EsyZOcSilG9RgZxC/FYCbPl999RVXXXUVkZGRANSvX59bb72VqVOnAjB16lSGDx8OQLdu3YiPj+edd96hcmWr3caY31JVPli6hfOfXcj8lB3cd2E7Zr10K6dlrIZNm0q1OIAViN8KwE0fVf3d46l9+/Zl06ZNLFy4kPz8fLp06QLA7NmzueOOO1i2bBk9e/YkLy+vxN9rjAktaZmHuP71H3ngo0TaNa7JnLvP5s7z2lKlcuB+jVuBKKyomzuncNNnwIABfPDBB2RmZgKwZ88eAG688UaGDh3669lDQUEBW7ZsoX///jz11FPs27ePgwcPlvh7jTGhIS+/gMmL1nPRc4tYsSWLRy/vwn8S+tCmUc2Af7ddxyhswoTf3oOAU77p07lzZ0aPHs25555LpUqVOP3003njjTeIj49nzJgxDB06FID8/Hyuv/56srKyUFXuuece6tate4o/kDGmPEvJyGLk9CSStmZxfsdGPHp5F5rWqV5m328ForBj1+9Gj3YuK0VHl8pNn2HDhjFs2LDfzPvuu++46qqrfi0C4eHhfPfdd6f0PcaY0JCTm8/zC9YyedEG6kWE89J1p3NJ16Zl3puCFYjjxceX+o2e4/3tb3/j888/Z86cOQH9HmNM+fPDhkxGzUhi4+5DXN0zitGXdKRuRBVPsliB8MCLL77odQRjTJDJOpzLxM9X8d5Pm4muH8E7t/TmD20jPc0UsAIhIlOAS4GdqtrFx/L7gWP/Va8MdAQaquoeEdkEHADygTxVjTuVLL6eJCqPVNXrCMaYAPgieTsPfZLM7oNHuO3sWO69oD3Vq1TyOlZAzyDeAF4C3vK1UFWfBp4GEJHLgHtUdU+hVfqr6u5TDVGtWjUyMzPLfZffx8aDqFatmtdRjDGlZOf+HMbNSuHz5O10bFqb14edQdeoOl7H+lXACoSqLhKRGD9XHwq8F4gcUVFRpKens2vXrkDsvkwdG1HOGFO+HWvwNmH2SnLyCnhgYHtuO7sV4ZWCq+WB5/cgRCQCGAjcWWi2AvNERIF/qerkYrZPABIAon20VwgPD7cR2IwxQWPT7kOMmpHE9xsy6RVbn4lXdKVVw8C3aSgJzwsEcBmw+LjLS31VNUNEGgHzRWSVqi7ytbFbPCYDxMXF2UV6Y0xQyssv4N/fbuS5L9dQpVIYjw/pyrVntCAsLHgvfQdDgbiW4y4vqWqG+75TRGYCvQCfBcIYY4Jd8tYsRkxPJCVjPxd2asyjl3ehce3gv5/oaYEQkTrAucD1hebVAMJU9YD7+ULgEY8iGmNMiR0+ms9zC9bw2rcbqV+jCpOu78HALk29juW3QD7m+h7QD4gUkXRgHBAOoKqT3NWGAPNU9VChTRsDM90njioD76rqF4HKaYwxgfDfdbsZNTOJtMxshvZqwchBHalTPdzrWCclkE8xDfVjnTdwHoctPG8D0D0wqYwxJrCysnOZMCeVD5amE9Mggndv681Zrb1t8FZSwXAPwhhjyj1V5fPk7Tz0SQp7s49y+7mt+fv5bakW7n2Dt5KyAmGMMadox/4cxn6czLzUHXRpXps3hp9Bl+bB0+CtpKxAGGNMCRUUKO8t2czEOavILSjgwYs7cHPfWCoHWYO3krICYYwxJbB+10FGzUjip417OKt1Ax4f0pWYyBpexypVViCMMeYk5OYXMHnRBp5fsJZqlcN46spuXB0XVa77eiuKFQhjjPHTii37GDE9kVXbD3BJ16aM+2MnGtUK/gZvJWUFwhhjTiD7aB7PzlvDlMUbaVirKpNv6MmFnZt4HSvgrEAYY0wxFq3ZxYMzk0jfe5jrekczclAHalcrXw3eSsoKhDHG+LD30FEem72S6T+n0yqyBv9JOJPerRp4HatMWYEwxphCVJVPE7cxflYKWYdzubN/G+48r025bvBWUlYgjDHGlbHvMGM/TmbBqp10i6rD27f0plOz2l7H8owVCGNMhVdQoEz7MY0nv1hNfoEy5pKODO8bS6UgHquhLFiBMMZUaOt2HmDE9CSWpe3l7LaRPD6kKy3qR3gdKyhYgTDGVEhH8wp49Zv1vPz1OiKqVuKZq7tzRY/mIdngraSsQBhjKpyfN+9l5PRE1uw4yGXdm/HQpZ1oWKuq17GCjhUIY0yFcehIHv+ct5o3/ruJJrWr8fqwOAZ0bOx1rKBlBcIYUyF8s3ono2cms3XfYW7s05L7L2pPrQrS4K2krEAYY0LankNHefSzVGb+spXWDWvw0e19iIup73WsciFgnZaLyBQR2SkiyUUs7yciWSKy3H09VGjZQBFZLSLrRGRkoDIaY0KXqvLJ8q2c/+xCPkvM4K4BbZlz99lWHE5CIM8g3gBeAt4qZp1vVfXSwjNEpBLwMnABkA4sEZFZqpoaqKDGmNCSvjebMR8n883qXZzWoi5PXtmN9k1qeR2r3AlYgVDVRSISU4JNewHrVHUDgIi8DwwGrEAYY4qVX6C8/f0mnpq7GoBxl3Xixj4xFb7BW0l5fQ+ij4isADKA+1Q1BWgObCm0TjrQ24twxpjyY82OA4yYnsgvm/dxbruGTBjShah61uDtVHhZIH4GWqrqQRG5GPgYaAv4KvVa1E5EJAFIAIiOjg5ATGNMMDuSl88rX6/nlW/WUbNqZZ7702kMPq2ZNXgrBZ4VCFXdX+jzHBF5RUQicc4YWhRaNQrnDKOo/UwGJgPExcUVWUiMMaFnWZrT4G3tzoNcflozxl7aiQY1rcFbafGsQIhIE2CHqqqI9MJ5oioT2Ae0FZFYYCtwLXCdVzmNMcHn4JE8nv5iFW/9kEazOtWZOvwM+rdv5HWskBOwAiEi7wH9gEgRSQfGAeEAqjoJuAr4i4jkAYeBa1VVgTwRuROYC1QCprj3Jowxhq9X7WT0zCS27c9hWJ8Y7ruoPTWren07NTSJ8zs5NMTFxenSpUu9jmGMCYDdB4/wyKepzFqRQdtGNZl4ZTd6tqzndaxyT0SWqWqcr2VWdo0xQU1VmfHzVh6dncqhI3ncc347/tKvNVUqB6ydr3FZgTDGBK0te7J5cGYS367dTc+W9Zh4RVfaNrYGb2XFCoQxJujkFyhTF2/kmXlrCBN4ZHBnru/dkjBr8FamrEAYY4LKym37GTk9kRXpWZzXoRGPXd6FZnWrex2rQrICYYwJCjm5+bz01TomLVxPnerhvDD0dC7r1tQavHnICoQxxnM/bdzDyBmJbNh1iCt7RDHmko7Uq1HF61gVnhUIY4xnDuTkMvHzVUz7cTNR9arz1s29OKddQ69jGZcVCGOMJ+an7mDsx8nsPJDDLX+I5R8XtiOiiv1KCib2p2GMKVO7Dhzh4U9TmJ24jQ5NajHphp6c1qKu17GMD1YgjDFlQlX5cFk6E2av5PDRfO67sB1/Prc14ZWswVuwsgJhjAm4tMxDPDgzicXrMukVU5/Hr+hKm0Y1vY5lTsAKhDEmYPLyC5iyeCPPzl9D5bAwHru8C9f1irYGb+WEFQhjTECkZGQxcnoSSVuzOL9jYx67vAtN6lTzOpY5CVYgjDGlKic3n+cXrGXyog3Uiwjn5et6cHHXJtbgrRyyAmGMKTU/bMhk1IwkNu4+xDVxUTx4cUfqRliDt/LKCoQx5pRlHc5l4ucree+nLUTXj2Darb3p2ybS61jmFFmBMMacki+St/PQJ8nsPniEhHNacc/57ahepZLXsUwpsAJhjCmRnftzeOiTFL5I2U7HprV5fdgZdI2q43UsU4qsQBhjToqq8p8lW5gwZyVH8gp4YGB7bju7lTV4C0EBKxAiMgW4FNipql18LI8HRriTB4G/qOoKd9km4ACQD+QVNV6qMaZsbdx9iFEzEvlhwx56xdZn4hVdadXQGryFqhMWCBFpCzwBdAJ+fYhZVVudYNM3gJeAt4pYvhE4V1X3isggYDLQu9Dy/qq6+0T5jDGBl5tfwGvfbuS5L9dQpXIYT1zRlT/FtbAGbyHOnzOIqcA44P+A/sBw4IR/K1R1kYjEFLP8v4UmfwCi/MhijCljSelZjJieSOq2/Qzs3ITxgzvTuLY1eKsI/CkQ1VV1gYiIqqYBD4vItzhFo7TcAnxeaFqBeSKiwL9UdXJRG4pIApAAEB0dXYqRjKnYDh/N57kv1/DvbzcQWbMqk67vwcAuTb2OZcqQPwUiR0TCgLUiciewFWhUWgFEpD9OgfhDodl9VTVDRBoB80Vklaou8rW9WzwmA8TFxWlp5TKmIlu8bjejZiSxeU82Q3u1YOSgjtSpHu51LFPG/CkQfwcigLuAR3EuMw0rjS8XkW7Aa8AgVc08Nl9VM9z3nSIyE+gF+CwQxpjSk5Wdy4Q5qXywNJ3YyBq8d9uZ9GndwOtYxiMnLBCqugTAucKkw0vri0UkGpgB3KCqawrNrwGEqeoB9/OFwCOl9b3GmN9TVeYkbWfcrBT2Zh/lr/1ac9eAtlQLtwZvFZk/TzH1AV4HagLRItId+LOq/vUE270H9AMiRSQd555FOICqTgIeAhoAr7ideB17nLUxMNOdVxl4V1W/KNFPZ4w5oe1ZOYz9JJn5qTvo0rw2b958Bp2bWYM3498lpueAi4BZAKq6QkTOOdFGqjr0BMtvBW71MX8D0N2PXMaYU1BQoLy3ZDMT56wit6CABy/uwM19Y6lsDd6My6+Gcqq65biuevMDE8cYUxbW7zrIqOlJ/LRpD33bNODxIV1p2aCG17FMkPGnQGwRkbMAFZEqODerVwY2ljEmEHLzC5i8aAPPL1hLtcphPHVVN67uGWVjNRif/CkQtwPPA82BdGAecEcgQxljSt+KLfsYMT2RVdsPcEnXpoz7Yyca1bIGb6ZoxRYIEakEPKeq8WWUxxhTyrKP5vHsvDVMWbyRhrWqMvmGnlzYuYnXsUw5UGyBUNV8EWkoIlVU9WhZhTLGlI5v1+7iwZlJbNlzmPje0YwY1IHa1azBm/GPP5eYNgGLRWQWcOjYTFV9NlChjDGnZu+hozw2eyXTf06nVWQNPvhzH3rF1vc6liln/CkQGe4rDKgV2DjGmFOhqnyWuI3xn6awLzuXO/u34c7z2liDN1Mi/rSkHl8WQYwxpyZj32HGfpzMglU76R5Vh7dv6U3HprW9jmXKMX9aUjcEHgA689vxIM4LYC5jjJ8KCpRpP6bx5BeryS9QxlzSkeF9Y6lkYzWYU+TPJaZpwH9wRoe7Haejvl2BDGWM8c+6nQcYMT2JZWl7ObttJI8P6UqL+hFexzIhwp8C0UBVXxeRu1V1IbBQRBYGOpgxpmhH8wp49Zv1vPz1OiKqVuKZq7tzRY/m1uDNlCp/CkSu+75NRC7BuWFto78Z45GfN+9l5PRE1uw4yKXdmjLuss40rFXV61gmBPlTIB4TkTrAP4AXgdrAPQFNZYz5nUNH8nh67mre/H4TTWpX4/VhcQzo2NjrWCaE+fMU02fuxyycwYKMMWXs69U7GTMzmYysw9xwZkvuv6g9tazBmwkwv3pzNcZ4I/PgER79LJWPl2fQplFNPvxzH+JirMGbKRtWIIwJQqrKx8u38sinqRw8ksddA9pyR//WVK1sDd5M2bECYUyQSd+bzeiZySxcs4vTWtTlySu70b6JdWJgyl6RBUJE7i1uQ+uLyZjSlV+gvPX9Jp6euxqAhy/rxA19YqzBm/FMcWcQx/7L0h44A3fIUeAyYNGJdiwiU3Aa1+1U1S4+lgvOOBMXA9nATar6s7tsoLusEvCaqk7066cxppxas+MAD3yUyPIt+zi3XUMmDOlCVD1r8Ga8VWSBONYHk4jMA3qo6gF3+mHgQz/2/QbwEvBWEcsHAW3dV2/gVaC3OwbFy8AFOAMULRGRWaqa6sd3GlOuHMnL5+Wv1/PqN+uoWbUyz/3pNAaf1swavJmg4M89iGig8FgQR4GYE22kqotEpLj1BgNvqaoCP4hIXRFp6u57napuABCR9911rUCYkLIsbQ8jpiexbudBhpzenDGXdKRBTWvwZoKHPwXibeAnEZkJKDCEos8KTkZzYEuh6XR3nq/5vYvaiYgkAAkA0dHRpRDLmMA6eCSPp75Yxds/pNGsTnXeGH4G/do38jqWMb/jT0O5CSLyOXC2O2u4qv5SCt/t6xxai5lfVL7JwGSAuLi4ItczJhh8tWoHo2cms31/DjedFcN9F7anRlV7mNAEJ3//ZkYA+1V1qjsEaayqbjzF704HWhSajsLp56lKEfONKbd2HzzC+E9T+XRFBu0a1+Tl+LPoEV3P61jGFMuf8SDGAXE4TzNNBcKBd4C+p/jds4A73XsMvYEsVd0mIruAtiISC2wFrgWuO8XvMsYTqsqMn7fy6OxUso/kc+8F7bj93NZUqRzmdTRjTsifM4ghwOnAzwCqmiEiJ2y1IyLvAf2ASBFJB8bhFBdUdRIwB+cR13U4j7kOd5flicidwFycx1ynqGrKyf1Yxnhvy55sHpyZxLdrd9OzZT0mXtGVto2twZspP/wpEEdVVUVEAUSkhj87VtWhJ1iuwB1FLJuDU0CMKXfyC5SpizfyzLw1hAk8Orgz8b1bEmYN3kw540+B+EBE/gXUFZHbgJuB1wIby5jyaeW2/YycnsiK9CzO69CIxy7vQrO61b2OZUyJ+PMU0z9F5AJgP859iIdUdX7AkxlTjuTk5vPSV+uYtHA9daqH88LQ07msW1Nr8GbKNX9uUj+pqiOA+T7mGVPh/bRxDyNnJLJh1yGu7BHFmEs6Uq9GFa9jGXPK/HmU4gIf8waVdhBjypsDObmMnpnENf/6nqN5BbwVc5BnHhhMvVrVICYGpk3zOqIxp6S43lz/AvwVaCUiiYUW1QIWBzqYMcFsfuoOxn6czM4DOdzyh1j+sfcXIm5PgOxsZ4W0NEhIcD7Hx3sX1JhTIM7DRD4WOONQ1wOeAEYWWnRAVfeUQbaTFhcXp0uXLvU6hglhOw/kMH5WKrOTttGhSS0mXtmN01rUdc4Y0tJ+v0HLlrBpUxmnNMZ/IrJMVeN8LSuuN9csIEtExgDbVfWIiPQDuonIW6q6LxBhjQlGqsqHS9N5bHYqOXkF3H9RexLOaUV4Jfcq7ebNvjcsar4x5YA/j7lOB+JEpA3wOk4L6HdxGrkZE/LSMg/x4MwkFq/LpFdMfZ64siutG9b87UrR0b7PIKwDSVOO+VMgCtzWzVcAz6nqiyJSGp31GRPU8vILeP27jfzfl2sIDwtjwpAuDD0j2neDtwkTnHsOx+5BAEREOPONKaf8KRC5IjIUuBFnNDlwu8wwJlSlZGQxYnoiyVv3c0Gnxjw6uAtN6lQreoNjN6JHj3YuK0VHO8XBblCbcsyfAjEcuB2YoKob3U703glsLGO8kZObz/ML1jJ50QbqRVThlfgeDOrSxL8Gb/HxVhBMSPGnJXUqcFeh6Y2AjRFtQs736zMZNSORTZnZXBMXxYMXd6RuhDV4MxWXPy2pN+JjwB5VbRWQRMaUsazDuTwxZyXvL9lCdP0Ipt3am75tIr2OZYzn/LnEVPj52GrA1UD9wMQxpmx9kbyNsZ+kkHnwCH8+pxV/P78d1atU8jqWMUHBn0tMmcfNek5EvgMeCkwkYwJvx/4cxn2Swhcp2+nUtDZThp1B16g6XscyJqj4c4mpR6HJMJwzChv1xJRLqsr7S7bw+JyVHM0rYMTADtx6duz/GrwZY37lzyWmZwp9zgM2AdcEJI0xAbRx9yFGzUjkhw176B1bn4lXdiM20q/xr4ypkPy5xNS/LIIYEyi5+QW89u1GnvtyDVUqhzHxiq5cE9fCRngz5gSK6831elV9R0Tu9bFYgT3ALFXdW8w+BgLP44wt/ZqqTjxu+f3AsQfHKwMdgYaqukdENgEHgHwgr6jOpIwpTlK60+Atddt+BnZuwvjBnWlcu5gGb8aYXxV3BnHs3Luo+w2xwF+AM30tFJFKwMs440mkA0tEZJbbrgIAVX0aeNpd/zLgnuN6iu2vqrv9+UGMKezw0Xye+3IN//52Aw1qVmXS9T0Y2KWp17GMKVeK6831X+77+KLWEZFHitl3L2Cdqm5w130fGAykFrH+UOC9EwU25kT+u243o2YmkZaZzbVntGDUxR2pU916hzHmZBV3iemF4jZU1btUtbhHXZsDWwpNpwO9i/iuCGAgcGfhrwDmiYgC/1LVyUVsmwAkAERbz5kVWlZ2LhPmpPLB0nRiGkTw7m29Oau1NXgzpqSKu8S0zH3vC3QC/uNOX11oWXF83QH0PTqR0wng4uMuL/VV1QwRaQTMF5FVqrrodzt0CsdkcAYM8iOXCTGqyufJ23nokxT2Zh/lL/1ac/eAtlQLtwZvxpyK4i4xvQkgIjfh3AvIdacnAfP82Hc60KLQdBSQUcS613Lc5SVVzXDfd4rITJxLVr8rEKZi256Vw9hPkpmfuoMuzWvz5s1n0LmZNXgzpjT40w6iGc6N6mP/u6/pzjuRJUBbt/fXrThF4LrjV3KHNj0XuL7QvBpAmKoecD9fCBR3v8NUMAUFyntLNjNxzipyCwoYNagDt/whlsrW4M2YUuNPgZgI/CIiX7vT5wIPn2gjd5ChO4G5OI+5TlHVFBG53V0+yV11CDBPVQ8V2rwxMNPtYrky8K6qfuFHVlMBrN91kFHTk/hp0x7Oat2AJ67oSssG1uDNmNImqie+bC8iTfjfDeYfVXV7QFOVUFxcnC5dutTrGCZAjuYVMHnRel74ah3VKocx5tJOXN0zyr+xGowxPonIsqLamflzBgFwBNiG05trOxFp5+uGsTGBsmLLPkZMT2TV9gNc0rUp4/7YiUa1rMGbMYHkT2d9twJ349xkXo7TMO574LyAJjMGyD6axzPz1jB18UYa1qrK5Bt6cmHnJl7HMqZC8OcM4m7gDOAHVe0vIh2AIhvPGVNaFq3ZxYMzk0jfe5jrz4zmgYEdqF3NGrwZU1b8KRA5qpojIohIVVVdJSLtA57MVFh7Dx3l0dmpzPh5K60a1uCDP/ehV6yNUWVMWfOnQKSLSF3gY5wGa3spuj2DMSWmqnyauI3xs1LIOpzL385rwx3921iDN2M84k9330Pcjw+7j7rWAeyRU1OqMvYdZuzHySxYtZPuUXV459bedGxa2+tYxlRo/j7FBICqLgxUEFMxFRQo7/yYxpOfr6JAYeylnbjprBgq2VgNxnjupAqEMaVp3c4DjJiexLK0vZzdNpLHh3SlRf0Ir2MZY1xWIEyZO5pXwKvfrOflr9cRUbUSz17TnSGnN7cGb8YEGSsQpkz9vHkvI6cnsmbHQS7r3oxxl3UismZVr2MZY3ywAmHKxKEjeTw9dzVvfr+JJrWr8fqwOAZ0bOx1LGNMMaxAmID7evVOxsxMJiPrMDee2ZL7B3agZlX7q2dMsLN/pSZg9hw6yiOfpvDx8gzaNKrJR7f3oWdLa/BmTHlhBcKUOlXlk+UZPPJZKgdycrlrQFvu6N+aqpWtwZsx5YkVCFOq0vdmM3pmMgvX7OK0FnV58sputG9Sy+tYxpgSsAJhSkV+gfLW95t4eu5qAMZd1okb+1iDN2PKMysQ5pSt2XGAEdMT+WXzPs5t15AJQ7oQVc8avBlT3lmBMCV2JC+fl79ez6vfrKNm1co896fTGHxaM2vwZkyICOgI7yIyUERWi8g6ERnpY3k/EckSkeXu6yF/tzXeWpa2h0te+I4XFqzlkq5N+fLec7ncWkMbE1ICdgYhIpWAl4ELgHRgiYjMUtXU41b9VlUvLeG2powdyMnl6bmrefuHNJrVqc7U4WfQv30jr2MZYwIgkJeYegHrVHUDgIi8DwwG/PklfyrbmgBZsHIHYz5OZvv+HIb1ieH+i9pTwxq8GROyAvmvuzmwpdB0OtDbx3p9RGQFziBE96lqyklsi4gkAAkA0dHRpRDbHG/3wSOM/zSVT1dk0K5xTV6OP4se0fW8jmWMCbBA3oPwdTFaj5v+GWipqt2BF3FGrfN3W2em6mRVjVPVuIYNG5Y0qyls2jSIiUHDwvio3584/4l5zE3ezr0XtOOzv51txcGYCiKQBSIdaFFoOorjhipV1f2qetD9PAcIF5FIf7Y1ATJtGiQksGVfDjdePZ77zryR1lvWMLvVPu4a0JYqlQP6XIMxJogE8hLTEqCtiMQCW4FrgesKryAiTYAdqqoi0gunYGUC+060rQmM/DFjmNrpAp45+wYqaQGPznuF+F8+J2xRNNwy1Ot4xpgyFLACoap5InInMBeoBExR1RQRud1dPgm4CviLiOQBh4FrVVUBn9sGKqtxrNy2n5Hn3MWKpu0YsO4nHp33Cs0O7HYWbt7sbThjTJkT5/dxaIiLi9OlS5d6HaPcycnN58Wv1vKvhRuok72fh+e+wqWrvv3tjaCWLWHTJo8SGmMCRUSWqWqcr2X2jGIF9+OGTEbNSGLD7kNc1TOK0dn7qPfvZb9dKSICJkzwJqAxxjNWICqo/Tm5PPn5Kqb9uJkW9avz9i29OLttQ6A7VFYYPdq5rBQd7RSH+HivIxtjypgViApoXsp2xn6SzK4DR7j1D7Hce2E7IqoU+qsQH28FwRhjBaIi2Xkgh/GzUpmdtI0OTWox+YY4ureo63UsY0yQsgJRAagqHy5N57HZqeTkFXD/Re1JOKcV4ZWsTYMxpmhWIEJcWuYhHpyZxOJ1mfSKqc8TV3aldcOaXscyxpQDViBCVF5+AVMWb+TZ+WsIDwvjscu7cF2vaMJshDdjjJ+sQISglIwsRkxPJHnrfi7o1JhHB3ehSZ1qXscyxpQzViBCSE5uPs8vWMvkRRuoF1GFV+J7MKhLExvExxhTIlYgQsT36zMZNSORTZnZXBMXxYMXd6RuRBWvYxljyjErEOVc1uFcnpizkveXbCG6fgTTbu1N3zaRXscyxoQAKxDl2BfJ2xj7SQqZB4/w53Na8ffz21G9SiWvYxljQoQViHJox/4cHvokmbkpO+jUtDZThp1B16g6XscyxoQYKxDlSEGB8p+lW3h8zkqO5hUwYmAHbj071hq8GWMCwgpEObFh10FGzUjix4176B1bn4lXdiM2sobXsYwxIcwKRJDLzS9g8qINPL9gLVUrhzHxiq5cE9fCGrwZYwLOCkQQS0zfx4jpSazctp9BXZow/o+daVTbGrwZY8qGFYggdPhoPv/35Rpe+3YDkTWrMun6ngzs0sTrWMaYCiagdzdFZKCIrBaRdSIy0sfyeBFJdF//FZHuhZZtEpEkEVkuIqE1jui0aRATA2Fhzvu0ab8u+m7tbi56bhGTF23gT2e0YP6951pxMMZ4ImBnECJSCXgZuABIB5aIyCxVTS202kbgXFXdKyKDgMlA70LL+6vq7kBl9MS0aZCQANnZznRaGiQksC9PmBDRmQ+XpRMbWYP3E87kzFYNvM1qjKnQAnmJqRewTlU3AIjI+8Bg4NcCoar/LbT+D0BUAPMEh9Gj/1ccAAXmtDidcT8Le2ts5a/9WnPXgLZUC7cGb8YYbwWyQDQHthSaTue3ZwfHuwX4vNC0AvNERIF/qepkXxuJSAKQABAdHX1KgcvE5s2/ftxWqwFjL/gLX7Y9k67b1/LmiIvp3MwavBljgkMgC4Sv5zDV54oi/XEKxB8Kze6rqhki0giYLyKrVHXR73boFI7JAHFxcT73H1SioylI28y7pw1kYr/h5IWF8eDXr3PzzuVUbvZ3r9MZY8yvAlkg0oEWhaajgIzjVxKRbsBrwCBVzTw2X1Uz3PedIjIT55LV7wpEebP+oScY9XU6PzXvRN9Ny3l87ku0PLofJvs8QTLGGM8E8immJUBbEYkVkSrAtcCswiuISDQwA7hBVdcUml9DRGod+wxcCCQHMGvAHc0r4KWv1jJofV1WR3fgqR/f5p0PxtKyTlWnOMTHex3RGGN+I2BnEKqaJyJ3AnOBSsAUVU0Rkdvd5ZOAh4AGwCvuoDZ5qhoHNAZmuvMqA++q6heByhpoK7bsY8T0RFZtP8Al3Zoy7rJONKp1GfC+19GMMaZIohr8l+39FRcXp0uXBk+TieyjeTwzbw1TF2+kUa1qPHp5Fy7o1NjrWMYY8ysRWeb+x/x3rCV1gCxas4sHZyaRvvcw158ZzQMDO1C7WrjXsYwxxm9WIErZ3kNHeXR2KjN+3kqrhjX48PY+nBFT3+tYxhhz0qxAlBJV5dPEbYyflULW4Vz+dl4b7ujfxhq8GWPKLSsQpSBj32HGfpzMglU76R5Vh3du7U3HprW9jmWMMafECsQpKChQ3vkxjSc/X0WBwthLO3HTWTFUsrEajDEhwAqEv6ZNc/pR2rwZoqNZO/YJRh5tybK0vZzdNpLHh3SlRf0Ir1MaY0ypsQLhj0I9sB4Nq8yrzfvw8qrqRFTJ5NlrTmfI6c1x22wYY0zIsALhD7cH1mXNOjBq4N9Y07All6UuZNzqOUROKNcNvI0xpkhWIPxwcPsu/jkggTd7XkqTA5m8/tF4BqxfAnbWYIwJYVYgTuDr1TsZc9skMiLqcuPPs7l/0VvUPHrYWVgeuhc3xpgSsgJRhMyDR3jks1Q+WZ5Bm/p1+Oidh+i5Yfn/VoiIgAkTPMtnjDGBFtAxqcuF48aH1nemMfOXdM5/diFzkrZx94C2zB57CT0fuQ9atnQuK7VsaT2wGmNCXsU+gzhufOj0vdmMnr2Rhcl1OT26Lk9e2Y12jWs568bHW0EwxlQoFbtAuE8n5UsYb/a4lH+ecwMADy/7gBsef8MavBljKrSKXSA2byarag2GXTOe5c060G/9UibMfZnmB3dD2JtepzPGGE9V7AIRHU3ttDRa7t3O8KWf8seVC52BtFu29DqZMcZ4rmLfpJ4wAYmI4PnP/sngY8XBnk4yxhigoheI+HjnaSR7OskYY34noAVCRAaKyGoRWSciI30sFxF5wV2eKCI9/N221MTHw6ZNUFDgvFtxMMYYIIAFQkQqAS8Dg4BOwFAR6XTcaoOAtu4rAXj1JLY1xhgTQIE8g+gFrFPVDap6FHgfGHzcOoOBt9TxA1BXRJr6ua0xxpgACmSBaA5sKTSd7s7zZx1/tgVARBJEZKmILN21a9cphzbGGOMIZIHw1cpM/VzHn22dmaqTVTVOVeMaNmx4khGNMcYUJZDtINKBFoWmo4AMP9ep4se2xhhjAiiQZxBLgLYiEisiVYBrgVnHrTMLuNF9mulMIEtVt/m5rTHGmAAK2BmEquaJyJ3AXKASMEVVU0Tkdnf5JGAOcDGwDsgGhhe37Ym+c9myZbtFJO242ZHA7lL6sUpbMGeD4M5n2UommLNBcOcL1WxFdh0hqj4v7YcMEVmqqnFe5/AlmLNBcOezbCUTzNkguPNVxGwVuyW1McaYIlmBMMYY41NFKBCTvQ5QjGDOBsGdz7KVTDBng+DOV+Gyhfw9CGOMMSVTEc4gjDHGlIAVCGOMMT6FbIEos+7CT4KIbBKRJBFZLiJL3Xn1RWS+iKx13+uVUZYpIrJTRJILzSsyi4iMco/lahG5yINsD4vIVvfYLReRiz3K1kJEvhaRlSKSIiJ3u/OD5dgVlc/z4yci1UTkJxFZ4WYb7873/NgVk83z41bo+yqJyC8i8pk7Hfjjpqoh98JpXLceaIXTbccKoFMQ5NoERB437ylgpPt5JPBkGWU5B+gBJJ8oC06X6yuAqkCse2wrlXG2h4H7fKxb1tmaAj3cz7WANW6GYDl2ReXz/Pjh9LFW0/0cDvwInBkMx66YbJ4ft0LfeS/wLvCZOx3w4xaqZxDlqbvwwcCb7uc3gcvL4ktVdRGwx88sg4H3VfWIqm7Eafneq4yzFaWss21T1Z/dzweAlTg9DQfLsSsqX1HKLJ86DrqT4e5LCYJjV0y2opTpn6uIRAGXAK8dlyGgxy1UC4Tf3YWXMQXmicgyEUlw5zVWp/8p3PdGnqUrOkuwHM87xRl5cEqh02nPsolIDHA6zv82g+7YHZcPguD4uZdJlgM7gfmqGjTHrohsEATHDXgOeAAoKDQv4MctVAuE392Fl7G+qtoDZ6S8O0TkHK8D+SkYjuerQGvgNGAb8Iw735NsIlITmA78XVX3F7eqj3le5AuK46eq+ap6Gk4Pzb1EpEsxqwdDNs+Pm4hcCuxU1WX+buJjXomyhWqB8Ker8TKnqhnu+05gJs5p3w5xRtHDfd/pXcIis3h+PFV1h/sPuAD4N/87ZS7zbCISjvPLd5qqznBnB82x85UvmI6fm2cf8A0wkCA6dsdnC5Lj1hf4o4hswrlcfp6IvEMZHLdQLRBB1124iNQQkVrHPgMXAslurmHuasOAT7xJCMVkmQVcKyJVRSQWZwzxn8oy2LF/CK4hOMeuzLOJiACvAytV9dlCi4Li2BWVLxiOn4g0FJG67ufqwPnAKoLg2BWVLRiOm6qOUtUoVY3B+V32lapeT1kct0DedffyhdON+BqcO/ijgyBPK5wnC1YAKccyAQ2ABcBa971+GeV5D+eUORfnfxy3FJcFGO0ey9XAIA+yvQ0kAYnuP4CmHmX7A87peiKw3H1dHETHrqh8nh8/oBvwi5shGXjoRP8GgiCb58ftuJz9+N9TTAE/btbVhjHGGJ9C9RKTMcaYU2QFwhhjjE9WIIwxxvhkBcIYY4xPViCMMcb4ZAXCmDIiIjEicl0JtqsrIn8tNN1MRD4q3XTG/J4VCGPKTgzgs0CISOVitqsL/FogVDVDVa8q1WTG+GAFwoQ8EbnR7WxthYi8LSItRWSBO2+BiES7670hIi+IyH9FZIOIXFVoHw+IM5bHChGZ6M5rLSJfuJ0vfisiHU6wn4nA2e64AveIyE0i8qGIfIrTiWNNN8/P7ncNLrRda3e7p90zkWT3u6qJyFR3/V9EpL87/yYRmeHmWysiT5XJwTahpSxa/9nLXl69gM44rUkj3en6wKfAMHf6ZuBj9/MbwIc4/3HqhNNlPDidK/4XiDi2D/d9AdDW/dwbpwuE4vbTD7cVrDt9E05L8WP7qwzUdj9H4nTTLDhnHoXHxvh1GvgHMNX93AHYDFRz970BqONOpwEtvP7zsFf5ehV3WmtMKDgP+EhVdwOo6h4R6QNc4S5/G2fglWM+VqdjtlQRaezOOx/nl3B2oX3UBM4CPnS6PwKcAVqK248v81X12NgXAjzu9vJbgNNFc3HbgtO1xoturlUikga0c5ctUNUsABFJBVry226gjSmWFQgT6oQTd3VcePmR47Ytah9hwD51uof2xdd+fDlU6HM80BDoqaq5bu+d1YrZ9kT7LpwhH/v3bk6S3YMwoW4BcI2INABnHF+cy0XXusvjge9OsI95wM0iEnFsH+qMsbBRRK5254mIdD/Bfg7gDANalDo4/f7nuvcSWvqx3SL3Z0BE2gHROJfUjDllViBMSFPVFGACsFBEVgDPAncBw0UkEbgBuPsE+/gCpyfPpeKMOHafuygeuMXdbwonHtY2Echzb3Tf42P5NCBORJa6+17lfn8msFhEkkXk6eO2eQWoJCJJwH+Am1T1CMaUAuvN1RhjjE92BmGMMcYnKxDGGGN8sgJhjDHGJysQxhhjfLICYYwxxicrEMYYY3yyAmGMMcan/weiSrr8GCpb+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r-squared cys =  0.9988046024941847\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArsElEQVR4nO3dd3xUVfrH8c+TQgIECCVABELoLRQhFGFVsIFdYN1dRRQLWF8iNsoqYu+LUdfCqlh+uja6CIgooqIoqCQh9B5aaKGFhCTz/P6YCRswhCFkcqc879crr7lzp9znAPly58w554qqYowxJnSEOV2AMcaYimXBb4wxIcaC3xhjQowFvzHGhBgLfmOMCTERThfgjTp16mhiYqLTZRhjTEBZsmTJLlWNO35/QAR/YmIiixcvdroMY4wJKCKysaT91tVjjDEhxoLfGGNCjAW/McaEmIDo4y9Jfn4+mZmZ5ObmOl1KuYmOjqZhw4ZERkY6XYoxJogFbPBnZmZSrVo1EhMTERGnyzltqsru3bvJzMykSZMmTpdjjHGYqh6TbcffPx0B29WTm5tL7dq1gyL0AUSE2rVrB9UnGGNM2YybP44Rc0ZQtIimqjJizgjGzR9XLu8fsMEPBE3oFwm29hhjTp2qkp2bTcqilKPhP2LOCFIWpZCdm015rKgcsF09xhgTjESE8X3HA5CyKIWURSkADO8+nPF9x5fLCWJAn/E7bcOGDSQlJTldhjEmyBQP/yLlFfoQQsF//McjuwCNMcZfFXXvRLoagyeqivf5n66QCH5fflFSWFjI0KFDadeuHRdddBGHDx9mzZo1XHDBBXTs2JHOnTuzdu1aBg8ezLRp046+btCgQUyfPv20j2+MCS6qyl0zH+CD7ws4I+/fzP7rVoZ3H35Mn//pCvrg9/UXJatXr+bOO+9k2bJlxMbGMmnSJAYNGsSdd97J0qVLWbhwIfHx8dxyyy1MnDgRgH379rFw4UIuueSS8miiMSaIfJWxg29+7UW1wn4MPbsJZ7eMY3zf8QzvPpzY6Nhy6e4J+i93ff1FSZMmTejUqRMAXbp0Yf369WzZsoX+/fsD7klZAOeeey533nknWVlZTJ48mYEDBxIREfR//MYYL2Xtz+WR6cuYlb6dNvF1eGZAEh0b1Tz6eHn28YdE8hSFf1HoQ/n9IUZFRR3dDg8PZ+/evSd87uDBg/nwww/5+OOPeeedd0772MaYwKeqfLp4M0/OXE5ugYsH+rZi2DlNiQw/tkOmPId7B31XD/yvT7+48vyipLjq1avTsGFDpk6dCkBeXh45OTkADBkyhJdeegmAdu3alfuxjTGBZcOuQ1z7n0WMnJRG6/jqzB5+Nnf2af6n0C9vQR/8xfv0h3cfjmusq9y/KDneBx98wMsvv0yHDh3o2bMn27dvB6BevXq0adOGG2+8sdyPaYwJHAWFLt74bi19X1pA+pZ9PNW/PR8P7UHTuJgKOX7Qd/WICLHRscf06Rf1+Z/uFyWJiYmkp6cfvX///fcf3f7mm2/+9PycnBxWr17NNddcU+ZjGmMCW/qWfYyclMqyrfu5qG09Hr8qiXrVoyu0hqAPfoBxvccds8BRUfhX5BIJX3/9NTfddBP33nsvNWrUqLDjGmP8w+Ejhbw0bxVvfb+eWlUr8cZ1nemXFO9ILSER/PDnL0Yqel2cCy64gE2bNlXoMY0x/mHhml2MnpLGxt05XNOtEaMubkONys4tvx7QwV+ey5T6A5tNbExw2ZeTz5NfZvDp4kwSa1fho6Hd6dmsjtNlBW7wR0dHs3v37qBZmrloPf6icf/GmMClqsxK387YacvYm3OE23s3Y/j5LYiODHe6NCCAg79hw4ZkZmayc+dOp0spN0VX4DLGBK7t+3J5eFo6czN2kNSgOu/e2JWkBv71vV7ABn9kZKRdqcoY4zdcLuWjXzbx7KwV5LtcjLmkNTf1akKEj8fkl0XABr8xxviLtTsPMnpSGr9s2EPPZrV5ekB7Gteu6nRZJ2TBb4wxZZRf6OLN79by8jdriI4I47m/duDqLg39/ntHnwW/iDQC3gfqAy5ggqqmiMg4YChQ1Dk/RlW/9FUdxhjjC0s3ZzNyUiorth/g0vbxPHJFW+pWC4zBGb484y8A7lPV30SkGrBEROZ6Hhuvqi/48NjGGOMTOUcKePGrVUz8cT1x1aKYMLgLF7Wr73RZp8Rnwa+q24Btnu0DIrIcaOCr4xljjK8tWLWTMVPSyNx7mOt6JPBgv9ZUj3ZuIlZZVUgfv4gkAmcCi4BewF0icj2wGPengj+tZSwiw4BhAAkJCRVRpjHGlGjvoSM8PjODyb9toWlcVT699Sy6NanldFllJr6eLSoiMcB3wJOqOllE6gG7cF9J8nEgXlVvKu09kpOTdfHixT6t0xhjjqeqzEjdxqPTl7HvcD63927GnX2a+81ErJMRkSWqmnz8fp+e8YtIJDAJ+FBVJwOo6o5ij/8H+MKXNRhjTFlsyT7Mw1PT+WZFFh0bxfLhwPa0rl/d6bLKhS9H9QjwNrBcVf9VbH+8p/8foD+QXtLrjTHGCS6X8sHPG3lu9gpcCg9f1pYhPRMJD/PvIZqnwpdn/L2AwUCaiPzh2TcGuEZEOuHu6tkA3OrDGowxxmurdxxg5KRUftuUzTkt43jyqiQa1aridFnlzpejen4ASvov0sbsG2P8ypECF6/NX8O/v11DTFQE4//ekas6NfD7iVhlZTN3jTEhbcnGvYyalMrqrINc2ekMxl7WltoxUU6X5VMW/MaYkHQwr4DnZ6/g/Z83El89molDutKndV2ny6oQFvzGmJDzzYodPDQlnW37c7nhrETu79uKmKjQicPQaakxJuTtOpjHYzMymL50Ky3qxvD5bT3p0rim02VVOAt+Y0zQU1Um/7aFx2dmcCivgHsuaMHtvZsRFREYE7HKmwW/MSaobd6Tw5gpaXy/ehedE2J5ZmAHWtar5nRZjrLgN8YEpUKXMvHH9bz41SrCBB67sh3XdW9MWBBNxCorC35jTNBZsX0/Iz9PZWnmPvq0iuOJ/u1pEFvZ6bL8hgW/MSZo5OYX8uo3a3jju7XUqBzJy9ecyeUd4oN2IlZZWfAbY4LCL+v3MGpyKut2HmJA5wY8dGlbalWt5HRZfsmC3xgT0Pbn5vPsrBV8uGgTDWtW5v2bunFOyziny/JrFvzGmIA1N2MHD09NJ+tALrf8pQn3XtSSKpUs1k7G/oSMMQEn60Auj07PYGbaNlrXr8Ybg7vQqVGs02UFDAt+Y0zAUFU+W5zJEzMzyC1w8UDfVgw9uymVIsKcLi2gWPAbYwLCxt2HGD05jYVrd9OtSS2eHtCeZnExTpcVkCz4jTF+raDQxds/rOdfc1dRKTyMp/q35x9dG9lErNNgwW+M8VvpW/YxclIqy7bu58K29Xj8yiTq14h2uqyAZ8FvjPE7ufmFjP96FW99v55aVSvx2qDOXJxU3yZilRMLfmOMX1m4dhejJ6excXcOf09uxJhL2lCjSqTTZQUVC35jjF/Yl5PPU18u55PFm2lcuwof3dKdns3rOF1WULLgN8Y4blbaNsZOX8aeQ0e49dym3HN+SypXCs218iuCBb8xxjE79ucydlo6c5btoN0Z1Zk4pCtJDWo4XVbQs+A3xlQ4l0v5+NfNPP3lco4Uuhh1cWtu+UsTIsJtIlZFsOA3xlSodTsPMnpyGovW7+GsprV5ekB7EutUdbqskGLBb4ypEPmFLiYsWEfKvNVER4Tx3MAOXJ3c0IZoOsCC3xjjc6mZ2Tz4eSorth/gkvb1GXd5O+pWt4lYTrHgN8b4TM6RAsbPXcXbP6wnrloUbw7uQt929Z0uK+RZ8BtjfOKH1bsYPSWVzXsOc233BEZd3Jrq0TYRyx/4LPhFpBHwPlAfcAETVDVFRGoBnwCJwAbgb6q611d1GGMq1t5DR3hi5nIm/ZZJ0zpV+WRYD7o3re10WaYYX57xFwD3qepvIlINWCIic4EhwDxVfUZERgGjgJE+rMMYUwFUlS9St/HojGVk5+RzV5/m3HVec6IjbSKWv/FZ8KvqNmCbZ/uAiCwHGgBXAr09T3sPmI8FvzEBbWv2YR6ems68FVl0aFiD92/qTtszqjtdljmBCunjF5FE4ExgEVDP858CqrpNROqe4DXDgGEACQkJFVGmMeYUuVzK/y3ayLOzVuBSeOjSNtzYqwnhtla+X/N58ItIDDAJuEdV93s7ZldVJwATAJKTk9V3FRpjymL1jgOMmpzGko17ObtFHZ7q355Gtao4XZbxgk+DX0QicYf+h6o62bN7h4jEe87244EsX9ZgjClfRwpcvD5/Lf/+dg1VosJ54eqODOzcwCZiBRBfjuoR4G1guar+q9hD04EbgGc8t9N8VYMxpnwt2biX0ZNTWbXjIFd0PIOxl7elTkyU02WZU+TLM/5ewGAgTUT+8OwbgzvwPxWRm4FNwNU+rMEYUw4O5hXwwpyVvPfTBuKrR/POkGTOa13P6bJMGflyVM8PwIk++53vq+MaY8rXtyuzeGhKOlv3Heb6Ho15oF9rYqJs7mcgs789Y0yJdh/M47EvMpj2x1aa143h89vOokvjWk6XZcqBBb8x5hiqytQ/tvDYjAwO5hUw/PwW3NGnGVERNhErWFjwG2OO2rwnh39OTWfBqp2cmRDLswM70LJeNafLMuXMgt8YQ6FLeXfhBl6YsxIRGHd5WwaflWgTsYLUSYNfRFoATwNtgaMLaKtqUx/WZYypICu272fkpDSWbs6md6s4nrgqiYY1bSJWMPPmjH8i8AgwHugD3MiJR+sYYwJEXkEhr36zhtfnr6V65UhS/tGJKzqeYROxQoA3wV9ZVeeJiKjqRmCciHyP+z8DY0wA+nXDHkZNSmXtzkMMOLMBD13WllpVKzldlqkg3gR/roiEAatF5C5gC1DiwmrGGP92IDef52av5IOfN9IgtjLv3dSNc1vGOV2WqWDeBP89QBXgbuBx3N09N/iwJmOMD3ydsYOHpqaTdSCXm//ShHsvbElVm4gVkk76t66qvwK4e3r0Rt+XZIwpTzsP5DFuxjJmpm6jdf1qvDG4C50axTpdlnGQN6N6zsK92FoMkCAiHYFbVfUOXxdnjCk7VeXzJZk8MXM5h48Ucv9FLRl2TjMqRYQ5XZpxmDef814C+uJeVRNVXSoi5/iyKGPM6dm0O4fRU1L5cc1uuiXW4qkB7WleN8bpsoyf8KqDT1U3HzfEq9A35RhjTkdBoYuJP27gxbkriQgL44mrkri2WwJhNhHLFONN8G8WkZ6Aikgl3F/yLvdtWcaYU7Vs6z5GTUojbcs+LmhTjyeuSqJ+jeiTv9CEHG+C/zYgBfeF0jOBr4A7fVmUMcZ7ufmFpMxbzYQF66hZpRKvDerMxUn1bSKWOaFSg19EwoGXVHVQBdVjjDkFP63dzejJqWzYncPfkhsy5pI2xFaxiVimdKUGv6oWikiciFRS1SMVVZQxpnT7Dufz9JfL+fjXzSTUqsKHt3SnV/M6TpdlAoQ3XT0bgB9FZDpwqGjncdfRNcZUkNnp23h42jJ2H8zj1nOacs8FLalcydbKN97zJvi3en7CAFuY2xiH7Nify9hp6cxZtoO28dWZOKQrSQ1qOF2WCUDezNx9tCIKMcaUzOVSPlm8mae+XM6RAhcj+7XmlrObEBluE7FM2XgzczcOeBBox7Hr8Z/nw7qMMcC6nQcZPTmNRev30KNpLZ4e0IEmdao6XZYJcN509XwIfAJchnto5w3ATl8WZUyoyy90MWHBOlLmrSYqIoxnBrTn710b2RBNUy68Cf7aqvq2iAxX1e+A70TkO18XZkyoSsvcx4OTUlm+bT8XJ9Xn0SvaUbe6TcQy5ceb4M/33G4TkUtxf9Hb0HclGROaco4UMH7uKt7+YT11YqJ447ou9Euq73RZJgh5E/xPiEgN4D7gFaA6MMKnVRkTYn5YvYvRU1LZvOcw13RLYNTFralROdLpskyQ8mZUzxeezX24L8JijCkn2TlHeGLmcj5fkkmTOlX5eFgPejSt7XRZJsjZ5XeMcYCqMjNtG+OmL2NvTj539G7G3ee3IDrSJmIZ37PgN6aCbdt3mIenpvP18izaN6jB+zd1p+0Z1Z0uy4QQnwW/iLyDewholqomefaNA4byv+GgY1T1S1/VYIw/cbmUDxdt5NnZKylwuXjo0jYM6ZlIhE3EMhXshMEvIveW9kIv1up5F3gVeP+4/eNV9QWvqjMmSKzJOsCoSWks3riXs1vU4cmr2pNQu4rTZZkQVdoZf9G6PK2ArnguvQhcDiw42Rur6gIRSTyt6owJcEcKXLzx3Vpe/WYNlSuF88LVHRnYuYFNxDKOOmHwF63RIyJfAZ1V9YDn/jjgs9M45l0icj2wGLhPVfeW9CQRGQYMA0hISDiNwxnjjN827WXUpFRW7TjIZR3ieeTydsRVi3K6LGPwpnMxASi+Fv8RILGMx3sdaAZ0ArYBL57oiao6QVWTVTU5Li6ujIczpuIdyivg0RnLGPj6Qg7kFvD2Dcm8em1nC33jN7z5cvcD4BcRmQIo0J8/99t7RVV3FG2LyH+AL0p5ujEBZ/7KLP45JZ2t+w5zfY/GPNCvNTFRNnjO+BdvJnA9KSKzgLM9u25U1d/LcjARiVfVbZ67/YH0sryPMf5mz6EjPDZjGVP/2ErzujF8fttZdGlcy+myjCmRt6ciVYD9qjrRcynGJqq6vrQXiMh/gd5AHRHJBB4BeotIJ9yfHDYAt5a1cGP8gaoy7Y+tPPZFBgdy8xl+fgvu6NOMqAibiGX8lzfr8T8CJOMe3TMRiAT+D+hV2utU9ZoSdr9dhhqN8UuZe3N4aGo681fu5MyEWJ4d2IGW9ewidcb/eXPG3x84E/gNQFW3ioj96zYhq9ClvP/TBp6fsxKAcZe3ZfBZiYSH2RBNExi8Cf4jqqoiogAiYpf/MSFr1Y4DjJyUyu+bsundKo4nrkqiYU2biGUCizfB/6mIvAnEishQ4CbgLd+WZYx/ySso5N/fruX1+WuIiYrgpb934spOZ9hELBOQvBnV84KIXAjsx93PP1ZV5/q8MmP8xJKNexg5KY01WQfpf2YDHrq0DbVjbEy+CVzefLn7rKqOBOaWsM+YoHUgN5/n56zkg583ckaNyrx7Y1d6t6rrdFnGnDZvunouBI4P+YtL2GdM0Ji3fAcPTU1n+/5cbuyZyH0XtaKqZyKWqloXjwlopa3OeTtwB9BURFKLPVQN+NHXhRnjhF0H83h0RgYzlm6lVb1qdG7zG3si51Kl0njAHfoj5owgNjqWcb3HOVusMWVU2lo9H+FeiXO657bop4uqXlcBtRlTYVSVz5dkcsG/vmNO+nbuu7Al0+/qRWT0FlIWpTBizoijoZ+yKIXs3GxU1emyjSmT0lbn3AfsE5GHgO2qmicivYEOIvK+qmZXTInG+NbmPTmMmZLG96t30TWxJk8P6EDzujEAjO/rPtNPWZRCyqIUAIZ3H874vuOtu8cELDnZWYuI/IF75m4iMAf3J4BWqnqJr4srkpycrIsXL66ow5kQUVDo4t2FG3jxq1WEhwkjL27NoG4JhB03EUtVCXvsfx+OXWNdFvomIIjIElVNPn6/N8syu1S1ABgAvKSqI4D48i7QmIqUsXU/A15fyBMzl9OreW3m3nsOg3s0LjH0R8wZccy+om4fYwKVN6N68kXkGuB63H384F6vx5iAk5tfyMvzVvPmgnXUrBLJq9eeyaXt40s8gy/ep1/UvVN0H7DuHhOwvAn+G4HbgCdVdb2INMG9SJsxAeXndbsZPTmN9bsOcXWXhvzz0jbEVql0wueLCLHRscf06Rf1+cdGx1rom4B10j5+f2B9/OZ07DuczzOzVvDfXzaRUKsKT/Vvz19a1PH69ceP27dx/CZQnKiP35uZu+txr59/DFVtWk61GeMzs9O3M3ZaOrsO5jH07Cbce2ErKlc6tbXyjw95C30T6Lzp6in+v0U0cDVglxYyfm3H/lwembaM2cu20ya+Om/f0JX2DWs4XZYxfsGbRdp2H7frJRH5ARjrm5KMKTtV5eNfN/PUl8vJK3DxYL9WDD27KZHh3gxgMyY0eNPV07nY3TDcnwDsQizG76zfdYjRk1P5ed0eujepxdMD2tM0LsbpsozxO9509bxYbLsA97Vy/+aTaowpg/xCF//5fh0vfb2aqPAwnh7Qnr8nN/rTmHxjjJs3XT19KqIQY8oiLXMfIyelkrFtP33b1eOxK5OoVz3a6bKM8Wulrc55nar+n4jcW8LDCuwBpqvqXp9VZ8wJHD5SyPivV/HW9+uoExPFG9d1pl+STSg3xhulnfEXXVv3RP35TYDbgR7lWpExJ/Hjml2MnpzGpj05XNOtEaMubkONyjaZ3BhvlbY655ue20dP9BwRecwXRRlTkuycIzw5czmfLckksXYV/ju0B2c1q+10WcYEnNK6el4u7YWqereq2pBO43Oqypdp23lk+jL25hzhjt7NuPv8FkRHntpELGOMW2ldPUs8t72AtsAnnvtXF3vMGJ/atu8wD09dxtfLd9C+QQ3eu6kr7c6wiVjGnI7SunreAxCRIUAfVc333H8D+KpCqjMhy+VSPvplE8/MWkGBy8U/L2nDjb0SibCJWMacNm/G8Z+B+wvePZ77MZ59xvjEmqyDjJ6cyq8b9tKreW2e7t+BhNpVnC7LmKDhTfA/A/wuIt967p8LjPNZRSZkHSlw8eZ3a3nlmzVUrhTO83/twF+7NLRF0YwpZ95M4JooIrOA7p5do1R1+8leJyLvAJcBWaqa5NlXC/d3BYl4ZgDbPAAD8MfmbEZNSmXF9gNc1iGeRy5vR1y1KKfLMiYoedthmgdsA/YCLUXkHC9e8y7Q77h9o4B5qtoCmOe5b0LYobwCHpuRQf/XfiQ7J5//XJ/Mq9d2ttA3xoe8WaTtFmA40BD4A/eErZ+A80p7naouEJHE43ZfCfT2bL8HzAdGnkK9Joh8t2onYyansSX7MNf1SGBkv9ZUi7aJWMb4mjd9/MOBrsDPqtpHRFoDJ5zUdRL1VHUbgKpuE5G6J3qiiAwDhgEkJCSU8XDGH+05dIQnvshg8u9baBZXlc9uO4uuiXaJB2MqijfBn6uquSKCiESp6goRaeXrwlR1AjAB3Jde9PXxjO+pKtOXbuXRGRnsP5zP3ec1544+zW0iljEVzJvgzxSRWGAqMFdE9gJby3i8HSIS7znbjweyyvg+JsBsyT7MQ1PS+HblTjo2iuXZge1pXb+602UZE5K8GdXT37M5zjOkswYwu4zHmw7cgHuI6A3AtDK+jwkQhS7lg5828NyclajC2MvackPPRMJtrXxjHOPNGf9Rqvqdt88Vkf/i/iK3johkAo/gDvxPReRmYBPu5R9MkFq14wAjJ6Xy+6ZszmkZx5NXJdGolk3EMsZppxT8p0JVrznBQ+f76pjGP+QVFPLat2t5bf4aYqIiGP/3jlzVqYFNxDLGT/gs+E1oWrJxL6MmpbI66yBXdjqDsZe1pXaMjck3xp9Y8JtycTCvgOdnr+D9nzcSXz2aiUO60qf1CUfrGmMcZMFvTts3K3bw0JR0tu3P5YazErm/bytiouyfljH+yn47TZntOpjHozMymLF0Ky3rxTBpUE86J9R0uixjzElY8JtTpqpM/m0Lj8/MICevkHsvbMlt5zajUoStlW9MILDgN6dk854cxkxJ4/vVu+jSuCbPDGhPi3rVnC7LGHMKLPiNVwpdysQf1/PiV6sIE3j8ynYM6t6YMJuIZUzAseA3J7V8235GTUplaeY+zm9dl8evSuKM2MpOl2WMKSMLfnNCufmFvPLNat78bh01KkfyyjVnclmHeJuIZUyAs+A3JVq0bjejJ6exbtchBnZuyEOXtqFm1UpOl2WMKQcW/OYY+3PzeWbWCj5atIlGtSrzwc3dOLtFnNNlGWPKkQW/OWrOsu2MnZbOzgN5DD27CSMubEmVSvZPxJhgY7/VIURVj+mfL7qftT+XR6YvY1b6dtrEV+c/1yfToWGsc4UaY3zKgj9EjJs/juzcbMb3HY+IoKrcM3sEO3c3J31NC3ILXDzQtxXDzmlKZLhNxDImmFnwhwBVJTs3m5RFKQCM7zueoVPG8MXiukS7mtC9SXWeHtCepnExDldqjKkIFvwhQEQY33c8ACk/v8LEHzKpUXAtMeFhPNY/iX90TbCJWMaEEAv+ECEi3Jw0jk++a0KUNicn7Cd+efCf1K9hE7GMCTXWmRsCDh8p5Kkvl3PFq98TrrXYWekpdkY9yTM/jUZVnS7PGFPB7Iw/yC1cs4vRk9PYuCeHA+Fz+XuPSvz7sh8YMWfEMX3+NhvXmNBhwR+k9uXk8+SXGXy6OJPE2lXo23UFlauGMb7v88f0+cdGx1roGxNiLPiDjKoyK307Y6ctY2/OEW47txn3XNCC6Mg+x4zjLwp/C31jQo8FfxDZvi+Xh6elMzdjB0kNqvPujV1JalDj6OPHh7yFvjGhyYI/CLhcyn9/3cQzX64g3+VizCWtualXEyJsIpYxpgQW/AFu7c6DjJ6Uxi8b9tCzWW2eHtCexrWrOl2WMcaPWfAHqPxCFxMWrCNl3mqiI8J47q8duLpLQ+u+McaclAV/AFq6OZuRk1JZsf0Al7aP55Er2lK3WrTTZRljAoQFfwDJOVLAi1+tYuKP64mrFsWEwV24qF19p8syxgQYR4JfRDYAB4BCoEBVk52oI5AsWLWTMVPSyNx7mOt6JPBgv9ZUj450uixjTABy8oy/j6rucvD4AWHvoSM8PjODyb9toWlcVT699Sy6NanldFnGmABmXT1+SlWZvnQrj83IYN/hfO7q05y7zmtOdGS406UZYwKcU8GvwFciosCbqjrBoTr80pbswzw8NZ1vVmTRsVEsHw5sT+v61Z0uyxgTJJwK/l6qulVE6gJzRWSFqi4o/gQRGQYMA0hISHCixgpX6FL+7+eNPDd7BS6Fhy9ry5CeiYTbWvnGmHLkSPCr6lbPbZaITAG6AQuOe84EYAJAcnJy0K8dvGrHAUZNSuW3Tdmc3aIOT/VvT6NaVZwuyxgThCo8+EWkKhCmqgc82xcBj1V0Hf4ir6CQ175dy2vz11A1KoIXr+7IgM4NbCKWMcZnnDjjrwdM8QRbBPCRqs52oA7HLdm4l1GTUlmddZArOp7B2MvbUicmyumyjDFBrsKDX1XXAR0r+rj+5GBeAS/MWcl7P20gvno07wxJ5rzW9ZwuyxgTImw4ZwX7dkUWD01NZ+u+w1zfozEP9GtNTJT9NRhjKo4lTjkrfrGT4vd3H8zj0RkZTF+6lRZ1Y/j8trPo0tgmYhljKp4FfzkaN38c2bnZR69sparcM3sEe/Y0I31tKw7mFXDPBS24vXczoiJsIpYxxhkW/OVEVcnOzT7mAubDpo5m+q+1qexqSueEqjw7sAMt6lVzuFJjTKiz4C8nxS9gnvLzK7zz43pi86+nWngED1/RlsE9EgmziVjGGD9gwV+ORIRhHR7lk+8aE6UtyQn7hYX3j6RhTZuIZYzxH3ZR1nKSm1/I83NWcOkrC4jQuuyMfI6dlR7jhUVjUA36icfGmABiZ/zl4Jf1exg1OZV1Ow9xMPwb/tYjgn9fNp8Rc0Yc0+dvs3GNMf7Agv807M/N59lZK/hw0SYa1qzMRckrqRKjjO/73DF9/rHRsRb6xhi/YcFfRnMzdvDw1HSyDuRy81+acN9FLalS6bxjxvEXhb+FvjHGn1jwn6KsA7k8Oj2DmWnbaF2/Gm8O7kLHRrFHHz8+5C30jTH+xoLfS6rKZ4szeWJmBrkFLh7o24ph5zQlMty+HzfGBBYLfi9s3H2I0ZPTWLh2N90Sa/H0wPY0i4txuixjjCkTC/5SFBS6ePuH9Yz/ehWRYWE8cVUS13ZLsIlYxpiAFvLBf6JF1dK37GPU5FTSt+znwrb1ePzKJOrXiHawUmOMKR8hHfwlLap295f3sT6zLcs3nkHNKpV4bVBnLk6qb1/SGmOCRsgGf0mLql338WPMT21DpMbzt+QG/POSttSoEulwpcYYU75CNviLT7B6+ee3+OB7F9UK+xIbncOE67rRq3mcwxUaY4xvhGzwgzv8+zYYyee5XQinBvsiPidj9H+oYlfEMsYEsZAdhL5jfy63frCY2z/8jULZw/aoe8mOfJcx395vi6oZY4Ja0J7anmi0jsulfPzrZp6etZxDR/LYG/E+N/RsRMrFq21RNWNMSAjK4C9ptM6IOSOQwnpszezNovV7OKtpbeIbfIsrvCHj+/7LFlUzxoSMoAv+kkbrDJ91L+8t3ErtwvOIidrPcwM7cHVyQ0R62KJqxpiQE3TBf8wlEBel8PpPM6l9ZDg19QL6JtXn0SvaUbd69DHPP/71xhgTzIIu+OF/4T/xh23EFgyikGzevK4zfZPinS7NGGMcF5Sjeor69AtkOwfDv2Jr9O3MynzWRusYYwxBGPxFoZ+yKIWhvc5k1xOvcHePoaQsSmHEnBEW/saYkBd0XT0iQmx0LMO7Dz/6Ra2N1jHGmP8RJ86ARaQfkAKEA2+p6jOlPT85OVkXL158Ssc40Th+Y4wJFSKyRFWTj99f4V09IhIO/Bu4GGgLXCMibX1wnFLvG2NMqHKij78bsEZV16nqEeBj4EoH6jDGmJDkRPA3ADYXu5/p2XcMERkmIotFZPHOnTsrrDhjjAl2TgR/SX0uf/qiQVUnqGqyqibHxdkSycYYU16cCP5MoFGx+w2BrQ7UYYwxIcmJ4P8VaCEiTUSkEvAPYLoDdRhjTEhyajjnJcBLuIdzvqOqT57k+TuBjcftrgPs8kmBzgm2Nll7/F+wtSnY2gOn16bGqvqnvnJHgr88iMjiksanBrJga5O1x/8FW5uCrT3gmzYF3ZINxhhjSmfBb4wxISaQg3+C0wX4QLC1ydrj/4KtTcHWHvBBmwK2j98YY0zZBPIZvzHGmDKw4DfGmBATkMEvIv1EZKWIrBGRUU7X4w0ReUdEskQkvdi+WiIyV0RWe25rFntstKd9K0WkrzNVn5iINBKRb0VkuYgsE5Hhnv2B3KZoEflFRJZ62vSoZ3/AtgncK+KKyO8i8oXnfsC2R0Q2iEiaiPwhIos9+wK2PQAiEisin4vICs/v01k+b5OqBtQP7klfa4GmQCVgKdDW6bq8qPscoDOQXmzfc8Aoz/Yo4FnPdltPu6KAJp72hjvdhuPaEw909mxXA1Z56g7kNgkQ49mOBBYBPQK5TZ467wU+Ar4Ign93G4A6x+0L2PZ46nwPuMWzXQmI9XWbAvGMPyCXdVbVBcCe43ZfifsvHc/tVcX2f6yqeaq6HliDu91+Q1W3qepvnu0DwHLcq6wGcptUVQ967kZ6fpQAbpOINAQuBd4qtjtg23MCAdseEamO+6TwbQBVPaKq2fi4TYEY/F4t6xwg6qnqNnAHKVDXsz+g2igiicCZuM+QA7pNnm6RP4AsYK6qBnqbXgIeBFzF9gVyexT4SkSWiMgwz75Abk9TYCcw0dMd95aIVMXHbQrE4PdqWecAFzBtFJEYYBJwj6ruL+2pJezzuzapaqGqdsK9amw3EUkq5el+3SYRuQzIUtUl3r6khH1+0x6PXqraGfcV/O4UkXNKeW4gtCcCdxfw66p6JnAId9fOiZRLmwIx+INpWecdIhIP4LnN8uwPiDaKSCTu0P9QVSd7dgd0m4p4Pm7PB/oRuG3qBVwhIhtwd4meJyL/R+C2B1Xd6rnNAqbg7uYI2PbgrjHT88kS4HPc/xH4tE2BGPzBtKzzdOAGz/YNwLRi+/8hIlEi0gRoAfziQH0nJCKCu19yuar+q9hDgdymOBGJ9WxXBi4AVhCgbVLV0araUFUTcf+efKOq1xGg7RGRqiJSrWgbuAhIJ0DbA6Cq24HNItLKs+t8IANft8npb7TL+C34JbhHkawF/ul0PV7W/F9gG5CP+3/tm4HawDxgtee2VrHn/9PTvpXAxU7XX0J7/oL7I2Yq8Ifn55IAb1MH4HdPm9KBsZ79AdumYnX25n+jegKyPbj7w5d6fpYV/e4HanuK1dgJWOz5dzcVqOnrNtmSDcYYE2ICsavHGGPMabDgN8aYEGPBb4wxIcaC3xhjQowFvzHGhBgLfmNOk4gkisi1ZXhdrIjcUez+GSLyeflWZ8yfWfAbc/oSgRKDX0QiSnldLHA0+FV1q6r+tVwrM6YEFvwmYInI9SKS6lk//wMRaSwi8zz75olIgud574rIyyKyUETWichfi73Hg5713ZeKyDOefc1EZLZnIbDvRaT1Sd7nGeBszxrxI0RkiIh8JiIzcC8oFuOp5zfPsa4s9rpmntc97/nkkO45VrSITPQ8/3cR6ePZP0REJnvqWy0iz1XIH7YJLk7PWrMf+ynLD9AO98zFOp77tYAZwA2e+zcBUz3b7wKf4T7RaYt7WW9wL/S1EKhS9B6e23lAC892d9xLHZT2Pr3xzIr13B+Ce3Z20ftFANU923VwL6UruD8pFL8+w9H7wH3ARM92a2ATEO1573VADc/9jUAjp/8+7Cewfkr7GGqMPzsP+FxVdwGo6h4ROQsY4Hn8A9wXsygyVVVdQIaI1PPsuwB3uOYUe48YoCfwmXs5IsB90YvS3qckc1W16PoLAjzlWUnShXsZ3dJeC+4lMV7x1LVCRDYCLT2PzVPVfQAikgE05tileo0plQW/CVTCyZejLf543nGvPdF7hAHZ6l6auSQlvU9JDhXbHgTEAV1UNd+zWmZ0Ka892XsXr6EQ+z02p8j6+E2gmgf8TURqg/u6q7i7bf7heXwQ8MNJ3uMr4CYRqVL0Huq+psB6Ebnas09EpONJ3ucA7stPnkgN3Ovi53v66ht78boFnjYgIi2BBNxdW8acNgt+E5BUdRnwJPCdiCwF/gXcDdwoIqnAYGD4Sd5jNu5lbheL+6pb93seGgTc7HnfZZz80p6pQIHnC+IRJTz+IZAs7ouDD8K91DOquhv4UUTSReT5417zGhAuImnAJ8AQVc3DmHJgq3MaY0yIsTN+Y4wJMRb8xhgTYiz4jTEmxFjwG2NMiLHgN8aYEGPBb4wxIcaC3xhjQsz/A1N/NiW0tvXFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r-squared hcy =  0.9985190146227263\n"
     ]
    }
   ],
   "source": [
    "y1_values = []\n",
    "y2_values = []\n",
    "substring = 'STD'                                         #Include files with 'STD'\n",
    "path = os.getcwd()\n",
    "files = glob.glob(os.path.join(path, \"*.csv\"))\n",
    "for file in files:\n",
    "    file_name = file.split(\"\\\\\")[-1]\n",
    "    if re.match(substring, file_name, re.IGNORECASE):\n",
    "        calculate_standards(file_name)\n",
    "\n",
    "plt.scatter(x1_values, y1_values, marker='o', linestyle='None', color='r', label='cys') \n",
    "plt.xlabel('concentration')\n",
    "plt.ylabel('adjusted area')\n",
    "plt.legend()\n",
    "trend1 = np.polyfit(x1_values, y1_values, 1)\n",
    "poly1 = np.poly1d(trend1)\n",
    "new_x1 = np.linspace(x1_values[0], x1_values[-1])\n",
    "new_y1 = poly1(new_x1)\n",
    "plt.plot(new_x1, new_y1)\n",
    "slope1, intercept1, r_value1, p_value1, std_err1 = stats.linregress(x1_values, y1_values)\n",
    "plt.show()\n",
    "plt.clf()\n",
    "\n",
    "int_cys = intercept1\n",
    "slope_cys = slope1\n",
    "print('r-squared cys = ', r_value1**2)\n",
    "\n",
    "plt.scatter(x2_values, y2_values, marker='x', linestyle='None', color='g', label='hcy') \n",
    "plt.xlabel('concentration')\n",
    "plt.ylabel('adjusted area') \n",
    "plt.legend()\n",
    "trend2 = np.polyfit(x2_values, y2_values, 1)\n",
    "poly2 = np.poly1d(trend2)\n",
    "new_x2 = np.linspace(x2_values[0], x2_values[-1])\n",
    "new_y2 = poly2(new_x2)\n",
    "plt.plot(new_x2, new_y2)\n",
    "slope2, intercept2, r_value2, p_value2, std_err2 = stats.linregress(x2_values, y2_values)\n",
    "plt.show()\n",
    "\n",
    "int_hcy = intercept2\n",
    "slope_hcy = slope2\n",
    "print('r-squared hcy = ', r_value2**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c31f63b",
   "metadata": {},
   "source": [
    "# Get Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "761fb7ed",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Get sample data\n",
    "ids = []\n",
    "final_cysteine = []\n",
    "final_homocysteine = []\n",
    "substring = 'STD'                                         #Exclude files with 'STD'\n",
    "path = os.getcwd()\n",
    "files = glob.glob(os.path.join(path, \"*.csv\"))\n",
    "for file in files:\n",
    "    file_name = file.split(\"\\\\\")[-1]\n",
    "    if re.match(substring, file_name, re.IGNORECASE):\n",
    "        pass\n",
    "    else:\n",
    "        calculate_HPLC_results(file_name)\n",
    "        \n",
    "final_data = pd.DataFrame(list(zip(ids, final_cysteine, final_homocysteine)),\n",
    "               columns =['ID', 'cysteine','homocysteine'])\n",
    "final_data.to_excel('final_data.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fca496f8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}