diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml new file mode 100644 index 0000000..b7053c0 --- /dev/null +++ b/.github/workflows/draft-pdf.yml @@ -0,0 +1,24 @@ +name: Draft PDF +on: [push] + +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # This should be the path to the paper within your repo. + paper-path: paper/paper.md + - name: Upload + uses: actions/upload-artifact@v4 + with: + name: paper + # This is the output path where Pandoc will write the compiled + # PDF. Note, this should be the same directory as the input + # paper.md + path: paper/paper.pdf diff --git a/docs/source/index.rst b/docs/source/index.rst index c3d98d2..89ffaeb 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -46,6 +46,8 @@ The documentation features various notebooks that demonstrate the usage and inve notebooks/Ex3_Pipeline_with_larger_example_dataset.ipynb notebooks/Investigation_doublepeak_separation.ipynb notebooks/Investigation_noise_sigma.ipynb + notebooks/Processing_test_1_raw_data.ipynb + notebooks/Create_validation_plot_from_raw_data.ipynb API Reference diff --git a/docs/source/notebooks/Create_validation_plot_from_raw_data.ipynb b/docs/source/notebooks/Create_validation_plot_from_raw_data.ipynb new file mode 100644 index 0000000..5fc93d0 --- /dev/null +++ b/docs/source/notebooks/Create_validation_plot_from_raw_data.ipynb @@ -0,0 +1,274 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook contains the code for creating the validation plot depicted in the documentation chapter denominated as \"Validation of `PeakPerformance`\" from the processed data of test 1 and raw data of tests 2 and 3. \n", + "For the data processing of test 1, the reader is referred to the `Processing test 1 raw data` notebook. \n", + "The raw data files are located within the `PeakPerformance` repository under `/docs/source/notebooks`." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "import arviz as az\n", + "import json\n", + "import numpy as np\n", + "import pandas\n", + "import pymc as pm\n", + "from matplotlib import pyplot as plt\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 1) Preparation of evaluation of synthetic data (test 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "with open('test1_all_data.txt', 'r') as file:\n", + " all_data = json.loads(file.read())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 2) Prepartion of border-line cases normal vs. skew normal (test 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "df = pandas.read_excel(\"test2_summary.xlsx\")\n", + "df_normal = df[(df.loc[:, \"model\"] == \"normal\") & (df.loc[:, \"Unnamed: 0\"].isin([\"area\", \"height\"]))]\n", + "df_normal.reset_index(inplace=True)\n", + "df_skew = df[(df.loc[:, \"model\"] == \"skew_normal\") & (df.loc[:, \"Unnamed: 0\"].isin([\"area\", \"height\"]))]\n", + "df_skew.reset_index(inplace=True)\n", + "df_comparison = pandas.DataFrame()\n", + "df_comparison.loc[:, \"ratio_mean_normal_to_skew\"] = df_normal.loc[:, \"mean\"] / df_skew.loc[:, \"mean\"]\n", + "df_comparison[\"ratio_sd_normal_to_skew\"] = df_normal[\"sd\"] / df_skew[\"sd\"]\n", + "df_comparison[\"parameter\"] = df_normal[\"Unnamed: 0\"]\n", + "df_comparison[\"test_iteration\"] = df_normal[\"test_iteration\"]\n", + "df_comparison_area = df_comparison[df_comparison[\"parameter\"] == \"area\"]\n", + "df_comparison_height = df_comparison[df_comparison[\"parameter\"] == \"height\"]\n", + "comparison_dict = {}\n", + "comparison_dict[\"fraction of mean (normal / skew normal)\"] = [[df_comparison_area[\"ratio_mean_normal_to_skew\"].mean(), df_comparison_height[\"ratio_mean_normal_to_skew\"].mean()], [df_comparison_area[\"ratio_mean_normal_to_skew\"].std(), df_comparison_height[\"ratio_mean_normal_to_skew\"].std()]]\n", + "comparison_dict[\"fraction of standard deviation (normal / skew normal)\"] = [[df_comparison_area[\"ratio_sd_normal_to_skew\"].mean(), df_comparison_height[\"ratio_sd_normal_to_skew\"].mean()], [df_comparison_area[\"ratio_sd_normal_to_skew\"].std(), df_comparison_height[\"ratio_sd_normal_to_skew\"].std()]]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 3) Prepartion of comparison to MultiQuant (test 3)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "df_comparison_total = pandas.read_excel(\"test3_df_comparison.xlsx\")\n", + "df_comparison_single = df_comparison_total[~df_comparison_total[\"PP experiment\"].isin([23, 24])]\n", + "df_comparison_double = df_comparison_total[df_comparison_total[\"PP experiment\"].isin([23, 24])]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 4) Plotting in one graph (for PeakPerformance paper)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAHJCAYAAAArV+dUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACw+0lEQVR4nOzdeXxMZ///8dckSIi1ghDEdttFJEGD2EKCqtr3nbtVd6uotrrZen9V76qiqlq3ShXVKl2Uqlpiv5Wi9tYWIqjS2IKs5/dHfhkiiywzk8nk/Xw8PGRmzrnOZyZnrnnnzHXOZTIMw0BERERERGzGKbcLEBERERHJbxTCRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGbBLC9+/fj8lkwmQy0bNnT1tsUkRERETEbuUohNeuXdscrk0mE15eXhiGkWo5X19fnnjiCQBWrVrF9u3bc7JZEREREZE8LdshfPfu3fz+++8p7jt//jybN29Oc/nx48ebf37ttdeyu1kRERERkTwv2yE8NDQ0S/e3bt0aT09PAHbs2MGhQ4eyu2kRERERkTwtWyH83r17fPXVV+bbNWvWNP+8evVqbt26lXpDTk5069bNfHvRokXZ2bSIiIiISJ6XrRD+zTffcP36dfPtzz77DGdnZwDu3LmTIqA/KCAgwPzzpk2bsrNpEREREZE8L1sh/MEhJ02aNOHxxx+nTZs2aT7+oEaNGpl/Pnr0aIogLyIiIiKSXxTI6gqRkZFs3LjRfLtfv34A9OnTx3z/jh07OHXqFDVq1EixbtmyZVPcvnjxIiVLlsxqCVaTmJjIxYsXKVasGCaTKbfLERERERErMAyDW7duUaFCBZyccmfanCyH8CVLlpCYmAgkjfPu3bs3AD169GD06NHExcUBSUNU3nrrrRTrFi9ePMVtezsSfvHiRSpVqpTbZYiIiIiIDURERFCxYsVc2bbJSOvC3hmoXbu2+dKErVq1IiwszPxY586dWbt2LQCVK1cmPDw8xRHla9eu4e7ubr599OhR6tatm5P6LerGjRuULFmSiIiIVH8wiIiIiIhjuHnzJpUqVeL69euUKFEiV2rI0pHwh68NvnXr1nSHbSRfMzwoKMh835UrV1IsU6FChaxs3uqSn0vx4sUVwkVEREQcXG4OP87SIJj0TrjM7PL79+83/1yvXj27Gg8uIiIiImIrmT4S/vC1watVq0bjxo1TLffbb79x4sQJIOma4fPnz6dYsWJA0pH0ZA8eIRcRERERyU8yHcIfvjb49OnT6dOnT6rl1q1bxxNPPAHcv2b4iBEjSExM5JtvvjEvN3z48ByULSLiGKKjoylatCgAt2/fxs3NLZcrEhERW8j0cJTPPvvM/PNjjz1G165d01wuJCQkxVjv5CEpW7Zs4eLFiwC0aNGChg0bZqNcERH7EB0djclkwmQyER0dndvliIhIHpPpI+Hr16/P1HLOzs5ERkamun/WrFnmn6dPn57ZzYqIiIiIOBybXJ18//79rFu3Dki6nnhgYKAtNisiIiIiYpeyPFlPdvj6+pLFy5GLiIiIiDis3JmnU0REREQkH7PJkXARkfws+mJf+v9zJ+s3XuSrxYF0CvY0PxYXHX//56vDiYl2Zswre1n42Sk+fr8pQwdUT9FWQfdlNqtbRESsR0fCRUSsLL0A/rDERCPDAC4iIo5DR8JFRKwsMwEcYMKbvxK6/IwCuIhIPqAQLiJiZZkJ4ACLl53hk9kK4CIi+YGGo4hIvpFbE+w8aghKsg/e8VcAFxHJJxTCRURySWKiwavTDmAywSezm/LM8JoZLj/9vSM2qkxERKxNIVxEJBdk9STM6e8dYcqMQzaqTkRErE1jwkXEKqKjoylatCgAt2/fxs3NLZcrypq4tQMyfvzeA5cWXD+cONcMutOmKW9mN4BPmej9yLrzo7y+rz2Koz8/kfzKbo+Eb9u2jSeffJIKFSpgMpn49ttvH7nO1q1b8fPzw9XVlWrVqrFgwQLrF5pP5dbYWpG8LicB/LUX69uoShERsTa7DeHR0dE0bNiQefPmZWr5s2fP0qlTJwIDAzlw4ACvvfYaY8aMYdWqVVauVMRx6I8r61IAT0n7m+3pNc85R38NHf352RO7HY7SsWNHOnbsmOnlFyxYQOXKlZk9ezYAderUYd++fcycOZMePXpYqUrJKUt+zWqvbYmAAnh+pH7Efjj678Jen5+91mUv7DaEZ9Xu3bsJDg5OcV9ISAiLFi0iLi6OggULplonJiaGmJgY8+2bN29avU4Rsa6rQyek+1h0XNz95Z5+jbtp9AvJSvSyXE0K4HlHboWGjM5ByNL5B0DBJ5ZZrC6xLoXU/M1hQvjly5cpV65civvKlStHfHw8V69epXz58qnWefvtt5k6daqtSkyXc8vnLNbWyikTLdJOwtbfeG/ldH79Yy+v9H0T/1opzyy7F3vP/PN3039myc//ZcO+dYx+ahztfENSLNvi9IZ0t5OVUARQoteldB/LygfV3oJX6dhzC/Vql+CHr9pQrGjK7cZF32/r2N6edO23leLFC7JhdRAe5QqnWLage976wMtofzMS7v8+igWPx+Sc/u/jUfvavbt3zD9/FxaJa+Ei6S7bvW3FDNvK6/JrAH9U32ap/S0r+xpYbn/L6A8+yL0/+jzLF0m3v3qwb7twrB9dB27l5s04NnwbRM3qxVMsm1/7NrDc/mbJvs1SBxjAcp+l03dk3F89uL8t/rA1Yybu459DajD3ncY4OZlSLJvX9jdLcJgQDmAypfyFGoaR5v3JXn31VcaPH2++ffPmTSpVqmS9AvOQjAL4wz5dv4DNBzakGcDtVUYB/GFd+oRRsmShND/QRDIjPwZwe+bk5JRuf/XgAYbhTzzNV1uX07/tEHq17p9q2RZWrTL70gvgD+vSP4zbt+PTDOBiGaWLu1PE1Y1pQ9+hVLHHUj3+4P5Ws0pt7sXcZdqw/+Dpnjq82+P+lpX+6vlX9vH00LQDeH7lMCHcw8ODy5cvp7jvypUrFChQgNKlS6e5jouLCy4uLrYoL8/JbAAH2LT/J/7VdXy2ArhbwYL8NeTF7JSYI48K4G5uBTjyv84Ed92U6Q80R2ByLogp8F+5XYbDUQBPW27tb5k9YJBRAM+M3OrfMuqv3NwKcP5IN4K7b0r3CLgjyq19LaMA/rC79+7w1vB30wzgj5Jb+9qj+is3twJJE5GN25PuEfD8zGFCeEBAAGvWrElx34YNG/D3909zPLhk7FEBPDEx0fzz052ft4sj4G6uBYj9IXMflo86Av7H6Zv5LoCL9SiA25fM9le9W/XPdgDPTRn1V5f/vJvvAnhuelQAvxtzh4plKnPnXnS6R8BtLSufpY/qr0KXnVYAz4DdhvDbt29z6tQp8+2zZ89y8OBBHnvsMSpXrsyrr75KZGQkS5YsAWDUqFHMmzeP8ePH889//pPdu3ezaNEivvjii9x6CnnaowL4kp//i8lkylNDUB5kqQB++c+7VHK3RoV5n2vhIqzaFJHbZdi9rATw0GWn+ecLNiosn1qzazVAjo6AW1pWQlFGshLAb92O4zH1bWnKSt+WUQCPuvU3k0JfyVQAvxtzJ93H7FVWAvi6DZE8ZR9vN5uy2+uE79u3j0aNGtGoUSMAxo8fT6NGjZg0aRIAly5d4vz58+blq1atyrp16wgLC8PHx4e33nqLuXPn6vKEFpaYmMgna+elexLmw1aGLbdRZZaR1QAe3H2TDauTnEr+yvavIS/iZgffkGU1gD8zbo+NKstbkkPRqk0RjzwpMyMrw5azfPNnmQrgG/f/lO3t5IasBvDOvbfYsLr8J6sBfNqS121YXc5lNYD3HrbdhtXZD7s9Et66dWvziZVpCQ0NTXVfq1at2L9/vxWryt+yE8CXb/6MObkwTi07shPAb96MS3cZcWw5PTqZnQD+zyE1sr09yVhWA/j8795nyuDxGS5nL7ITwI+euGHDCvOX7ATw81fCbVdgDmUngHdoV8GGFdoPuz0SLvYluwG8f9shNqowZ7IbwDd8G2TDKsVRZDeAz32nsY0qzF+yE8CD/TvZqLqcyW4A//HrNjasMv/IbgCfPORtG1aZfdkN4MsXNrdhlfZDIVweKScB3F7GVGYkJwFcJzVZn6NNoZyTAK6TmiwvuwH86ScsN7+DteQkgDf21YBwS8tJAK9ZsbYNK82enATwQoWcbVip/VAIlwwpgN+nAC45pQBuX3ISwJ2c7PvjUwHcviiA36cAfp999yKSqxTA71MAl5xSALcvCuBJFMCtTwH8PgXwlOy7J5FcowB+nwK45JQCuH1x5AAOKIDbEQXw+xTAU8vR1VHOnDnDd999x969e/n7778pVaoUvr6+PPXUU9SsWdNSNUoucOQADiiAi00pgNsXRw7ggAK4HXHkAA4ogOdQtkL48ePHeeWVV1i3bh2GYaS4lOBXX33FxIkT6dChA++++y5169a1WLFiO44cwAEFcLEpBXD74sgBHFAAtyOOHMABBfAcylYI9/b2JiEhgapVq9KmTRvq1q1L8eLFuXnzJseOHWPLli38+OOPbNy4kZiYGEvXLDbgyAEcUAAXm1IAty+OHMABBXA74sgBHFAAz6FshfCOHTvyyiuv0Lx5+td13LlzJ//5z3+yXZjkLkcO4IDFArimdpac0tTO9iWvB/CMZCWAx8YmkPtzyuZ9lgrgcfH2OTGcpQL43v1XaRZsjQrtW7Z6l++//z7DAA7QvHlzvvvuu2wVJfZNUzsn0dTOklOa2tm+ZCWA7/t9jw0ry7msBvD+/9xpw+ryn6wG8PdWTrdhdTmX1QDesWf+/Cy12J/4P/30E9OmTbNUc2KnsnNVgbxCUzuLLWlqZ/uS1QD+zoq3bFhdzmQngK/feNGGFeYv2Qngv/6x14YV5kx2Ani92iVsWKH9sFgIX79+PVOnTrVUc2KHNLVzEk3tLDmlqZ3tS3YCuF/NxjasMPuyG8C/Whxowyrzj+wG8Ff6vmnDKrMvuwH8h6/y52ep4wx2E6vS1M5JdFKT5JSuq2tfshvAX+z1mg2rzJ6cBPBOwZ42rDR/yEkA96/V1IaVZk9OAnixovnzDIQcXSdc8gdHnthCAdw+fPvWBiDtiS3uxd4zL7dmxiZcC7kC6X+gtbB9+ZmmAG5fchLACxaw79CgAG5fFMDvUwC/z2IhvH///vj4+FiqObETCuBJFMCtTzPL3acAbn0K4EkUwK1PAfw+BfCULBbCGzduTOPGeWOMnGSOAngSBXDrUwC/TwHc+hw5gAMK4HZEAfw+BfDUspWUxo0bR2RkZIbLREZGMm7cuGwVJbnPkQM4oABuZxw5gIOmdrY3jhzAAQVwO+LIARxQAM+hbB0JnzNnDvPmzaNVq1a0bds21YyZmzZtYtu2bSQmJvL++3nnEnVynyMHcEAB3M5kJoADzFg+lQtXz+epAA6a2tneOHIABxTA7YgjB3BAATyHshXC165dy0svvcTmzZvZsiX1BdYNw6BOnTrMnDkzxwVK7nDkAA4ogNuZzARwgIi/zjFl6Iw8FcBBUzvbG0cO4IACuB1x5AAOKIDnULanre/YsSM///wz33zzDfv37ycqKoqSJUvi6+tL165dCQnJeNpzsW+OHMABiwVwTe1sGY8agpLstQFTNbVzPpza2dIcOYBnJCsBPDHRsGFljstSATwxMdEa5eWYpQL4H6dvUi8fHuvK0YmZ7du3p3379paqRfKIrJ7U1CIPXQgzO1cV+G6tDQvMZ+7G3GHG8vuTgNXwrJXusskfaG2aBNiiNIvIzklNN27asMB8JisB/I8LJ+z6cpgPy2oAH/PKXj5ebMMC85msBvBP1s6jZd288w1gVgN4cNdNRF6yYYF2Im8dvpRcp6mdk2hqZ+tLPgkz4q9zj1xWUztLTmU1gE/97FUbVpcz2QngCz87ZcMK85fsBPAN+9bZsMKcyU4AL17ccb5xygqFcMk0Te2cRFM7W9+DV0F5bcDUDJfV1M6SU9kJ4JXLVrFdgTmQ3QD+8ft5a2xyXpHdAD76qbxxtbnsBvANq4NsWKX9UAiXTNHUzkl0UpP1PXwZwswMQclLJzXpurr2JbsBfNLg/7NhldmTkwA+dEB1G1aaP+QkgLfztf/z7HISwD3KFbZhpfZDIVweyZEntlAAty+a2OI+BXDry0kAL+xSxIaVZp0CuH1RAL9PAfw+i58yt3z5crp160bhwvn3RXUkCuBJFMCtL70A7lrIlW+m/ZRiWQVwySkF8CQK4NanAH6fAnhKFjsSnjyD5sCBA/Hw8GDYsGFpXkNc8g4F8CQK4NanI+D3KYBbnyMHcEAB3I4ogN+nAJ6aRUL41q1bGTVqFAC///47r732GsePH6d9+/Z4eXkxefJk/v77b0tsSmzEkQM4oABuZxw5gIOmdrY3jhzAAQVwO+LIARxQAM8hi4Tw/fv307JlSwD+8Y9/EBwcTJMmTShatCiXL19m8eLFBAdrhom8xJEDOKAAbmccOYCDpna2N44cwAEFcDviyAEcUADPIYuE8HHjxlGgQNLwcm9vb/z9/Vm7di0TJ07kwoULHDp0iP3791tiU2IjjhzAAQVwO+PIARw0tbO9ceQADiiA2xFHDuCAAngOWezEzOHDhwNQrVo1/vOf/xASEoLJdH+q5rNnz1pqU2IDjhzAAYsFcE3tbBma2llTO9uSIwfwjCiA254jB3DAYgH88p93qZQP+zaLnZjZsGFDgoKCaN26NTVr1kwRwAG8vLwstSmxAUsF8D8unLBGeVaTnasKiPVk56SmvCQ7JzWJ9WQlgEdevWDDynIuqwF8+ntHbFRZ/pTVAL4ybLmNKrOMrAbw4O75s2+zWAgPDw/n3Xff5caNG/Tu3Zv69eszceJEdu3ahWFk72jh/PnzqVq1Kq6urvj5+bF9+/Z0lw0LC8NkMqX6d+JE3gqB9k5TOyfR1M7Wp6md78vvUzvbQlYD+KTFL9uwupzJTgCfMuOQjarLf7ITwJdv/sxG1eVcdgL4zZtxNqzQflh0sh5fX18mT57Mvn372LBhA9WqVWPGjBm8//77WW7ryy+/ZOzYsbz++uscOHCAwMBAOnbsyPnz5zNc7/fff+fSpUvmf//4xz+y+3TkIZraOYmmdrY+Te18n6Z2tr7sBPAirm42rDD7shvAp0z0tlGF+Ut2A3j/tkNsVGHOZDeAb/g2f/ZtFp8xc/ny5dy7d48KFSrw9NNP8/333zN+/PgstzNr1ixGjBjByJEjqVOnDrNnz6ZSpUp89NFHGa5XtmxZPDw8zP+cndMfixkTE8PNmzdT/JO0aWrnJBpTaX26ru59OqnJ+rIbwKcNfceGVWZPTgL4ay/Wt1GV+UdOAniv1v1tVGX25SSA16xe3IaV2g+Lh/CBAwdSrlw5hg0bRlhYWLbaiI2N5ddff011WcPg4GB27dqV4bqNGjWifPnyBAUFPXKyoLfffpsSJUqY/1WqVClb9To6R57YQgHcviiA36cAbn05CeClij1mw0qzTgHcviiA36cAfp/FQ/iDk/W0a9cuW5P1XL16lYSEBMqVK5fi/nLlynH58uU01ylfvjyffPIJq1atYvXq1dSqVYugoCC2bduW7nZeffVVbty4Yf4XERGR6RrzCwXwJArg1qcAfp8CuPUpgN+nAG5dCuD3KYCnZLFLFCZLnqwnMjKSEydOmCfrWbt2Lfv27ctSWw9fYcUwjFT3JatVqxa1atUy3w4ICCAiIoKZM2eaJxJ6mIuLCy4uLlmqKT9RAE+iAG59CuD3KYBbnyMHcEAB3I4ogN+nAJ6axY+EW2KyHnd3d5ydnVMd9b5y5Uqqo+MZefzxxzl58mSml5f7HDmAAwrgdsaRAzhoamd748gBHFAAtyOOHMABBfAcsngIr1atGmvXruXUqVNMnDiRMmXKULJkySxN1lOoUCH8/Pz4+eefU9z/888/06xZs0y3c+DAAcqXL5/p5eU+Rw7ggAK4nXHkAA6a2tneOHIABxTA7YgjB3BAATyHsj0cZe7cuTRq1IjAwEDOnDnDhQsXaNmyJd9++22ay2d1sp7x48czaNAg/P39CQgI4JNPPuH8+fOMGjUKSBrPHRkZyZIlSwCYPXs2VapUoV69esTGxrJ06VJWrVrFqlWrsvsU8zVHDuCAAridceQADpra2d44cgAHFMDtiCMHcEABPIeyHcLHjh3L2LFjCQwM5IMPPmDu3LkkJCRYrLA+ffpw7do1pk2bxqVLl6hfvz7r1q0zh/lLly6luGZ4bGwsEyZMIDIyksKFC1OvXj3Wrl1Lp06dLFZTfuLIARxQALczjhzAQVM72xtHDuCPogBuW44cwAGLBfBbt+N4LB/2bRYfjmJJo0ePJjw8nJiYGH799dcUJ1iGhoamuATiyy+/zKlTp7h79y5///0327dvVwDPAUsFcE3tLDmhqZ3vy89TO9tKVgJ41K3MX/HLXmQlgIcuO22jqvKvrATwjft/slFVlpHVAN65d8aXlHZUdh3Cxf5oauf7NLWzdWlq5/vy+9TOtpDVAD4p9BUbVpdzWQ3gz4zbY6PK8qesBvD532V95vHckp0AfvTEDRtWaD9ydInCw4cP88knn3D48GEAFi5ciGEY5seffvrpnFUndkVTO9+nqZ2tS1M736epna0vOwH8zr1oG1aYM9kJ4P8cUsNG1eU/2Qngwf5545v97AbwH79uY8Mq7UeOQvjmzZvZvHmz+XbySZPJFMIdR46mdr6StevD25omtrAvOb6u7ukNNqo0e3RdXfuS3QA+bdh/4MYxG1aaPdkN4HPfaWyjCvOX7Abwp594Ds5utFGV2ZOTAN7YNx8OCCcHIbxy5crpTpwjjiXHE1tcsWGxWaQAbl80scV9CuDWl5MA7ule0e5DeE4CuJOTPt8tLScB3MnJvkcPK4BnT7ZDeHh4uAXLEHvlyDPLKYDbFwXw+xTArS/HAdzOKYDbFwXwJArgKdn3b1ZylQL4fQrg1qUAfp8CuPUpgN+nAG59CuBJFMBTs+/fruQaRw7ggAK4nXHkAA6a2tneOHIABxTA7YgjB3BAATyH7P83LLnCkQM4oABuZxw5gIOmdrY3jhzAAQVwO+LIARxQAM+hvPFbFptz5AAOKIDbGUcO4KCpne2NIwdwQAHcjjhyAAcUwHMo7/ymxaYcOYADCuB2xpEDOFh2amfJOUcO4I+iAG5bjhzAAYsF8NjYBGuUZ/fy1m9bbMZSAVxTO0tOaWrnJPl5amdbyUoAvxtzx4aVWUZWAvi6DZE2rCx/ykoA3/d73prBNKsBvP8/d9qwOvuhEC5Zoqmd79PUztanqZ2T5PepnW0hqwF82pLXbVhdzmU1gPcett2G1eU/WQ3g76x4y4bV5Ux2Avj6jRdtWKH9UAiXTNPUzvdpamfr09TOSTS1s/VlJ4CfvxJuuwJzKDsBvEO7CjasMH/JTgD3q5k3ZjDNbgD/anGgDau0Hwrhkik5mto5D9DUzvYlR1M72zldV9e+ZDeATx7ytg2rzL7sBvDlC5vbsMr8I7sB/MVer9mwyuzJSQDvFOxpw0rth0K4PJImtrhPJzVZnyNfV1cB3L7kJIDXrFjbhpVmT04CeKFCzjasNH/ISQAvWKCgDSvNOgXw7LHvTyzJdQrg9ymAW58CeBIFcOtTAL9PAdz6FMCTKICnZN+fWpKrFMDvUwC3PgXwJArg1qcAfp8CuPUpgCdRAE/Nvj+5JNcogN+nAG59jhzAQVM72xtHDuCAArgdceQADiiA55D9f3pJrnDkAA4ogNsZRw7goKmd7Y0jB3BAAdyOOHIABxTAcyhvfIKJzTlyAAcUwO2MIwdw0NTO9saRAzigAG5HHDmAAwrgOZR3PsXEphw5gAMK4HbGkQM4aGpne+PIARxQALcjjhzAAYsF8MREwxrl2b289UkmNmOpAK6pnSWnNLVzkvw8tbOtZCWAx8XH2bAyy8hKAN+7/6oNK3Nclgrgf1w4YY3yrCarAXzMK3ttWJ39UAiXLNHUzvdpamfr09TOSfL71M62kNUA/t7K6TasLueyGsA79txiw+ryn6wG8KmfvWrD6nImOwF84WenbFih/VAIl0zT1M73aWpn69PUzkk0tbP1ZSeA//pH3jlyl50AXq92CRtWmL9kJ4BXLlvFdgXmQHYD+MfvN7VhlfZDIVwyRVM736epna1PUzsn0UlN1pfdAP5K3zdtWGX2ZTeA//BVGxtWmX9kN4BPGvx/Nqwye3ISwIcOqG7DSu2HQrg8kia2uE8nNVmfI19XVwHcvuQkgPvXsv8jdzkJ4MWK2vd7KS/KSQAv7FLEhpVmnQJ49iiES4YUwO9TALc+BfAkCuDWpwB+nwK49SmAJ1EAT0khXNKlAH6fArj1KYAnUQC3PgXw+xTArU8BPIkCeGoK4ZImBfD7FMCtz5EDOGhqZ3vjyAEcUAC3I44cwAEF8BxSCJc0OXIABxTA7YwjB3DQ1M72xpEDOKAAbkccOYADCuA5pBAuaXLkAA4ogNsZRw7goKmd7Y0jB3BAAdyOOHIABxTAc8iuQ/j8+fOpWrUqrq6u+Pn5sX17xhOjbN26FT8/P1xdXalWrRoLFiywUaWOx5EDOKAAbmccOYCDpna2N44cwAEFcDviyAEcUADPIbsN4V9++SVjx47l9ddf58CBAwQGBtKxY0fOnz+f5vJnz56lU6dOBAYGcuDAAV577TXGjBnDqlWrbFy5Y7BUALfXqZ0tFcA1tbNlaGpnTe1sD7ISwBMTE21YmWVkJYD/cfqmDStzXJYK4JFXL1ijPKvJagCf/t4RG1VmX+w2hM+aNYsRI0YwcuRI6tSpw+zZs6lUqRIfffRRmssvWLCAypUrM3v2bOrUqcPIkSMZPnw4M2fOtHHljk1TO9+nqZ2tT1M7J8nvUzvbQlYD+Cdr59mwupzLagAP7rrJhtXlP1kN4JMWv2zD6nImOwF8yoxDNqrOvhTI7QLSEhsby6+//srEiRNT3B8cHMyuXbvSXGf37t0EBwenuC8kJIRFixYRFxdHwYKpO5yYmBhiYmLMt2/cuAHAzZu2PQJgxMdarK070bcs0o7TvehU992Nvcs7X0zjwl/nmdh/ChXdK3EnjeUA4hLi+eCbmRw8tZ9bPr4WqQnAdMdCR9bTeJk2bLrIkNG7aNe6PAtmNeFeTCL3YtI+2rX/4DV6DN5G7ZrFbb6/5JSl9jdL7WuQ9v524NQ+3v/6P/jU8OXZLmOJi48lLp3aT108yYzlU6hYpjK3YmPSXCarLLavQar9LTY2gRHP/4+NYZf4bH4zWgSU5eattLeXmGjw8qT9hC4/w5wZfnlqf7PHvg1S728P9lfjer5MXa/66fZtiYmJhG74hE37N3CrRlWL1WTN/e3B/mr5ohYYBunub6fO3qRrv60UK1YwT+1rYJ/7W1p924P91YQ+r2MYRrr728VrF5m+7E0Ku7hZrG8D632WPtxfde9SOd19DWDmB8eY8f5RJo6rZ/P9LXl7hpGLw/wMOxQZGWkAxs6dO1Pc/3//939GzZo101znH//4h/F///d/Ke7buXOnARgXL15Mc53JkycbgP7pn/7pn/7pn/7pn/7lw38RERGWCa/ZYJdHwpOZTCnH7RqGkeq+Ry2f1v3JXn31VcaPH2++nZiYyN9//03p0qUz3I5k3s2bN6lUqRIREREUL148t8sRB6f9TWxF+5rYkvY3yzMMg1u3blGhQoVcq8EuQ7i7uzvOzs5cvnw5xf1XrlyhXLlyaa7j4eGR5vIFChSgdOnSaa7j4uKCi4tLivtKliyZ/cIlXcWLF1fHITaj/U1sRfua2JL2N8sqUaJErm7fLk/MLFSoEH5+fvz8888p7v/5559p1qxZmusEBASkWn7Dhg34+/unOR5cRERERCS32GUIBxg/fjz//e9/+fTTTzl+/Djjxo3j/PnzjBo1CkgaSjJ48GDz8qNGjeLcuXOMHz+e48eP8+mnn7Jo0SImTJiQW09BRERERCRNdjkcBaBPnz5cu3aNadOmcenSJerXr8+6devw8vIC4NKlSymuGV61alXWrVvHuHHj+PDDD6lQoQJz586lR48eufUUhKQhP5MnT0417EfEGrS/ia1oXxNb0v7mmEyGkZvXZhERERERyX/sdjiKiIiIiIijUggXEREREbExhXARERERERtTCBcRh9a6dWtMJhNTpkyxeNthYWGYTCZN7iUiZsl9Qnh4eLbbGDp0aJb7reysI7lLIVxEbCo0NBSTyUTr1q1zuxQREbsUHBzMCy+8wOOPP27RdqdMmYLJZGLo0KEWbVeyx24vUSgiYi8SEhJwcnLSEW8BIC4uTpPAiVX179+f/v3753YZYmU6Ei7pSv5KbcaMGdSvX58iRYowePBgIiMjeeqppyhSpAi+vr4cPnzYvM6lS5cYPnw4Xl5eFC1aFH9/f1avXm1+fOvWrTRu3JhSpUpRsGBBypYtS//+/fnrr79SbXf+/Pn4+flRrFgx2rZtS2RkpE2fv+TM1KlTqVKlCi4uLpQuXZrmzZvTrl07hg0bBiTtCw8O5fjggw+oU6cOxYoVo2DBglSpUoXXXnuNxMRE4P7QjypVqvDOO+9Qo0YNSpUqxciRI4mPjzdv991336VSpUo89thjTJw4kYevwpqVffDDDz+kVq1aFCpUiBs3bnDmzBnatm2Lm5sbjRs35uDBg1Z+FcUWhg0bhpeXF4ULF8bV1ZWGDRuybNky4P6Rw169ejFo0CCKFi3KSy+9BCTN4tyyZUtKlSpFuXLl6NWrF2fPnjW3+/rrr1O9enXc3NwoVKgQtWrV4v3338+V5yi2t3nzZho0aECRIkVSfYadPHmSnj174unpSYkSJQgMDGTr1q3mxx8eWhIfH89LL71E2bJl8fDwYM6cOVSpUgWTyURoaGiK7f7999/07t0bNzc3qlatyqpVq8xtTp06FYDPPvvM3J9KLjJE0gEYgFGiRAlj6NChRuHChQ3AKFmypNGxY0ejbt26BmC0bNnSMAzDuH37tuHl5WUARmBgoDFgwACjVKlSBmD88MMPhmEYxvLly43WrVsbI0aMMJ555hmjevXqBmD07Nkz3e2WKVPGAIyhQ4fmyusgWbdp0yYDMMqUKWM888wzRv/+/Y3atWsboaGhRvv27Q3A8PT0NF544QXjhRdeMAzDMMaOHWs8+eSTxqhRo4whQ4YYJUuWNABj3rx5hmEYxpYtW8z7Ro0aNYyBAwcazs7OBmCEhoYahpG0fwGGk5OT0a9fP6NJkyaGk5OTARiTJ082L5PZfdDFxcXo37+/0adPH+P69evmfb5WrVrG0KFDjaJFi5qXlbyrZcuWRp8+fYzRo0cb3bt3N5ycnAxnZ2fj4MGDxuTJkw3AMJlMRoMGDYxnnnnGmD9/vvH9998bJpPJKFasmNGnTx/zfl21alXj1q1bhmEYRu/evY3u3bsbzz77rNG/f3/D1dU1RX8ojim5TyhZsmSan2ERERFGyZIlDZPJZHTo0MHo06ePUaRIEaNgwYLGgQMHDMMwjCFDhqTot6ZPn24ARqFChYzBgwcbtWvXNvdtixcvTrGOyWQyOnXqZAQEBBiAUaxYMePmzZvGsmXLjKZNmxqAUadOHeOFF14wpk6dmguvkCTTkXB5pGnTprF48WKCg4MB8PT0ZN26dcyYMQOAX3/9FYDVq1dz7tw5SpQoga+vL+7u7tStWxeAefPmAdC3b1/efPNN89Gh+vXrA0lHlB720UcfsXjxYiZOnAjA3r17rftExWJiYmIAqFy5Mt27d2f69OkcP36cgQMHmr9irVGjBrNnz2b27NkAvP322wwYMAAPDw9KlSpF9erVAdiwYUOKtp2cnNiyZQuff/45nTp1Au7vG4sXLwbgmWeeYfny5Wzfvp3SpUunWD8r++CcOXNYtmwZK1as4OjRoxw7dsy8/cWLFzNt2jRLvFySy1auXEmbNm0oWbIklSpVokyZMiQkJLB582bzMpUrV2bfvn0sWLCAZ599ltmzZ2MYBvXr18fDw4O6detSqlQpzp49y9q1awFYuHAhTz75JO7u7pQtW5ZKlSoBqfdpcUzz589P8zNs8eLFXL9+nYoVK1KrVi08PDyoXr06cXFxfPzxx2m2ldy3TZ48mc8++4ywsDCcnNKOcO3bt2ft2rX88MMPANy6dYvff/+d/v3706FDBwCaNGnC7NmzmTRpkkWfs2SNxoTLIyUH6ZIlSwJQp04dAIoVKwZAdHQ0AOfPnwfgxo0bzJkzJ0UbJ0+eBGDMmDHmQP6gGzduEBMTk2JKXn9/fwBKlSoFwO3bty3yfMT6QkJCeP755/n8888JCQkBwMvLiyVLlqS5fFxcHM2aNePAgQOpHrty5UqK2x4eHlSsWBFIvW9EREQAULt2bQAKFSpE9erVUww1yco++ODJo8ltlyhRgvLlywP33xuSd50+fRp/f3+uX7+e6rErV66Y94eAgAAKFSpkfiy5v9u9eze7d+9Osd7JkyeJiorCx8fHvNzD7YrjS+8zLHmfiIiISPez8mEXLlwAMB80KFeuHO7u7ly+fPmR231w22JfdCRcHqlAgZR/qzk7O6e5XHIwqly5Mvfu3cMwDAzDIDY2li1btgCwdOlSAP7zn/8QHx/PF198YV7feGjsbvKJTzoZLu+Jj49nzpw5REVFce7cOSZPnsy5c+eYOnWqef9JHusNcOzYMXMA3717N4mJiTzzzDNA+vsFpN43ko80Hj9+HIDY2FhOnz6dYpms7IMPBvLktq9fv86lS5fMdUve9sMPP3D9+nUqVarEn3/+SWJiovlAw4P7w4P7Atzv715++WVzX2cYBhcvXuTFF19k27ZtnD9/HldXV86cOYNhGOY/SB/ez8QxpfcZlrzvPP744yQmJpr3nejoaHP/9LDk/ufEiRNA0h9yV69ezdJ2gTT7X8k9OhIuFtO9e3def/11zp8/T+PGjWnWrBlXrlxhx44djB49milTpuDh4cH169cJDQ3lxIkT5q/LxLHs2rWLQYMG0axZM9zd3dmzZw+QdGSmcuXKAOzbt49nn32WatWq0b9/f5ydnUlISOCNN97A3d2d77//PsvbHTp0KD///DOffPIJN27c4PTp06k+qLK7DzZt2pQ6depw/Phx2rRpQ0BAACtXrsxyjWJfPDw8gKSTyseNG0dkZCSnTp165HrPP/88YWFhvPfee/z++++ULVuWkydPsmPHDk6ePGlu9969e7z44ovExcWlGN4i+dfQoUN57733+N///kfz5s1p2LAhkZGRbN26lTlz5qR5+cBhw4bx6quvMmXKFI4fP86ePXuyFaST+99169bx3HPP0ahRI0aMGJHTpyTZpCPhYjHFihVj9+7dDB06lJs3bxIaGsrevXsJDAw0j0NbvHgx3t7enD59msOHDzN58uRcrlqsoWLFitStW5etW7eycOFCwsPD6dq1K++//z4tW7Zk8ODBFCxYkAULFrBs2TI8PT35+OOP8fT0ZNeuXcTHxzNq1Kgsb7d///68/fbbeHh4sG7dOlq0aEHz5s1TLJPdfdDZ2Znvv/+eVq1acf78eX777Tftvw6gV69ejB49mqJFi7J+/XoCAwNp1qzZI9fr3r0769ato1mzZmzfvp0vvviCK1eu8Nxzz+Hu7k7Tpk156623KFOmDBs3bqRChQr06tXLBs9I7J2Xlxd79uyhR48enDt3jsWLF3P48GG6dOmS7nXBJ0yYwIQJEyhatCg//vgjTz/9NOXKlQNSf1udkV69evHEE09w9+5dPvzwQ9asWWOR5yTZYzL0vZiIiIiI3bp37x4mk8k8LCoiIgIvLy8Mw2D79u20aNEilyuU7NBwFBERERE7Fh4eTocOHejRowdFihThiy++wDAM89BPyZt0JFxERETEjl2+fJmePXty5MgRYmJi8PLyomvXrkycONF85TLJexTCRURERERsTCdmioiIiIjYmEK4iFhVaGgoJpOJKlWqWKS9KVOmYDKZ0ryMV1aZTCZMJhPh4eE5bktERCQrFMJFJE95/PHHeeGFFwgODs70OkOHDsVkMjFlypQU97/wwgu88MILFC9e3MJVioiIZExXRxGRPKVDhw7m687n1OzZsy3SjoiISFbpSLiIWNT+/ftp3LgxRYoUISgoKNVQj5MnT9KzZ088PT0pUaIEgYGBbN26FYCtW7diMpnw8PAgISEBgJs3b1K4cGFMJhO///57quEox48fp3nz5ri7u1OwYEEee+wxnnzySfN09a1bt+azzz4DYOrUqZhMJlq3bg2kHo5y/fp1XnjhBapXr46bmxv169dnzpw55lqSh9a0aNGCl156icqVK+Pu7s5rr71mxVdUREQckY6Ei4jF3Lp1iw4dOvDXX3/h5+dH+fLlmT59uvnxCxcu0KRJE27cuEFISAglSpRgzZo1tG/fnl9++YVWrVpRo0YNTp06xcaNGwkJCWH16tXcu3eP5s2bU6tWrVTbvHbtGomJiTz55JMUKVKEnTt38sMPP3D58mX27t1Lz549uXLlCsePH6dp06Y8/vjj1KhRI1U7hmHQrVs3wsLC+Mc//kHfvn1ZvXo1Y8eO5e+//2bq1KnmZXfu3Mnt27dp2rQpX3/9NW+//TYhISG0atXKOi+siIg4HB0JFxGLWbNmDX/99RelSpVix44dLF26lGeffdb8+OLFi7l+/ToVK1akVq1aeHh4UL16deLi4vj4448BzEe4ly9fnuL/4cOHp7nNFi1aMGfOHOrUqUORIkXw9vYGYN++fURFRfHcc8/RpEkTIGkoy+zZs3nuuedStfPrr78SFhaGyWRiy5YtLFq0iPnz5wPw/vvv8+DVXEuVKsXOnTtZuXIlDRo0AGDv3r3Zft1ERCT/0ZFwEbGYiIgIAKpUqYKrqysAdevWNT9+/vx583Jz5sxJse7JkyeBpBA+efJkvvnmG6ZNm8bmzZspWrQovXv3TnOb7733HhMmTEjzsStXrlCqVKlM1X727FkASpQogaenZ4rab926xdWrV83L1qlTBzc3NwBz+7dv387UdkREREBHwkXEgipVqgQkTbF87949AI4dO2Z+vGLFikDSFU4SExMxDAPDMIiOjmbp0qUAeHp6EhwczK1btxgyZAgJCQn07t2bokWLprnN5PWef/55YmJi2L17t/mx5KPXzs7OACQmJqZbe9WqVQG4ceMGly5dApLGmwMULVoUd3d387IFCxY0/2wymTJ+UURERNKgEC4iFtO5c2fc3d2JiooiMDCQgQMHmod0QNJR7mLFivG///2P5s2b8+yzz9KlSxfKly/P+vXrzcslDz1JPmEzvaEoAB4eHgB8//33jB49mn79+qVapnLlygAsWbKEMWPG8M0336Raxs/Pj8DAQAzDoHXr1owcOdI8lGbs2LEK2yIiYlEK4SJiMcWLF2fdunX4+flx5MgRIiIieOWVV8yPe3l5sWfPHnr06MG5c+dYvHgxhw8fpkuXLjz++OPm5bp06WI+8lyrVi2aN2+e7jZnz55Ns2bN+PPPP9m1axfTpk1LtczTTz9NixYtuHz5Mh988IE53D/IZDLx7bff8q9//YvY2FiWL19OhQoVeO+995g8eXJOXhYREZFUTMaDZxuJiIiIiIjV6Ui4iIiIiIiNKYSLiIiIiNiYQriIiIiIiI3pOuGS5yUkJBATE5PbZYiISD7m4uJivhyqSGYohEuedvPmTU6dOoXOLxYRkdxkMpmoUaMGxYsXz+1SJI/Q1VEkz0pISOC3336jWLFilC9fHicnja4SERHbS0xM5NKlS9y6dYvatWtTpEiR3C5J8gAdCZc8KyYmBsMwKF++fLqzKYqIiNhC+fLluXnzJl999RU1atSgWbNmOjgkGdLeIXmeOjkREcltyZ9FiYmJ7Ny5k4MHD+ZuQWL3lF5ERERELKREiRKYTCYuXryY26WInVMIF3FQYWFh+Pv722x74eHh5qnmM2LruiRvyK/7RZUqVThy5IjNtte6dWt++OGHRy5n67ocjbOzM3Fxcbldhtg5hXAROxMfH5/bJYhIBvQeFRFLUAgXsSCTycQ777xD06ZNqVq1KosXLzY/tm/fPgICAvD29qZJkybs3LkTuH8Eedq0aQQGBvLBBx/QunVrXnrpJVq2bEmlSpV49913WbFiBc2aNcPLy4sVK1aY2x04cCD+/v54e3vTuXNnrly58sg6hw4dyqhRowgKCsLLy4sXXniBLVu20LJlS6pUqcKsWbMeWTfAhx9+SI0aNQgMDOS///1vim389NNPtGjRAj8/P5o2bcq2bduy/bqK47h79y59+vShbt26NGzYkODg4FTL3Lhxg/bt2/PWW28B6e9L/fr144svvgBg7ty5uLi4EB0dDUBgYCDbt29P1fbQoUMZPXo07dq1o2bNmnTv3p3Y2FgAbt++zfDhw6lfvz7169dn6tSp5vVat27N66+/TlBQECEhIYSGhhIcHEzfvn2pXbs2bdu25ejRozzxxBPUrFmTvn37kpiYCMDy5ctp2rQpjRo1wsfHh3Xr1j3ydQoLC8PHx4dRo0bRoEEDfH19OXLkiPm1a9++Pbdv335k3ceOHaNp06b4+voyYMAA7t27Z37s8uXL9O7dmyZNmuDt7c2kSZMeWZeIWJAhkkdFR0cb+/btM6Kjo3O7FDPAmD17tmEYhnHs2DGjaNGiRlxcnBETE2NUqlTJWL9+vWEYhrF9+3bDw8PDuH37tnH27FkDMJYtW2Zup1WrVkbv3r2NhIQEIzIy0nB1dTVef/11wzAMY8+ePUb58uXNy/7111/mn99++23jX//6l2EYhrFlyxbDz88vzTqHDBliNG/e3Lh3754RHR1tlClTxhg2bJiRkJBgXLhwwXBzczNu3bqVYd2//fabUb58eePy5cuGYRjGs88+a5QuXdowDMM4ffq0ERAQYNy4ccMwDMM4efKkUaFCBSM2NjZVXR07djT27t2bg1dd8pLVq1cb7du3N9++du2aYRj399dz584Zvr6+xueff24YRsb70n//+19j2LBhhmEYRpcuXYyAgADjxx9/NG7dumWUKlXKiI2NTbX9IUOGGAEBAcadO3eM+Ph4o1mzZsby5csNwzCMl19+2RgwYICRkJBg3L592/Dx8TG++uorwzCS3pOdOnUyt7l48WKjZMmSRkREhGEYhvHEE08Y3t7exvXr1424uDijQYMGxk8//WQYhmFcvXrVSExMNAzDMM6ePWuUL1/e3I6Xl5dx+PDhVHVu2bLFKFCggHHgwAHDMAxj9OjRhqenp3l7HTt2ND7++ONH1u3r62uEhoYahmEYu3fvNpycnIw1a9YYhmEYwcHBxtatWw3DMIy4uDgjJCTEWL16daq6PvroI+PNN9981K8230v+TPr666+NmTNnml9LkfToEoUiFjZgwAAA6tSpQ4ECBbh8+TJRUVEUKlSIkJAQAFq0aEHZsmU5dOgQ5cuXx9XVlX79+qVop1evXjg5OVGhQgXc3d3p2rUrAH5+fly6dIl79+7h6urKsmXL+Pzzz4mJieHu3bt4eHhkqs6uXbvi4uICQK1atejUqRNOTk54enpSqlQpLly4QFxcXLp17927lyeeeIJy5coB8PTTT/PVV18BsH79ek6dOkXLli1TbDMiIiJVHZk5KiiOo2HDhpw4cYLRo0fTqlUrOnXqZH7s0qVLtGrVisWLF9O6dWsg432pffv2TJ06lYSEBI4fP8706dPZuHEjCQkJBAQEULBgwTRr6N69O4ULFwagSZMmnD59GoCNGzcyZ84cnJyccHNzY/DgwWzcuJFevXoBMGjQoBRtNm/enIoVKwLQqFEjqlSpQokSJczP88yZMwCcPXuWAQMGcOHCBQoUKMDVq1c5d+4cNWrUyPC1qlWrFj4+PgD4+vpy7tw58/b8/PzM7adXd0hICEeOHGHQoEEAPP744zRo0ACA6OhoNm/ezJ9//mne3u3btzlx4kSqOkaNGpVhnSKSPQrhIhbm6upq/tnZ2Zn4+HgMw8BkMqVaNvk+Nze3VI8/3E7y7eRpkePj49mxYwfz5s1j165dlClThu+//55p06Zlq86s1m1kMM+XYRh06NCBJUuWpHrs/PnzmapPHFO1atU4duwYmzdvZuPGjbz88svmS7mVKlUKLy8vfvjhB3MIz2hfgqSpwpcuXYq/vz9BQUHMmDGDhIQE2rVrl24Nae3rydt6eH9/8PbD8xFk5j0E0LdvX2bOnGn+Q/qxxx5LMSwks3U+fPvu3buPrDut9y8kXUbPZDKxd+/edP9YERHr0phwERuoXbs2MTExbN68GYBdu3Zx5coV81Gp7IqKiqJ48eI89thjxMbG8vHHH1uiXLOM6m7Tpg3r1q0zj0FftGiReb3g4GDWr1+f4uoKv/zyi0Vrk7zpwoULmEwmunTpwsyZMzEMw/wNiaurK99++y3nzp1j1KhRJCYmPnJfateuHZMnT6Zdu3aUKlUKZ2dnVq9enWEIT0/79u1ZuHAhhmEQHR3N0qVLs9XOw6KioqhSpQoAS5cuJSoqKsdtPii9uosXL079+vVZtmwZkPS6HT58GIBixYoRGBjIjBkzzO1cvHiRCxcuWLQ2EUmfQriIDRQqVIhVq1bx+uuv4+3tzdixY1m5ciVubm45ardjx47UqFGD2rVrExISYv7q2lIyqtvb25vXXnuNZs2a0aJFCypUqGBe7x//+AdLly5l5MiRNGzYkDp16jBnzpw0t9GpUyf27dtn0brFfh0+fJhmzZrh7e2Nr68vgwYNwtvb2/x4wYIFWbFiBTExMQwePJiqVatmuC+1b9+ec+fOmcNyUFAQ9+7dy9YfuG+++SYmk4kGDRrQtGlTunTpQs+ePXP8nOfMmUO3bt1o0aIFv/32G5UrV85xmw/KqO4lS5Ywb948fH19+eSTT2jatKl5vWXLlnH8+HEaNGhAgwYN6NGjB9euXUvV/oIFC3TSpogVmIyMvlMWsWN37tzh+PHj1KlThyJFiuR2OSIiko8lfyaFh4cTHh5OtWrV6NatW26XJXZMR8JFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGNGOmOBTnls9Zpd2EbfOs0q5Ibrg6dIJV2nUPnWmVdkVEHJGOhIs4qLCwMPz9/TO17NChQ5k379F/aMyePds8Q6aIJWVlf3UkVapUSTEbaHoy+/qEhYWxYcMGS5QmIlamEC5iZ+Lj43O7hHQphIvY93tUIVwk71AIF7Egk8nEO++8Q9OmTalatSqLFy82P7Zv3z4CAgLw9vamSZMm7Ny5E4Dw8HDc3d2ZNm0agYGBfPDBB7Ru3ZqXXnqJli1bUqlSJd59911WrFhBs2bN8PLyYsWKFeZ2Bw4ciL+/P97e3nTu3DlTITkyMpKgoCC8vb156qmnuHr1qvmx5cuX07RpUxo1aoSPjw/r1q0DYNq0aVy8eJGePXvi4+PDwYMH2bRpEwEBATRq1Ij69euneL4iabl79y59+vShbt26NGzYkODg4FTL3Lhxg/bt2/PWW28B8NNPP9GiRQv8/Pxo2rQp27ZtA6Bfv3588cUXAMydOxcXFxeio6MBCAwMZPv27anaHjp0KKNHj6Zdu3bUrFmT7t27ExsbC8Dt27cZPnw49evXp379+kydOtW8XuvWrXn99dcJCgoiJCSE0NBQgoOD6du3L7Vr16Zt27YcPXqUJ554gpo1a9K3b18SExOB9N9Tj/LGG29Qo0YNWrVqxQ8//GC+//Lly7Rp0wY/Pz/q1avHmDFjMAyDgwcPsmDBApYsWYKPjw/Tpk0jPj6ekJAQ/P39qVevHgMGDODOnTuZ2r6IWJfGhItYmKurK3v27OH48eM0adKEQYMGkZiYSPfu3Vm4cCEhISHs2LGDnj17curUKQCuXbtGjRo1mDRpEgDfffcd58+fJywsjMuXL1O9enVefPFFdu3axS+//ELXrl3p27cvkHR02t3dHYAZM2Ywbdq0Rw4tGTNmDC1btmTy5MmcOXOGhg0b0qFDBwBCQkLo168fJpOJ8PBwmjVrxrlz55g0aRKffvopX3/9NfXr1wcgKiqKHTt24OzszN9//42vry8dOnSgfPnyVnltJe9bv349UVFRHDt2DIC///47xePnz5+nW7dujBs3joEDB3LmzBmmTp3K+vXrKV68OKdOnaJVq1aEh4fTrl07fv75Z/r168emTZvw8/Nj+/bttGjRgqNHj/L444+nWUPyH5CFChWiZcuWrFq1in79+vHWW28RGxvLoUOHuHv3Li1atKBu3br06tXLvN769espWLAgoaGh7N27l8OHD1OxYkU6d+5M//792bZtG25ubvj6+rJx40aCg4PTfU8VLFgw3ddpzZo1fP/99xw8eJDChQunmP68ZMmSrFmzhqJFi5KQkMBTTz3FqlWr6NmzJ6NGjeL27dvMnJk0Pt8wDJYvX07p0qUxDIPRo0czf/58JkywznkBIpJ5CuEiFjZgwAAA6tSpQ4ECBbh8+TJRUVEUKlSIkJAQAFq0aEHZsmU5dOgQ5cuXx9XVlX79+qVop1evXjg5OVGhQgXc3d3p2rUrAH5+fly6dIl79+7h6urKsmXL+Pzzz4mJieHu3bt4eHg8ssYtW7Ywd+5cAKpVq0ZQUJD5sbNnzzJgwAAuXLhAgQIFuHr1KufOnaNGjRqp2rl27RojRozgjz/+MC979OhRhXBJV8OGDTlx4gSjR4+mVatWdOrUyfzYpUuXaNWqFYsXL6Z169ZAUmg/deoULVu2TNFOREQE7du3Z+rUqSQkJHD8+HGmT5/Oxo0bSUhIICAgIN2Q2717dwoXLgxAkyZNOH36NAAbN25kzpw5ODk54ebmxuDBg9m4caM5hA8aNChFm82bN6dixYoANGrUiCpVqlCiRAnz8zxz5gyQtfdUsi1bttCnTx+KFi0KwPDhw/n3v/8NQGJiIq+88go7duzAMAyuXLmCj48PPXv2TNWOYRi8//77rF27lvj4eG7cuJHqtRSR3KHhKCIW5urqav7Z2dmZ+Ph4DMPAZDKlWjb5Pjc3t1SPP9xO8m1nZ2cgaVzqjh07mDdvHj/++COHDx9m1qxZ3Lt3L0f19+3bl1GjRnHkyBEOHjxI0aJF021z1KhRtGrVisOHD3Pw4EFq1qyZ4+2LY6tWrRrHjh2jQ4cO7Ny5k/r16xMVFQVAqVKlqFu3boqhF4Zh0KFDBw4ePGj+FxkZSbVq1ahcuTIuLi4sXboUf39/goKCCAsLY+PGjbRr1y7dGtJ6jyZv6+H34YO3kwNxeu2k125W3lMPPu/0zJo1i2vXrrFnzx4OHTpE//79021v+fLlbN26lW3btnH48GEmTJig96iInVAIF7GB2rVrExMTw+bNmwHYtWsXV65coUGDBjlqNyoqiuLFi/PYY48RGxvLxx9/nKn12rZty6effgokjUnftGlTijarVKkCwNKlS80BCaB48eLcuHEjxbJeXl6YTCa2bdvGb7/9lqPnI47vwoULmEwmunTpwsyZMzEMg4iICCAp1H777becO3eOUaNGkZiYSHBwMOvXr09xBZFffvnF/HO7du2YPHky7dq1o1SpUjg7O7N69eoMQ3h62rdvz8KFCzEMg+joaJYuXZqtdh6W0XsqPUFBQXz11VdER0eTkJBAaGhoivY8PDxwdXXlzz//ZOXKlebH0nqPli5dmmLFinHr1q0U7YhI7tJwFHEo9no970KFCrFq1SrGjBlDdHQ0rq6urFy5Ejc3N/76669st9uxY0eWLl1K7dq1qVixIs2aNeOnn3565Hpz5sxh8ODBrFy5kpo1a6YIGnPmzKFbt254enoSEBBA5cqVzY+NGTOGYcOGUaRIEUJDQ5kxYwajR49mxowZ1K1bl6ZNm2b7uYjt5Ob1vA8fPszEiRMxDIPExEQGDRqEt7c3YWFhABQsWJAVK1YwcuRIBg8eTGhoKEuXLmXkyJHcvXuX2NhYfH19WbZsGZAUnBcsWGDeh4OCgli0aFG2/sB98803ef75583r9urVK80hHlmV0XsqPZ07d2b37t00bNgQT09PWrVqxYULF4Ck92GvXr3w8fHB09Mzxfu3W7dufP755/j4+NC9e3deeOEFvvvuO+rWrYunpyeBgYFERkbm+DmJSM6ZjIy+8xKxY3fu3OH48ePUqVOHIkWK5HY5IiKSjyV/JoWHhxMeHk61atVSnFAr8jANRxERERERsTGFcBERERERG1MIFxERERGxMYVwEREREREbUwgXEREREbExhXARERERERvTdcLFoazefMEq7XZvW9Eq7VpTWFgYEyZMYN++fY9cdujQofj7+/Pcc89luNzs2bPp378/ZcuWzXI9d+7coXnz5mzbto1ixYpRpUoVfvjhB+rXr59q2ZEjRzJkyBACAwMtUretTJo0iXr16tGnT59srd+zZ0/Gjx9Ps2bNLFxZSnFrB1il3YJPLMv2ulnZXx1JRu+DB2X29QkLCyM2Npbg4OBs1TN+/HiaNGlC3759mTJlCrdv32bmzNTXlf/+++/Zvn077777rkXqtpV9+/bx/vvvm68zn1Vr1qzhhx9+yPTEaCIZ0ZFwETuTPNW1PZo9ezZXrlzJ1rrz5s2jW7duFCtW7JHL/ve//31kALcUS77e06ZNy3YAB3jttdd47bXXLFaPWIc9v0fDwsLYsGFDttaNjIzkxx9/zNQ+3KVLl0cGcEuy1Gvu7++f7QAO8OSTT7J3715Onz5tkXokf1MIF7Egk8nEO++8Q9OmTalatSqLFy82P7Zv3z4CAgLw9vamSZMm7Ny5E0iaNt7d3Z1p06YRGBjIBx98QOvWrXnppZdo2bIllSpV4t1332XFihU0a9YMLy8vVqxYYW534MCB+Pv74+3tTefOnTMVkiMjIwkKCsLb25unnnqKq1evmh9bvnw5TZs2pVGjRvj4+LBu3TogKWBevHiRnj174uPjw8GDB9m0aRMBAQE0atSI+vXrp3i+D/vkk08YMCDlEdhVq1bRrFkzqlatyr///W/z/a1bt+aHH35IUWu9evXo3LkznTt3Zt68+zOjHjt2jHbt2lGzZk26d+9ObGwsAHFxcUycOJEmTZrg4+ND3759uX79OpB0BH3MmDF06NCBhg0bpqp1ypQp9O/fnyeffJK6devStm1b/v77bwASEhKYMGEC9evXp379+jz//PPmbQ4dOtRc25o1a/D29sbHx4f69evz3XffAXD58mV69+5NkyZN8Pb2ZtKkSebt+vr6cvnyZU6ePJnRry9Pu3v3Ln369KFu3bo0bNgwzSO2N27coH379rz11lsA/PTTT7Ro0QI/Pz+aNm3Ktm3bAOjXrx9ffPEFAHPnzsXFxYXo6GgAAgMD2b59e6q2hw4dyujRo9PcZ27fvs3w4cPNv9upU6ea12vdujWvv/46QUFBhISEEBoaSnBwMH379qV27dq0bduWo0eP8sQTT1CzZk369u1LYmIikP576lHeeOMNatSoQatWrczvB0jah9q0aYOfnx/16tVjzJgxGIbBwYMHWbBgAUuWLMHHx4dp06YRHx9PSEgI/v7+1KtXjwEDBnDnzp00t/fpp5/Ss2dPTCaT+b6LFy+m+T4IDQ1NMZvo66+/To0aNWjatCkvvfQS/v7+5sfi4+MZPXo0DRs2pF69eimOiqf3uw0LC8PHx4cxY8YQEBDAN998k6LW5H5z0qRJ+Pn5UaNGjRSv6/r16/H19cXb25tWrVpx7Ngxc7vJtf31118EBwfToEEDvL29GTZsmHn9mTNn0qRJE3x9fenUqRMRERHmx3r37s2nn376qF+fyCMphItYmKurK3v27GHdunWMGTOG+Ph4YmNj6d69O1OmTOHQoUPMmjWLnj17mgPDtWvXqFGjBtu3b2fcuHEAnD9/nrCwMPbs2cOkSZM4cuQIu3btYuXKlYwfP968vdmzZ7Nv3z4OHTpEixYtmDZt2iNrHDNmDC1btuTQoUO8//77bN261fxYSEgI//vf/zhw4ADffvstI0eOJC4ujkmTJlGhQgW+/vprDh48iI+PD76+vuzYsYMDBw6wbds2pk6dyqVLl1JtLyIigps3b1K9evUU91+/fp1du3bxyy+/8O6776Y5nfaYMWNo06YNR48eZf78+eYP6WQHDx5kzZo1HD9+nD///JNVq1YB8O6771K0aFF++eUXDh48SL169Zg8ebJ5vR07dvD1119z9OjRNF+jPXv28Nlnn3Hs2DHKli1r/vr5k08+4ddff+XXX3/l4MGDnD59mjlz5qRa/4033mDBggUcPHiQQ4cO0apVKwCGDBnCc889xy+//ML+/fv55ZdfUgSMZs2asWnTpjRrcgTr168nKiqKY8eO8dtvv6X4gxKS9vu2bdsyZMgQ3nzzTc6cOcPUqVNZt24dv/76K8uWLaNfv37ExcXRrl07fv75ZwA2bdqEn58f27dv5/bt2xw9epTHH388zRrS22feeustYmNjOXToEHv27OHbb79l5cqVKdZbv369+fezd+9eZs6cyYkTJyhSpAj9+/dn+fLlHDt2jGPHjrFx40Yg/fdURtasWcP333/PwYMH2bx5M3/88Yf5sZIlS7JmzRp+/fVXDh06xJkzZ1i1ahU+Pj6MGjWKwYMHc/DgQSZNmoSzszPLly9n3759HDlyhOLFizN//vw0txkWFpZqKFR674OHa/3hhx/47bff2L17d6qjxEePHmX48OH89ttvPP/887z++usAGf5uAQ4dOkTv3r3ZvXs3vXr1SrXda9eu4efnx6+//sq8efPMfeeVK1cYOHAgn332GYcOHeLpp5+md+/eqdZfunQpVapU4fDhwxw6dIj33nsPSPqj6Y8//mD37t3s37+ffv36pRjy5ujvUbEdjQkXsbDko7116tShQIECXL58maioKAoVKkRISAgALVq0oGzZshw6dIjy5cvj6upKv379UrTTq1cvnJycqFChAu7u7nTt2hUAPz8/Ll26xL1793B1dWXZsmV8/vnnxMTEcPfuXTw8PB5Z45YtW5g7dy4A1apVIygoyPzY2bNnGTBgABcuXKBAgQJcvXqVc+fOUaNGjVTtXLt2jREjRvDHH3+Ylz169Cjly5dPsdyFCxdS3ffga1WmTBmqVavG2bNn8fT0TLfWypUrp6gVoHv37hQuXBiAJk2amAPAt99+y82bN/n6668BiI2NTfFHQO/evSlatGi6r1HHjh157LHHAAgICODw4cMAbNy4kREjRuDi4gLAP//5TxYsWMBLL72UYv2goCDGjh1Lz549CQ4OxsfHh+joaDZv3syff/5pXu727ducOHHCfNvDw4MLF6xzboM9aNiwISdOnGD06NG0atWKTp06mR+7dOkSrVq1YvHixbRu3RpICu2nTp2iZcuWKdqJiIigffv2TJ06lYSEBI4fP8706dPZuHEjCQkJBAQEULBgwTRrSG+f2bhxI3PmzMHJyQk3NzcGDx7Mxo0bzQFw0KBBKdps3rw5FSsmnS/SqFEjqlSpQokSJczP88yZM0DW3lPJtmzZQp8+fcz76PDhw83fFiUmJvLKK6+wY8cODMPgypUr+Pj4pDgyncwwDN5//33Wrl1LfHw8N27cSPVaJrtw4UKq/iO998HDtfbu3Rs3Nzcg6Q/N5G8xAGrVqmU++hwQEGAeY57R7xagZs2atGjRIt3XyM3NjaeeesrcbvLvcc+ePfj4+NCgQQMgqZ/517/+leoAweOPP87777/Piy++SKtWrcz987fffsu+ffvw8/MDkr79cnZ2Nq/n6O9RsR2FcBELc3V1Nf/s7OxMfHw8hmGk+Io3WfJ9bm5uqR5/uJ3k28kfBvHx8ezYsYN58+axa9cuypQpw/fff5+pI+EZ6du3LzNnzjSH/scee4x79+6lueyoUaN48sknWbVqFSaTCV9f3zSXLVKkCHfv3k11f1qvVVrSeu0e1YZhGMyfP5+2bdumuV5GAfxR7T5cT1r1zZo1i6NHj7JlyxaGDBnCgAEDePbZZzGZTOzduzfdgHjv3j1Kly6dYW15WbVq1Th27BibN29m48aNvPzyyxw8eBCAUqVK4eXlxQ8//GAO4YZh0KFDB5YsWZJmey4uLixduhR/f3+CgoKYMWMGCQkJtGvXLt0asvu7fXifSe89+nC7WXlPJTMMI93HZs2axbVr19izZw+urq6MHz8+3faWL1/O1q1bzSdEz507N9W3ScnSep9m5j2aXv/2qDYy+t2eP38+y+/RhISEDOt5+L6AgAAOHjzIxo0bWbVqFW+88QYHDhzAMAzeeOMNhg8fnuZ27927Z/4jTiQnNBxFxAZq165NTEwMmzdvBmDXrl1cuXLFfKQmu6KioihevDiPPfYYsbGxmT5jv23btuYxjeHh4Sm+Wo2KiqJKlSpA0te1UVFR5seKFy/OjRs3Uizr5eWFyWRi27Zt/Pbbb2lur1atWvz555+PDB5pad26NaGhoUDSEbLk1/BRunTpwqxZs8zjX+/cuZPu0JOsaN++PaGhocTGxhIfH8+iRYvSDHwnTpygXr16PPfcczz77LP873//o1ixYgQGBjJjxgzzchcvXkxxVO348eNpjlN3FBcuXMBkMtGlSxdmzpyJYRjmI5+urq58++23nDt3jlGjRpGYmEhwcDDr16/nyJEj5jZ++eUX88/t2rVj8uTJtGvXjlKlSuHs7Mzq1aszDOHpad++PQsXLsQwDKKjo1m6dGm22nlYRu+p9AQFBfHVV18RHR1NQkKC+T2Q3J6Hhweurq78+eefKYbMpPUeLV26NMWKFePWrVsp2nmYt7d3im9lMqtNmzasXLmSO3fukJiYyOeff56p9R71u82u5HB9/PhxAFasWEHFihVTHeU/e/YsRYsWpXfv3nzwwQf88ccf3L59my5dujB//nzz+Pe4uDgOHDhgXs/R36NiOzoSLg7FXi8lWKhQIVatWsWYMWOIjo7G1dWVlStX4ubmxl9//ZXtdjt27MjSpUupXbs2FStWpFmzZvz000+PXG/OnDkMHjyYlStXUrNmzRRBY86cOXTr1g1PT08CAgKoXLmy+bExY8YwbNgwihQpQmhoKDNmzGD06NHMmDGDunXr0rRp0zS35+rqSrt27di0aRNPPPFElp5jcq1ffvklNWvWpHnz5uav/DMyceJEpk6dStOmTc1HwF555RXq1auXpe0/7Omnn+b06dP4+voCSX8kjBkzJtVyr776Kn/88QeFChWiSJEifPTRRwAsW7aM8ePHm/8AK1q0KAsWLKBixYpER0dz9OjRdI/eW0pOLiWYU4cPH2bixIkYhkFiYiKDBg3C29ubsLCwpNoKFmTFihWMHDmSwYMHExoaytKlSxk5ciR3794lNjYWX19f8xUu2rdvz4IFC8z7cFBQEIsWLcrWH7hvvvkmzz//vHndXr16pTnEI6syek+lp3PnzuzevZuGDRvi6elJq1atzH+sjRkzhl69euHj44Onp2eK92+3bt34/PPP8fHxoXv37rzwwgt899131K1bF09PTwIDA9M89wKSLpG5ZMmSdI8Ap6dLly7s2rWLhg0bUqFCBR5//PFM/aHxj3/8I8PfbXaVKVOGzz//nAEDBpCQkEDJkiX56quvUi0XFhbGrFmzzEfR3333XUqUKMGgQYO4du0arVu3xmQyER8fz4gRI2jUqBGQNIymR48eOapRBMBkZPSdl4gdu3PnDsePH6dOnToUKVIkt8uRR9izZw9vvfVWiqs8ZMbdu3cpWLAgBQoU4NKlSzRu3JhNmzZRq1YtK1WaexYsWEBkZGSK8bQitpKYmEjjxo357rvvzGPdM+vWrVsUK1aMxMRERo4cSYUKFVJc8chRXL16laCgIPbu3UuhQoVSPJb8mRQeHk54eDjVqlWjW7duuVSp5AU6Ei4iNtG0aVO6du1q/rDOrJMnTzJ48GAMwyAuLo7Jkyc7ZAAHcHJyYuLEibldhuRTTk5OfPzxx4SHh2c5hA8ePJjw8HDu3r2Lr68vL7/8spWqzF2nT5/mo48+ShXARbJDR8Ilz9KRcBERsRc6Ei5ZpRMzRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMZ0dRRxKN++tYGVYctZvvkz+rcdQq/W/TNcfuP+n5j/3fsE+3fi6Seew8kp7b9Lu74ZbI1yRXKFyWRiykRvXnuxfobLhS47zTPj9vDPITWY+05jnJzSnxVx3YZInuofZuFKRUQcl46Ei0OxRgDPqbCwMPz9/a3Stj2rUqVKipnw0pPZ1ycsLIwNGzZYorR8zxoBvPew7ZYu0yay8v4cOnQo8+bNe+Rys2fP5sqVKzktTUQcnEK4OBRrBPB9v++xdJl2Kz4+PrdLSJdCuOVYI4B3aFfB0mVmyJ73VYVwEckMhXBxKNYI4O+syNzshXfv3qVPnz7UrVuXhg0bEhycegjLjRs3aN++vXlGxJ9++okWLVrg5+dH06ZN2bZtGwD9+vXjiy++AGDu3Lm4uLgQHR0NQGBgINu3pz7qOHToUEaPHk27du2oWbMm3bt3JzY2FoDbt28zfPhw6tevT/369Zk6dap5vdatW/P6668TFBRESEgIoaGhBAcH07dvX2rXrk3btm05evQoTzzxBDVr1qRv374kJiYCsHz5cpo2bUqjRo3w8fFh3bp1mXqt3njjDWrUqEGrVq1SzKB5+fJl2rRpg5+fH/Xq1WPMmDEYhsHBgwdZsGABS5YswcfHh2nTphEfH09ISAj+/v7Uq1ePAQMGcOfOnUxtX9KX3QC+fGHzTG/DZDLxzjvv0LRpU6pWrcrixYvNj+3bt4+AgAC8vb1p0qQJO3fuBCA8PBx3d3emTZtGYGAgH3zwAa1bt+all16iZcuWVKpUiXfffZcVK1bQrFkzvLy8WLFihbndgQMH4u/vj7e3N507d85USI6MjCQoKAhvb2+eeuoprl69an4svX1/2rRpXLx4kZ49e+Lj48PBgwfZtGkTAQEBNGrUiPr166d4viKSf2lMuDgUawRwv5qNM7Xt9evXExUVxbFjxwD4+++/Uzx+/vx5unXrxrhx4xg4cCBnzpxh6tSprF+/nuLFi3Pq1ClatWpFeHg47dq14+eff6Zfv35s2rQJPz8/tm/fTosWLTh69CiPP/54mjUkf+AXKlSIli1bsmrVKvr168dbb71FbGwshw4d4u7du7Ro0YK6devSq1cv83rr16+nYMGChIaGsnfvXg4fPkzFihXp3Lkz/fv3Z9u2bbi5ueHr68vGjRsJDg4mJCSEfv36YTKZCA8Pp1mzZpw7d46CBQum+zqtWbOG77//noMHD1K4cOEUk1mULFmSNWvWULRoURISEnjqqadYtWoVPXv2ZNSoUdy+fZuZM2cCYBgGy5cvp3Tp0hiGwejRo5k/fz4TJkzI1O9LUstJAC9UyDlL23J1dWXPnj0cP36cJk2aMGjQIBITE+nevTsLFy4kJCSEHTt20LNnT06dOgXAtWvXqFGjBpMmTQLgu+++4/z584SFhXH58mWqV6/Oiy++yK5du/jll1/o2rUrffv2BZKOTru7uwMwY8YMpk2b9sihJWPGjKFly5ZMnjyZM2fO0LBhQzp06ACQ7r4/adIkPv30U77++mvq10/6xiEqKoodO3bg7OzM33//ja+vLx06dKB8+fJZes1ExLEohEu+kd0A/mKv1zLVfsOGDTlx4gSjR4+mVatWdOrUyfzYpUuXaNWqFYsXL6Z169ZAUmg/deoULVu2TNFOREQE7du3Z+rUqSQkJHD8+HGmT5/Oxo0bSUhIICAgIN2Q2717dwoXLgxAkyZNOH36dNJz37iROXPm4OTkhJubG4MHD2bjxo3mED5o0KAUbTZv3tw8bXWjRo2oUqUKJUqUMD/PM2fOAHD27FkGDBjAhQsXKFCgAFevXuXcuXPUqFEj3ddpy5Yt9OnTh6JFiwIwfPhw/v3vfwOQmJjIK6+8wo4dOzAMgytXruDj40PPnj1TtWMYBu+//z5r164lPj6eGzdupHotJfNsGcABBgwYAECdOnUoUKAAly9fJioqikKFChESEgJAixYtKFu2LIcOHaJ8+fK4urrSr1+/FO306tULJycnKlSogLu7O127dgXAz8+PS5cuce/ePVxdXVm2bBmff/45MTEx3L17Fw8Pj0fWuGXLFubOnQtAtWrVCAoKMj+WlX3/2rVrjBgxgj/++MO87NGjRxXCRfI5DUeRfCEnAbxggfSP6j6oWrVqHDt2jA4dOrBz507q169PVFQUAKVKlaJu3bophl4YhkGHDh04ePCg+V9kZCTVqlWjcuXKuLi4sHTpUvz9/QkKCiIsLIyNGzfSrl27dGtwdXU1/+zs7GweN2sYBiZTylD14O3kQJxeO+m127dvX0aNGsWRI0c4ePAgRYsW5d69exm+ToZhpPvYrFmzuHbtGnv27OHQoUP0798/3faWL1/O1q1b2bZtG4cPH2bChAmP3LakzdYBHNLeV9PaT+H+vurm5pbq8fT2VWfnpLri4+PZsWMH8+bN48cff+Tw4cPMmjUrx/tKVvb9UaNG0apVKw4fPszBgwepWbOm9lURUQgXx2eLAA5w4cIFTCYTXbp0YebMmRiGQUREBJAUFL799lvOnTvHqFGjSExMJDg4mPXr16e4gsgvv/xi/rldu3ZMnjyZdu3aUapUKZydnVm9enWGITw97du3Z+HChRiGQXR0NEuXLs1WOw+LioqiSpUqACxdutT8R0dGgoKC+Oqrr4iOjiYhIYHQ0NAU7Xl4eODq6sqff/7JypUrzY8VL16cGzdupFi2dOnSFCtWjFu3bqVoRzIvNwJ4emrXrk1MTAybN28GYNeuXVy5coUGDRrkqN2oqCiKFy/OY489RmxsLB9//HGm1mvbti2ffvopkDQmfdOmTSnaTG/fT2tf9fLywmQysW3bNn777bccPR8RcQwajiIOJa3reXclmHm8l7n1CebfS9/M1rYPHz7MxIkTMQyDxMREBg0ahLe3N2FhYQAULFiQFStWMHLkSAYPHkxoaChLly5l5MiR3L17l9jYWHx9fVm2bBmQFJwXLFhgDstBQUEsWrQoW4HkzTff5Pnnnzev26tXrzSHeGTVnDlz6NatG56engQEBFC5cuVHrtO5c2d2795Nw4YN8fT0pFWrVly4cAFIGoPbq1cvfHx88PT0TPGHQrdu3fj888/x8fGhe/fuvPDCC3z33XfUrVsXT09PAgMDiYyMzPFzyg8Kui8z//zPF5L+ZcZT/SEm49MucqRQoUKsWrWKMWPGEB0djaurKytXrsTNzY2//vor2+127NiRpUuXUrt2bSpWrEizZs346aefHrnenDlzGDx4MCtXrqRmzZop9seM9v0xY8YwbNgwihQpQmhoKDNmzGD06NHMmDGDunXr0rRp02w/FxFxHCYjo++GRezYnTt3OH78OHXq1KFIkSK5XY6IiORjyZ9J4eHhhIeHU61atRQnnos8TMNRRERERERsTCFcRERERMTGFMIlz0ueOEZERCS36LNIskonZkqe5eLigslk4tKlS5QvXz7Dq56IiIhYS2JiIhcvXsQwDOLi4nK7HMkjFMIlz3J2dqZGjRqcOnWKmzdv5nY5IiKSjxmGwYULF0hMTMxwPgSRZArhkqcVL16cypUr88033+Di4oKbm1tulyQiIvlQXFycOYDHxMTo80geSSFc8rxSpUrh5eXFr7/+irOzs4aliIhIrjAMg9jYWEqVKkXt2rVzuxyxc7pOuDiEuLg4Dhw4wNWrV81TqouIiNiam5sb//jHPzI1eZnkbwrhIiIiIiI2pu/tRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGFMJFRERERGxMIVxERERExMYUwkVEREREbEwhXERERETExhTCRURERERsTCFcRERERMTGFMJFRERERGzM5iF81qxZmEwmTCYT8+bNs/XmRUTEhu7evYuHhwcmk4nKlStz9+7d3C5JRMQuZDuEh4WFmcP0g/8KFChAiRIlaNiwIc899xx//PFHivVGjRpFmTJlAJgyZQrXr1/P0RMQERHbunTpElOnTqVly5aUK1eOQoUKUa5cOfz8/HjhhRfYvXu3ednChQszYcIEACIiIpg5c2ZulS0iYldMhmEY2VkxLCyMNm3aPHK5IkWKEBYWRuPGjc33TZo0ibfeeguAN954w/yziIjYtw8++ICXXnqJmJiYDJeLioqiZMmSANy+fZty5cpx584dihYtyoULFyhRooQNqhURsV8WG47Sp08f3n33Xd5++226detmvv/OnTv83//9X4plBwwYYP554cKFxMbGWqoMERGxkrfffpsxY8aYA3iBAgXo2rUrU6dOZerUqQwcOBB3d/dU6xUtWpQuXboASYH8s88+s2ndIiL2qIClGurQoQNDhw413/b29ubw4cMAnDhxIsWytWrVon79+hw5coQ///yTtWvXpgjuIiJiX44ePcqbb75pvl22bFl++uknfHx8UiwXExPDJ598QsGCBVPc37NnT1asWAHAokWLGDNmjNVrFhGxZxY/MTMhIYH//e9/nDt3znyfh4dHquUCAgLMP2/atMnSZYiIiAV98MEHJCQkmG9/9NFHqQI4gIuLC88//zxubm4p7n+wzz98+DB//fWX1WoVEckLLBbChw0bZj4xMyAggJs3byZtwMmJl156KdXyjRo1Mv+8c+dOS5UhIiJWsHnzZvPPpUqVomvXrllav0KFCpQtWxYAwzDYtWuXJcsTEclzLDYcJT3//ve/eeKJJ1Ldn9wZA0RGRlq7jExJTEzk4sWLFCtWDJPJlNvliEguMwyDW7duUaFCBZyc8ve0Cg/20zVr1szW61G2bFmuXLmSqj1bUP8uIg+yh/7dYiG8T58++Pv7k5CQwJEjR1ixYgXx8fG89tprxMXFMWnSpBTLFy9e3PyzvVym8OLFi1SqVCm3yxAROxMREUHFihVzu4w8Lzf7ffXvIpKW3OzfrXZiZvXq1Zk6dSoAb731FiNGjMDT09P8ePJwFcB8GavcVqxYMSDpF/Lgh4WI5E83b96kUqVK5r4hP/P09OTkyZMA/PHHHxiGkeUjyrnZ76t/F5EH2UP/brXhKE2aNDH/HB8fz969e1OE8OSvJIEU9+em5A+U4sWLq5MWETMNX4C2bduaQ3hUVBTfffddlseF52a/r/5dRNKSm/271QbB7N27N8XtB8+qB9i/f7/55+bNm1urDBERsYDnnnsuxbjJZ599lkOHDqVaLjY2lnnz5hEdHZ3i/osXL5pDuMlkolmzZtYtWETEzlnsSPj69eu5evUqCQkJHDt2jOXLl5sfc3Z2pmnTpimWf3Ba46CgIEuVISIiVlC/fn2mTp1qvlb45cuX8fPzo0uXLuZLFf7xxx/mz4KBAwemWP/Bq6E0aNCAMmXK2Kx2ERF7ZLEQ/uWXX/Lll1+m+dikSZNSDHo/ceIER48eBaBcuXJpXj1FRETsyxtvvIGbmxsTJ04kNjaW+Ph4Vq9ezerVqx+57tdff23+ecSIEdYsUyRN0dHRFC1aFEiaufXha9mL2JpVxoS7uLjg4eFB48aNGTlyJCEhISkeX7Zsmfnnf/7znxQqVMgaZYiIiIWNGzeO3r17s3DhQjZu3Mgff/zB9evXKVmyJJUqVaJ58+b07ds3xYmXt27d4vvvvweSprAfMmRILlUvImI/sh3CW7dujWEYWV7vzp07fPzxxwCULl2aF198MbsliIhILvD09GTKlClMmTIlU8svWLCAu3fvAvDKK69QokQJK1Zn33Q0VkSS2fzq5AsWLDBPVzxlyhS7uTyhiIhY3t27d3nvvfcAqFSpkg68iIj8f1afMfNh48ePZ/z48bberIiI5ILChQtz+fLl3C5DRMTu5O95mEVEREREcoFCuIiIiIiIjSmEi4iIiIjYmEK4iIiIiFhNdHQ0JpMJk8mUajbd/EwhXERERETExhTCRURERERsTCFcRERERMTGbH6dcBEREZGMOLd8zuJtGglx5p+LBY/H5FzQ4tsASNg2zyrtiuPRkXARERERERuz2yPh27Zt49133+XXX3/l0qVLfPPNN3Tt2jXDdbZu3cr48eM5evQoFSpU4OWXX2bUqFG2KVhERB4pL/XtOhorItZkt0fCo6OjadiwIfPmZa4jOXv2LJ06dSIwMJADBw7w2muvMWbMGFatWmXlSkVEJLPUt4uIJLHbI+EdO3akY8eOmV5+wYIFVK5cmdmzZwNQp04d9u3bx8yZM+nRo4eVqhQRkaxQ3y4iksRuQ3hW7d69m+Dg4BT3hYSEsGjRIuLi4ihYMPVXfjExMcTExJhv37x50+p12qPo6GiKFi0KwO3bt3Fzc8vlikRS0j6af2Wnbwf17yLZYY0hWKBhWOlxmBB++fJlypUrl+K+cuXKER8fz9WrVylfvnyqdd5++22mTp2ao+1aa4fNqry244lt2ct+unLKxGytd+/uHfPP34VF4lq4SLZr6N62YrbXFdvLTt8OlunfRUSsyWFCOIDJZEpx2zCMNO9P9uqrrzJ+/Hjz7Zs3b1KpUiXrFWhFqzdfyPa6Cjhi71wLF2HVpojcLkNySVb7dnCs/l1EHJPDhHAPDw8uX76c4r4rV65QoEABSpcuneY6Li4uuLi42KI8ERHJhuz07aD+XUTsn8OE8ICAANasWZPivg0bNuDv75/umEFJoqOMIpKWM2fO8N1337F3717+/vtvSpUqha+vL0899RQ1a9a0SQ3q20XEUdltCL99+zanTp0y3z579iwHDx7kscceo3Llyrz66qtERkayZMkSAEaNGsW8efMYP348//znP9m9ezeLFi3iiy++yK2nICKSJx0/fpxXXnmFdevWYRiGefgHwFdffcXEiRPp0KED7777LnXr1s1S2+rbRUSS2G0I37dvH23atDHfTh7bN2TIEEJDQ7l06RLnz583P161alXWrVvHuHHj+PDDD6lQoQJz587VJaxERLLI29ubhIQEqlatSps2bahbty7Fixfn5s2bHDt2jC1btvDjjz+ycePGFFcgyQz17SIiSew2hLdu3TrF0ZeHhYaGprqvVatW7N+/34pViYg4vo4dO/LKK6/QvHnzdJfZuXMn//nPf7Lctvp2EZEkdhvCRUQkd3z//fePXKZ58+Z89913NqjGsZicC2IK/Fdul5Ev6bUXe2OxEG4PJ/CIiIh1/PTTT+zZs4dJkybldikiIg4hxyHcmifwiIiIfVi/fj1z585VCBeRLNO3EGnLcQi35gk8IiIiIiKOKMch3Jon8IiIiIiIOKIch3CdwCMi4vj69++Pj49PbpchIuIwdHUUERF5pMaNG9O4cePcLkNExGE4WaPRn376iWnTplmjaRERsbJx48YRGRmZ4TKRkZGMGzfORhWJiDgeqxwJ11n0IiJ515w5c5g3bx6tWrWibdu2qU6437RpE9u2bSMxMZH3338/t8sVEcmTNBxFRERSWLt2LS+99BKbN29my5YtqR43DIM6deowc+bMXKhORMQxKISLiEgKHTt2pGPHjvz8889888037N+/n6ioKEqWLImvry9du3YlJCQkt8sUEcnTrBLCdRa9iEje1759e9q3b5/bZYiIOKQcn5iZ1gk8jRs3ZsiQIebbOoFHREREROS+HB8J1wk8IiIiIiJZk+MQrhN4RERERESyxiLT1usEHhERERGRzLPYiZk6gUdEREREJHOsMmOmiIg4liNHjhAQEECxYsXw8fFh27ZtuV2SiEiephAuIiKPNHz4cI4cOUKDBg24ePEiffr0ye2SRETyNE3WIyIij/Tbb7+xc+dO/P39iYqKoly5ckRGRuLp6ZnbpYmI5EkK4SIi8khxcXH8+uuvnDlzBoACBQrwzTffULZsWXr37p3L1YmI5D0K4SIikinPPvtsittjxozBZDIphIuIZINVQ/jp06dZvHgxUVFRGIaByWTiww8/tOYmRUTECtKaB0JERLLPqiG8Y8eOnDp1ynxbIVxEJG967LHHePrppzly5AjVq1dn7ty5tGzZMrfLEhHJs6x6dZTbt2+zdetWEhMTSUxMJCEhwZqbExERKxkxYoSujiIiYkFWPRLevHlz5s2bR2RkJE5OTphMJnr16mXNTYqIiBXo6igiIpZl1RC+atUqAFauXAmgEC4ikkfp6igiIpZl1RCuE3lERBzH6NGjMQzDfFtXRxERyT6rhvCzZ8/i7e3NH3/8wUcffcSYMWOsuTkREbESHVQREbEsq4bw119/nW+++YbXXnuNMmXKMGHCBHr06GHNTYqIiBX4+PhQokSJ3C5DRMRhWPXqKNeuXQMgKiqK5cuXc+nSJWtuTkRErKRhw4YEBQUxe/Zs87hwERHJPqseCa9Zsybt27cnICCAP//8k3LlyllzcyIiYiXh4eHs37+fNWvW0Lt3b+7du0fnzp3p0qULAQEBmEym3C5RRCRPseqR8K+//ppJkyaxYMECihYtyrx587K0/vz586latSqurq74+fmxffv2dJcNCwvDZDKl+nfixImcPg0REQF8fX2ZPHky+/btY8OGDVSrVo0ZM2bw/vvvZ7kt9e8ikt9Z5Uj4jRs3KFGiBDVr1uTFF1803+/t7Z3pNr788kvGjh3L/Pnzad68OR9//DEdO3bk2LFjVK5cOd31fv/9d4oXL26+XaZMmew9CRERMUtMTOT999/nv//9L+fPn8fT05Onn36a7777LstDDdW/i4hY6Ui4JcYOzpo1ixEjRjBy5Ejq1KnD7NmzqVSpEh999FGG65UtWxYPDw/zP2dn53SXjYmJ4ebNmyn+iYhIau+88w4vvfQS0dHR+Pj4EBMTwyuvvMJzzz1Hu3btstSW+ncRESuF8PDwcN59911u3LhB7969qV+/PhMnTmTXrl0prjGbntjYWH799VeCg4NT3B8cHMyuXbsyXLdRo0aUL1+eoKCgR15S6+2336ZEiRLmf5UqVXr0kxMRyYc+/fRTpk+fzvnz59m5cyfnzp1j+vTpLFiwgObNm2e6HfXvIiJJrDYmPK2xg++8806mxg5evXqVhISEVCdylitXjsuXL6e5Tvny5fnkk09YtWoVq1evplatWgQFBbFt27Z0t/Pqq69y48YN87+IiIisPUkRkXzizz//5Mknn0xxX+fOnXFzc2PhwoWZbkf9u4hIEqteHeX06dMsXryYqKgoDMPA09OT8ePHZ3r9h8+2Nwwj3TPwa9WqRa1atcy3AwICiIiIYObMmbRs2TLNdVxcXHBxccl0PSIi+VWDBg3417/+xZQpU6hYsSIXL17krbfeytK5Pg9S/y4i+Z1VQ3jHjh05deqU+bbJZGL+/PmPXM/d3R1nZ+dUR0WuXLmSpcscPv744yxdujTzBYuISJr+85//EBwcTFBQkPk+FxcXNmzYkKV21L+LiCSx6iUKb9++zdatW0lMTCQxMZGEhIRMrVeoUCH8/Pz4+eefU9z/888/06xZs0xv/8CBA5QvXz5LNYuISGrNmzfn2LFj/Pvf/2bUqFH83//9H8eOHaNFixZZakf9u4hIEqseCW/evDnz5s0jMjISJycnTCYTvXr1ytS648ePZ9CgQfj7+xMQEMAnn3zC+fPnGTVqFJA03i8yMpIlS5YAMHv2bKpUqUK9evWIjY1l6dKlrFq1ilWrVlnt+YmI5CdeXl68+uqrOW5H/buIiJVDeHIHuXLlSoAshfA+ffpw7do1pk2bxqVLl6hfvz7r1q3Dy8sLgEuXLnH+/Hnz8rGxsUyYMIHIyEgKFy5MvXr1WLt2LZ06dbLwsxIRyX+qVauW5v0mk4nTp09nqS317yIiVg7hj7qE1KOMHj2a0aNHp/lYaGhoitsvv/wyL7/8co62JyIiaQsPD8dkMtGkSRNKliyZ4/bUv4tIfmeRED537lwaNWpEYGAgZ86c4cKFC7Rs2ZJWrVpZonkREcllM2fOJDQ0lIMHD9K9e3eGDx+e4iRNERHJGoucmDl27Fi++eYbAD744APatGljiWZFRMROjB8/nkOHDrFt2zbOnDlDcHAwI0eOzO2yRETyLKsORxEREcdw9epVPv/8cxYvXsypU6fo27cvTz/9dG6XJSKSZymEi4jII1WsWJG4uDhatGjB7NmzKVmyJOHh4YSHh9O7d+/cLk9EJM+xWAg/fPgwn3zyCYcPHwZg4cKFGIZhflxHTERE8q7Y2FgAtm/fzo4dO4D7s1wqhIuIZJ3FQvjmzZvZvHmz+Xby9V6TKYSLiORdOb3alYiIpGSREF65cmVMJpMlmhIRETvw8FWvTCYTLVu2zO2yREQchkVCeHh4uCWaEREROzF27FjGjh1LYGAgH3zwAXPnziUhISG3yxIRcRgWuUShiIiIiIhknkK4iIiIiIiN6RKFIiKSJl31SkTEehTCRUQkTbrqlYiI9SiEi4hIKrrqlYiIdSmEi4hIKrrqlYiIdenETBERERERG1MIFxERERGxMYVwEREREREbUwgXEREREbExhXARERERERtTCBcRERERsTGFcBERERERG1MIFxERERGxMYVwEREREREbUwgXEREREbExhXARERERERtTCBcRERERsTGFcBERERERG1MIFxERERGxMYVwEREREREbUwgXEREREbExhXARERERERuz6xA+f/58qlatiqurK35+fmzfvj3D5bdu3Yqfnx+urq5Uq1aNBQsW2KhSERHJCvXvIpLf2W0I//LLLxk7diyvv/46Bw4cIDAwkI4dO3L+/Pk0lz979iydOnUiMDCQAwcO8NprrzFmzBhWrVpl48pFRCQj6t9FROw4hM+aNYsRI0YwcuRI6tSpw+zZs6lUqRIfffRRmssvWLCAypUrM3v2bOrUqcPIkSMZPnw4M2fOtHHlIiKSEfXvIiJQILcLSEtsbCy//vorEydOTHF/cHAwu3btSnOd3bt3ExwcnOK+kJAQFi1aRFxcHAULFky1TkxMDDExMebbN27cAODmzZuZrtWIj830stZ0J/pWbpcAZO21E9vRfnpfVvbR5GUNw7BWOfmO+nfbyOt9sV773JGXX3fIe/27XYbwq1evkpCQQLly5VLcX65cOS5fvpzmOpcvX05z+fj4eK5evUr58uVTrfP2228zderUVPdXqlQpB9XnjkFdFuZ2CSKPlFf301u3blGiRIncLsMhqH+3jRIl8uZ7zRHotc892Xntc7N/t8sQnsxkMqW4bRhGqvsetXxa9yd79dVXGT9+vPl2YmIif//9N6VLl85wO/bm5s2bVKpUiYiICIoXL57b5YikKS/up4ZhcOvWLSpUqJDbpTic/Ny/58X3gqPQa5977O21t4f+3S5DuLu7O87OzqmOily5ciXV0ZBkHh4eaS5foEABSpcuneY6Li4uuLi4pLivZMmS2S88lxUvXtwudmyRjOS1/VRHwC1L/ft9ee294Ej02ucee3rtc7t/t8sTMwsVKoSfnx8///xzivt//vlnmjVrluY6AQEBqZbfsGED/v7+aY4XFBER21P/LiKSxC5DOMD48eP573//y6effsrx48cZN24c58+fZ9SoUUDSV42DBw82Lz9q1CjOnTvH+PHjOX78OJ9++imLFi1iwoQJufUUREQkDerfRUTsdDgKQJ8+fbh27RrTpk3j0qVL1K9fn3Xr1uHl5QXApUuXUlxTtmrVqqxbt45x48bx4YcfUqFCBebOnUuPHj1y6ynYjIuLC5MnT0711auIPdF+Ksnye/+u90Lu0Wufe/Tap2YydO0tERERERGbstvhKCIiIiIijkohXERERETExhTCRURERERsTCFcRERERMTGFMJFJBWTyYTJZCI8PDzbbQwdOhSTycSUKVOsuo5IfhQWFobJZKJKlSoAhIeHm9+3+Ym1+ozWrVtjMpkIDQ21aLt5VWhoaIr9zRIy8xpXqVIFk8lEWFiYxbZrTxTCRcQqgoODeeGFF3j88cct2u6UKVMwmUwMHTrUou2KSN5jrX5GxBbs9jrhkjlxcXGaMU7sUv/+/enfv39ulyGSpxiGQWJiIs7OzrldSp6gfkbyMh0JtzPDhg3Dy8uLwoUL4+rqSsOGDVm2bBlw/whgr169GDRoEEWLFuWll14CkqZ8btmyJaVKlaJcuXL06tWLs2fPmtt9/fXXqV69Om5ubhQqVIhatWrx/vvv58pzlLxj8+bNNGjQgCJFitC2bVsiIyPNj508eZKePXvi6elJiRIlCAwMZOvWrebHH/6aOD4+npdeeomyZcvi4eHBnDlzzF81Pvx15N9//03v3r1xc3OjatWqrFq1ytzm1KlTAfjss88s/vWoyMMiIyMZMmQIlStXplixYvj6+rJ06VLOnj2Lk5MThQoV4u+//waSAnSlSpUwmUxs2LABeHTfnPweeOedd2jUqBEFCxbk8OHDfPDBB9SpU4dixYpRsGBBqlSpwmuvvUZiYmKuvA65berUqVSpUgUXFxdKly5N8+bN2bFjR6p+JvlzsmfPnowYMYJy5cpRoUIFPvzwQ3Nb0dHRDBs2jFKlSlGlShW+/PJL81Ce9IY9JCQkMH/+fLy9vXFzc6NatWpMmDCB27dv2+DZ297+/ftp3LgxRYoUISgoKNXQxN9//51u3bpRvnx5SpYsSbNmzVi/fr358eTfy4PfWKY3tOTcuXO0b9+ewoUL4+3tzebNm9Ot6969e7z11lvUqVMHNzc3atWqxfTp04mLi7PE07Y9Q+xKy5YtjT59+hijR482unfvbjg5ORnOzs7GwYMHjcmTJxuAYTKZjAYNGhjPPPOMMX/+fOP77783TCaTUaxYMaNPnz5G+/btDcCoWrWqcevWLcMwDKN3795G9+7djWeffdbo37+/4erqagDGDz/8kMvPWOwRYABGyZIljaFDhxplypQxAGPo0KGGYRhGRESEUbJkScNkMhkdOnQw+vTpYxQpUsQoWLCgceDAAcMwDGPIkCEGYEyePNkwDMOYPn26ARiFChUyBg8ebNSuXdtwcnIyAGPx4sUp1jGZTEanTp2MgIAAAzCKFStm3Lx501i2bJnRtGlTAzDq1KljvPDCC8bUqVNz4RWS/ODOnTtGrVq1DMBo1KiRMWTIEMPFxcUAjE8//dRo27atARgLFiwwDMMwtmzZYgBGpUqVjISEhEz1zV5eXgZgODk5Gd26dTMGDx5s/P7778bYsWONJ5980hg1apQxZMgQo2TJkgZgzJs3L8W2vLy8DMMwjLNnz5rft45m06ZNBmCUKVPGeOaZZ4z+/fsbtWvXNkJDQ1P1M8mfk4DRokULIyQkxAAMZ2dn4/Tp04ZhGMbTTz9tAEaJEiWMYcOGmX8HgLFlyxbDMAyjVatWKfqmF1980fy7HTJkiNGgQQMDMHr16pULr4h13bx509zn+/n5GQMGDDAKFixo3t/+/PNPw93d3QDMmcXJyckwmUzGpk2bDMO435cPGTLE3G7y6/zwa1ywYEGjb9++RsuWLQ3AcHV1Nc6dO5fmOj169DAAo1atWsbQoUON6tWrG4AxYcIEW75EFqMj4XZm5cqVtGnThpIlS1KpUiXKlClDQkJCir8MK1euzL59+1iwYAHPPvsss2fPxjAM6tevj4eHB3Xr1qVUqVKcPXuWtWvXArBw4UKefPJJ3N3dKVu2LJUqVQIwH60RScv8+fNZvHgxEydOBGDv3r0ALF68mOvXr1OxYkVq1aqFh4cH1atXJy4ujo8//jjNthYvXgzA5MmT+eyzzwgLC8PJKe0uqH379qxdu5YffvgBgFu3bvH777/Tv39/OnToAECTJk2YPXs2kyZNsuhzFkm2du1afv/9d0qWLMnOnTsJDQ1l8uTJALz33nsMHz4cgOXLl6f4f+jQoTg5OWWqb0728ssvs3r1aj777DNq1qzJ22+/zYABA/Dw8KBUqVJUr14dyJ99dkxMDJD02de9e3emT5/O8ePHGThwYLrr1K5dm61bt/Ljjz9SvHhxEhIS2L9/P4mJiXz22WdAUv/26aefsnr16gy3Hxsby/z58wF4/PHHKVmypHkM+sqVK/nzzz8t8TTtxpo1a/jrr78oVaoUO3bsYOnSpTz77LPmx5csWcLVq1epVasWYWFhrFixgpEjR2IYBrNmzcry9p599lm++OILwsLCqF27Nvfu3TOPAHhQRESE+VvRFi1aUKJECfz9/YGk32Ve/JZIY8LtyOnTp/H39+f69eupHrty5QouLi4ABAQEUKhQIfNj58+fB2D37t3s3r07xXonT54kKioKHx8f83IPtyuSnuQOrlSpUgDmr16T96WIiAjmzJmTYp2TJ0+m2daFCxcAqF+/PgDlypXD3d2dy5cvP3K7D25bxFaSh41UqVKFwoULA1C3bl3zYz169OC5555j+/btnDlzhlWrVmEymRg2bBjw6L75Qa1btzb/HBcXR7NmzThw4ECqmvJjnx0SEsLzzz/P559/TkhICABeXl4sWbIk3XV8fX3Nf+SXLFmSmzdvcvv2ba5evWoO9cl9Ub169TLc/l9//cXdu3eBpND9sJMnT1KuXLmsPzE7FRERASTt966ursD9/R7uvy/q1KljvhrPg++L9MTHx6d5f/K6JpOJ2rVrc+LECXMND3owwyxatCjFY3fu3OHixYtUrFgx4ydnZ3Qk3I788MMPXL9+nUqVKvHnn3+SmJhInTp1gKSxhsmSw3iy5J3u5ZdfxjAM87+LFy/y4osvsm3bNs6fP4+rqytnzpzBMAxzR/ZguyIPSz7p9+HLniXvc48//jiJiYnmfS46OpqlS5em2Vbyty8nTpwAksLE1atXs7RdwHzCWl486iF5S9WqVYGky/8lh7Djx4+bH3N1daVfv34YhsHw4cP5+++/adOmjXm9R/XND3qwXz927Jg5gO/evZvExESeeeYZIH/22fHx8cyZM4eoqCjOnTvH5MmTOXfunPn8kLQ8eMGCB/sRd3d3c7BM7ouOHTuW4fYfXGfdunUpfpenTp2iRYsW2X5u9ii5rw4PD+fevXtAytcoef9Ofv0g5fsCoGjRogDmg4p//fVXut8YJLdtGIa5neQaHvRgwD527FiK38Pp06fzXAAHHQm3Kx4eHgBcunSJcePGERkZyalTpx653vPPP09YWBjvvfcev//+O2XLluXkyZPs2LGDkydPmtu9d+8eL774InFxcRme+CDyKEOHDuW9997jf//7H82bN6dhw4ZERkaydetW5syZk+blA4cNG8arr77KlClTOH78OHv27MlWkK5cuTKQ9GH43HPP0ahRI0aMGJHTpySSyhNPPME//vEPTp48SYsWLWjQoAFffvklAOPGjQNg+PDhfPTRR+aTkpOHqMCj++b0Tip2d3fH2dmZhIQE3njjDdzd3fn++++t+2Tt2K5duxg0aBDNmjXD3d2dPXv2ACm/KcssJycnhgwZwscff8zo0aP56aefHvl56OLiwj//+U8++OAD+vbtS9euXQE4dOgQ169fz/Dob17UuXNn3N3duXr1KoGBgdSqVcu83wMMGjSIGTNmcOLECVq3bk358uX56quvMJlMjB07FgA/Pz8A1q9fz4svvkhYWFi6R8I/+ugjrly5wsWLF/n9999xdXVlwIABqZbz8vLiySefZM2aNbRp04bOnTtz9+5d9u3bR/ny5fPktcR1JNyO9OrVi9GjR1O0aFHWr19PYGAgzZo1e+R63bt3Z926dTRr1ozt27fzxRdfcOXKFZ577jnc3d1p2rQpb731FmXKlGHjxo1UqFCBXr162eAZiaPy8vJiz5499OjRg3PnzrF48WIOHz5Mly5d0r1e74QJE5gwYQJFixblxx9/5OmnnzZ/hVugQOaPB/Tq1YsnnniCu3fv8uGHH7JmzRqLPCeRhxUuXJhNmzYxcOBA/vzzT77++mtq167N4sWLzX/4+fv74+3tDUCJEiXo3r27ef1H9c3p8fT05OOPP8bT05Ndu3YRHx/PqFGjrPtk7VjFihWpW7cuW7duZeHChYSHh9O1a9dsX+HrvffeY+jQoSQkJLBx48YUR9TT64tmzZrFnDlz8PLyYtWqVaxZswZXV1dz6HQkxYsXZ926dfj5+XHkyBEiIiJ45ZVXzI+XK1eObdu28dRTT3HixAnWrVtHkyZNWLNmDe3atQNg4MCBDBs2DFdXV1atWkXXrl3NB1Ae9sYbb3D16lV++eUX6tevz/fff5/usitWrGDKlCmUKlWKZcuWsWHDBipUqMA///lPy78QNmAy8uN3WyJic/fu3cNkMpm/do+IiMDLywvDMNi+fbvDfaUrIvbp1q1buLm5mceM79y509z/RERE5MlhDZI3aTiKiNhEeHg4HTp0oEePHhQpUoQvvvgCwzBo3Lhxpr7xERGxhJ07dzJhwgQ6d+6MYRjmq6X06NFDAVxsSiFcRGyiZMmSVKxYkUWLFhETE4OXlxevvPIKEydOTPdShSIilubp6UnhwoX58MMPMQyDqlWrMnbsWIccWiL2TcNRRERERERsTIefRERERERsLFeHoyQkJJgvmi8iIpnj4uJivl66iFhGeHi4+TrXGiRgO/n5dc+1EH7z5k1OnTqV715wEZGcMplM1KhRg+LFi+d2KSIikk25EsITEhI4deoUxYoVo3z58jopS0QkkxITE7l06RKnTp2iYcOGOiIuIpJXGbkgOjra2Ldvn3Hr1q3c2LyISJ5269YtY9++fcY333xj/P7777ldjojVAQZgLFiwwKhdu7bh5uZmtG/f3jh79qx5mT/++MPo0aOHUaFCBaN48eJGixYtjLCwMPPjX375pdGgQQOjePHiRoECBYzy5csbo0aNMqKjow3DMIyzZ8+at2MYhhEREWFUq1bN4P+1d95RUV1dH/4hI2ABKwoqQbFQpjI0QaSJIhYMaggRjd2oUZdJNCq4jPpqJLGXKJbX3ivGElQQxAZqZGgqEQsWEBBBQEBg2N8fvHM/BoYiAcHkPGuxFveee/bZZ98z5+677ykAzZ07l4iINm/eTL169SJNTU1q1aoVWVhY0IkTJz6eIT4yzO71S4OGoFkEnMFgMD4cRd+ZnJyMoKAgPHr0qIE1YjA+Dr6+vrCzs0OnTp1w6dIleHh4oKSkBC9evIC1tTVOnjwJkUgEd3d33L17F/3794dMJgMAJCUloVOnTvD29sbYsWMBAAEBAfjpp58qlJOcnAwXFxc8fvwYixYtwq+//orHjx9j+vTpSE5OxpgxYzBs2DDI5XLExcV9TBM0CMzu9QNbJ5zBYDA+UTp06IC//voLSUlJ6N69e0Orw2DUOzt27ICnpyeSk5PRpUsXxMbGIjIyEsHBwcjKyoKBgQGMjY0BAN27d0dsbCy2bt2KLVu24LvvvoOJiQmio6Px5s0bmJiYICUlBRcvXsTKlSuVynFxccHDhw/h7+/PbdmuWEiiffv2GDZsGExMTNC9e3eUlJR8XCM0AMzu9UOjcsLVHWbUi1x5+KZ6kctgMBgNiZqaGng8HltlivGvwczMDADQqVMntGrVCllZWXj+/DmePXsGoHTb+fXr1yvlefjwIQDA09MTZ8+erSAzLS2twrmEhAS0b98ePj4+3DlTU1MsX74c69evx9ChQwGUvghv2rQJX3zxRd1UsJHC7F4/sPEgZTh9+jRMTU0hkUgQGxtbJzLDwsJw8eJF7jg5ORnOzs51IrsmREREQCgUwtzcHBcuXPho5daG5ORkWFtbN8q327CwMFhaWlaaXlRUhK5du6KwsLBG8saNG4dNm/5dL4dlbVhQUAALCwu8ffu2gbViMBifEvfu3QNQ+rxQ9B8GBgbcdvO9e/dGSUkJiAhEhHfv3mH//v3IysriHMEjR45ALpdjxYoVAFQvizdx4kS8fv0arq6unLMol8sxb948pKamIjk5GVu3bkVaWhrmz59f7/VuaJjd64dGFQlvaAICArB06VKVb1bFxcXg8T7cXGFhYcjNzcWAAQMAlL5FhoaG/m1da8qePXvw9ddfY+7cuR+tzNqybNkyzJgxo07nCtT2vn0oYWFhsLW1hYaGRr2X9bGpDxtqaWnBx8cHa9euxeLFi+tUNoPB+OcyZcoUnDt3DuHh4SAiCAQCWFtbo1OnTli9ejUiIiLQp08fiMVivHz5EleuXMH69evh4+MDbW1t5OTkYPXq1Th37hwCAwMrLWfr1q14/fo1Tp8+jf79+yMsLAxv376FlZUVHBwcoKenxzmmbdq0+Ui1bziY3esHFgn/H7NmzcLVq1cxb9482NnZASj91Lt69Wo4OTlhwYIFiI2NRd++fSGVSmFmZsa9zQHA27dvMWnSJAiFQojFYkyYMAEymQwBAQHYu3cvJBIJli5diqdPn6J9+/ZcvqCgIEilUohEIjg6OnKNKywsDBKJBNOnT4dYLAafz8edO3dU6p6amgpPT08IhUIIBAJs27YNAODv748jR45g/fr1kEgkyMrKUsq3e/duDBgwAN7e3jAxMYGLiwvi4+MxePBg9OrVC97e3lxUOicnB5MnT4a1tTVEIhGmTp2KoqIiAMCaNWtgZWUFc3NzWFtbIzIykitDTU0Nv/zyC2xsbNCtWzfs2rVLZR0KCgpw5MgRjBw5skZ579y5A1tbW4hEIlhbW+P69esAwNl36dKl6Nu3LzZu3AgnJyfMnTsXDg4OMDAwwMqVK3H48GHY2dnB0NAQhw8f5uSOHj0alpaWEIlEGDJkiMrPZaoIDAyEp6dnhfMRERGwsLCARCKBQCDAli1bKlxz9epV7v4WFRVh/vz5sLa2hkQigbe3N7KyspCbm4u2bdtyNpdKpdznusePH1c6Hrgx2/Crr77Cjh07amRfBoPBAIBff/0VN2/exMuXL9GvXz+cPn0a6urqMDQ0RGRkJEaMGIGkpCTs2rULsbGx8PDwQO/evdG0aVMcPHgQPXr0QHR0NF6+fIkff/yx0nLU1dVx+PBhODo6IiYmBu7u7lBTU4OdnR0iIyOxY8cOyGQy9O/fv9Ln2j8JZvd6oiGWZFEsUahYnkZBk77f1stfTXF0dKQzZ85wxwBo+fLl3HF2djYVFBQQEVFeXh5JJBK6ffs2ERGNGzeOZsyYQXK5nIiI0tLSiIjop59+oh9++IGT8eTJE2rXrh0REaWmplK7du0oJiaGiIj2799PfD6fiIhCQ0OJx+Nx8rds2UIDBgxQqbeXlxfNnz+fk9mlSxeKjIwkIqKxY8fSxo0bVebbtWsXtW7dmp4/f05ERIMHDyaRSERZWVlUVFREQqGQLly4QEREkydPpr179xIRUUlJCU2cOJHWrFmjVFciops3b3J1UNhw3bp1RER07949atmyJRUVFVXQJTw8nKytrZXOVZb3/fv3ZGBgQEFBQUREdPXqVdLT06Pc3FxuqaMDBw5wchwdHcnLy4vkcjm9fPmStLS0yM/Pj4iIIiMjSV9fn7s2PT2d+3/FihX07bel7Sc0NJQsLCxU2rGkpISMjIwoOzu7QpqHh4eSLm/evCGi/78vBw8eJIlEwi33tHz5cvrPf/7DXb906VKaNWsWERHZ29vTlStX6PXr1yQSiahnz55ERBQQEEBTpkxRqVtjt2G3bt3o/v37KnVnVI6iDz1+/Dj98ssvdP78+YZWicGoV/C/JezKLo3HqH+Y3esXNhylGiZMmMD9n5+fj+nTp0Mmk6FJkyZ4/vw5ZDIZLC0tcfbsWfz555/cUApdXd1qZUdGRkIikUAoFAIAfHx88O233yIlJQUAYGxszI2htbW1xapVq1TKCQ4ORnR0NIDSyQrDhw9HSEgIrK2tq9WhT58+3Jguc3NzdO3aFa1atQIAiMViPH78GEBppDciIgKrV6/mbKEYehEVFYXly5cjIyMDPB4P9+7dQ2FhIZeuiNiampqCx+Ph1atXXJkKXrx4AT09vQr6qcqbmZkJDQ0NuLm5AQDs7e3RoUMHxMTEQF9fH1paWvjqq6+U5HzxxRdo0qQJOnXqhPbt2+Pzzz8HAFhYWCAlJQUFBQXQ0tLCgQMHsG/fPrx//x75+fkqdSrP7du3YWJiAm1t7Qppzs7OWLZsGRITE+Hi4gJ7e3subdeuXWjatClCQ0PRunVrzs7Z2dk4fvw4AKCwsJCLcru6uiI4OBipqalwc3PD/fv3ERcXh+DgYHh5eVWqX2O2oZ6eHl68eAETE5NqrMxgMBgMxj8L5oRXQ8uWLbn/fX190bFjR0RFRYHH42H48OEoKCiotWwigpqaWoXzinNaWlrcOXV1dRQXF1cqq7wcVXJVUb6MysokIgQGBsLIyEgpf2FhIUaMGIGwsDBYWFggOzsbrVq1UnLCa1KP5s2bIz8/v1r9iouLq7VbixYtKqRXVk/FboPFxcW4du0aNm3ahBs3bkBXVxe///47li5dWqGc8pw6dYpzSMsze/ZseHh4ICQkBL6+vhAIBNi8eTOA0pec8PBwxMXFcc45EWHz5s1wcXGpIMvV1RVz585FWloahg8fjs6dO+PSpUu4cuUKAgICKtWvMduwoKAAzZo1qzSdwWAwGIx/Ko3KCW/sSwlmZmZCIBCAx+MhISEBly5d4pwlDw8PrFy5EuvXr0eTJk2Qnp4OXV1d6Ojo4OXLlyrl2draYuLEibh//z5MTU1x+PBhdOnSBXp6enjw4EGN9XJ1dcW2bduwZMkSpKen49SpU1wkta7w8PCAv78/Nm/eDB6Ph8zMTGRkZKBDhw4oKiqCgYEBAGDjxo21ki8Wi2tcZxMTE7x//x6XL1+Gi4sLbty4gbS0NAiFQqSnp9eqfKD0/uro6KBt27YoLCzE1q1ba5Tv999/x+XLl1WmJSQkwNjYGEZGRjAwMICvry+XJpVKMWfOHAwbNgybN29G//794eHhgTVr1qB3795o3rw58vLy8OTJE/D5fNjY2ODBgwdIS0vD2rVr0blzZwwePBgGBgZo167dB9W1MdhQLpfj8ePHEAgEtS6PwWD8OyAVK2kw6h9m9/qFTcz8ABYuXIgdO3bAysoKCxcuVIpWrl27Fnl5eRAIBJBIJJyz5enpiTt37nATM8uiq6uLffv2wcfHB2KxGFu2bMHRo0c/WK8NGzYgJiYGIpEIzs7O8PPzq9FQlA9h3bp14PF4kEgkEIlEcHV1xdOnT6Gjo4OlS5fC2toaDg4O0NTUrJX8rl27QldXF/Hx8dVeq6GhgRMnTsDPzw8ikQizZ8/GsWPH0KJFi1qVrcDd3R09evSAiYkJ3NzcIJFIqs2TkJCANm3aoGPHjirTN27cCD6fD3NzcyxcuJAbzqPAzMwMQUFBmDlzJgIDAzF//nxIJBLY2NhAJBKhd+/e3K5jPB4P9vb2+Oyzz9CsWTPw+XwUFRXB1dX1g+vaGGx47do12NjYcMOfGAwGg8H4N6FGDfCak5eXx0V/mzdv/rGLZzRSjh49iitXruC3335raFVqjL+/P3g8HubMmdPQqnxyeHt7Y9KkSbV6ifi3o+hDnz59ikePHkEoFMLd3b2h1WIwGAzGB9CohqMw/t14eXkhNTUVJSUldbpWeH3yT9gsoCEoKCiAk5MTc8AZKpHL5WwXUAaD8cmiqanJzZeqCuaEMxoVM2fObGgVGB8BLS0tTJ06taHVYDRCsrOzkZiYyMaiMhiMTxY1NTX06NEDOjo6VV7HnHAGg8FgNArkcjkSExOhra0NfX39T+aLGIPBYCgoKSlBSkoKEhMTIRaLq4yIMyecwWAwGI2C9+/fg4igr6+vtDwsg8FgfEro6+sjOzsbd+/ehVQqrdQRZ2EGBoPBYDQqWAScwWB8yij6sFu3biEsLAwlJSUqr2tUkfCTl1/Ui9zhLl2qv4jBYDAYDAaDwagjWrRogfj4eFhZWakcH87CDWU4ffo0TE1NIZFIEBsbWycyw8LCcPHiRe44OTkZzs7OdSK7JkREREAoFMLc3BwXLlyoUZ7AwEDcunWrznXZvXs3Ro4c+bflhIWFwdLSskbXOjk54ezZs7Uu686dO9y271Uhk8kqrPEukUhU7gL6d5DJZBg8eHCdyqwrqru/aWlpH7Qxz9+9d58iZW346tUr2NjYVLlT7r8BdYcZ9fLHYDAY9Y2WlhaKi4srXe2JOeFlCAgIwNKlSyGTySAUCpXSavsgLO+Ed+rUCaGhoX9Lzw9hz549+PrrrxEVFQU3N7ca5akvJ7w2NLQDYmlpiQMHDlR7nSonXCaT1fmW7L6+vpg3b16dyvxYNj59+jSGDh36Ucr62NSHDfX09GBjY4P9+/fXuWxGzWHBmfqHiNCnTx8kJSUBAMaNGwc1NTVERUVx1+Tm5qJly5Y1DsCoqakhNzcXQOlmc2lpaVxaQEAA1q5dyx3fv38fgwcPRvfu3dG9e3e4u7vX2b2ujJo8ZysLXOzevRutW7eGRCKBQCCAu7s7nj17VmtdFi9eXOleF127doWJiQkkEgnMzMxqvY/H1q1bOTkZGRm11rWxEB0djUGDBv1tOcwJ/x+zZs3C1atXMW/ePNjZ2QEo/RGvXr0aTk5OWLBgAWJjY9G3b19IpVKYmZlhxYoVXP63b99i0qRJEAqFEIvFmDBhAmQyGQICArB3715ux8ynT5+iffv2XL6goCBIpVKIRCI4Ojri3r17AEo7aYlEgunTp0MsFoPP5+POnTsqdU9NTYWnpyeEQiEEAgG2bdsGoHQjmSNHjmD9+vWQSCTIyspSyhcREQELCwvuh7xlyxacP38ev//+O/z9/SGRSLBjxw68evUKzs7OsLCwAJ/Px6xZs7jlwxYvXoxRo0Zh6NChMDMzg4uLC968eQMAKCwsxDfffINevXrB2dkZkZGRXNlV2XLcuHGYNWsWBg4cCLFYDKB0t9IePXrA0dGxyujovXv3YGNjA6lUCh8fHxQUFHBpr169gpeXF6ytrSESibBo0SIAwP79+5WcQyJCt27dEBMToxR1Ly4uhpubGywtLcHn8+Hj44O8vDykpaVh0aJFCA4OhkQi4ZbeK/sQuHPnDmxtbSESiWBtbY3r168DANceFi1aBAsLC/To0QPnz59XWbdnz54hPj4eDg4ONcpbXduaNWsWbG1tcerUKXTt2hWLFi2CnZ0dPvvsM+zfvx/r16+HtbU1unfvjrCwsCptUBMCAwPh6elZ4fyZM2cgEom4dnj69OkK1xw/fhwSiQSPHj1CTk4OJk+ezN3HqVOnoqioCAkJCTA2NubuYfv27eHn5wcACAkJQb9+/SrIbew2HDVqFLZv314j+zLqBxacUU1dvngeO3YMxsbGMDQ05M5ZWFhg586d3PHhw4dhampaK/nlnfCpU6fiu+++A1D6AuTo6AgfHx88evQIjx49wtdffw0nJyc8f/68ljWqnpo44VUFLlxdXSGTyRAXFwcTExOuPvXB8ePHIZPJcOHCBfj5+SEmJqbGeRXtZN26ddi3bx9kMhnatWtX4/xyufyD9f0YiMVi8Hg8rl+vNdQAvHv3ju7cuUPv3r1TOn8i5Hm9/NUUR0dHOnPmDHcMgJYvX84dZ2dnU0FBARER5eXlkUQiodu3bxMR0bhx42jGjBkkl8uJiCgtLY2IiH766Sf64YcfOBlPnjyhdu3aERFRamoqtWvXjmJiYoiIaP/+/cTn84mIKDQ0lHg8Hid/y5YtNGDAAJV6e3l50fz58zmZXbp0ocjISCIiGjt2LG3cuFFlPg8PDzpw4AB3/ObNG5V58vPzKScnh4iIiouLafDgwXTs2DGufkZGRpSRkUFERF9++SX9/PPPRES0YcMG6t+/PxUWFtK7d+/IwsKCRowYUa0tx44dS+bm5lyZv//+OwmFQsrJyaHi4mIaOnQoWVhYqKyTVCql3bt3ExHRzZs3qUmTJtw9HTBgAF25coWIiIqKisjNzY1OnjxJeXl51K5dO0pJSSEiosuXL5NUKuXug6KskpISev36Nff/1KlTaeXKlUREtGvXLq5uCgBQTk4OvX//ngwMDCgoKIiIiK5evUp6enqUm5tLT548IQAUGBhIRER//PEH9erVS2Xd9u7dS15eXtxxVXmra1tqamp09epVTpahoSHNmTOHiIhu3bpFzZo1o99++42IiI4cOUK2tra1soGC7OxsMjIyopKSkgppIpGIrl+/TkREcrmcMjMziej/f4+rVq0iBwcHro1NnjyZ9u7dy+kwceJEWrNmDRERGRgYUFJSEt29e5dsbW2pd+/eREQ0f/58rl2WpbHbsLCwkJo1a0a5ublKeiv60OPHj9Mvv/xC58+fV2n3Tw1Vz4Ymfb+tl7+aMHPmTGrRogV169aNu38AaNWqVeTo6Ehz5syhmJgYsre3J3NzczI1NVVqZ1lZWTRx4kQSCAQkEolo/PjxFBUVRR07diRdXV0Si8W0ZMkSpecCUWk7NDc3J6FQSA4ODhQfH09Epe1OLBbTtGnTSCQSkZmZGddvlufVq1f0+eefk0AgID6fT1u3biUiohUrVlCbNm2oc+fOJBaLud+bgqrqM3bsWJo5cya5ubmRmZkZEZX2S9bW1mRubk4ODg4UGxtbrZzyuLi4UHBwsFI5/v7+ZGhoyD0n7OzsaPPmzVx/XN5mOTk5VNalUfS/S5YsoaZNm5KxsTGJxWKKiopSei77+fnRF198UUGnL774gubNm8fpU/aZ+MMPP9BPP/1ERETBwcHUu3dvkkgkxOfzaefOndx1jo6O9OOPP5K9vT0ZGRnRN998Q0RE586dU7oH27dvV2mXQYMGcc/yspTvJ86dO8f1TUFBQdSnTx+SSqVkbW3NPfNSUlLIycmJpFIpmZmZ0cyZM7n+uKw94uPjSSAQcH2KoaEhd0+JiKysrOjYsWOUnZ1NkyZNIisrKxIKhfTNN99QYWEhV29fX19ycXEhFxcXGjFiBHcPFHrv3buXBAIBCYVCGjRoEL148YKrm5ubG40ZM4YsLCzoxo0bBIB+/vlnsrKyom7dutGlS5do/vz5JJFIyMzMjOLi4mpUx6+++oqGDBlCpqam5OzszD1TiIj8/f2536mNjQ3XB1XWvomIDh06RKNGjVJ57xR92aFDh2j16tWcT1ieRjUxszEyYcIE7v/8/HxMnz4dMpkMTZo0wfPnzyGTyWBpaYmzZ8/izz//5GbE6urqVis7MjISEomEi674+Pjg22+/RUpKCgDA2NiYi8La2tpi1apVKuUEBwcjOjoaANChQwcMHz4cISEhsLa2rrJ8Z2dnLFu2DImJiXBxcYG9vb3K60pKSjBv3jxcu3YNRIS0tDRIJBJu7Kq7uzvatm3L6an4jBcaGoqxY8eiadOmaNq0KUaPHo1r165Va0ugdPdMxRJloaGh+PLLL7njCRMmYNmyZRX0zM7ORlxcHMaMGQMA6N27N2fbd+/e4fLly0hNTeWuz83NxYMHD+Dp6YkRI0Zg//79mDNnDnbt2oXx48dXkE9EWLt2Lc6dO4fi4mK8ffuWi0pXRUJCAjQ0NLiIk729PTp06ICYmBjo6+ujRYsWGDZsGGe/R48eqZTz4sUL6OnpKZ2rLG91batXr14V7veXX34JAJBKpcjPz4eXlxeA0ojU48eP/5YN/vjjDwwcOBBqamoV0vr164fZs2dj5MiRGDBgACQSCZe2ePFidOrUCRcvXoSmpiaA0ghSREQEVq9eDaC0LWloaHCygoODkZGRgdGjR2Pbtm14+/YtgoODsXnzZpW6NWYbNm3aFK1bt0ZKSgp69OhRlYkZ9cCGDRsQExODOXPmYMiQIdz59+/fcxGwnJwcBAcHQ1NTE/n5+bCzs0P//v1haWmJ2bNno2XLloiOjkaTJk2Qnp4OXV1dTJ06Fbm5uVyf/vTpU052WloaRo8ejdDQUAiFQhw4cABeXl6Ii4sDAMTHx2PHjh3YvHkzAgIC4Ofnp3JIyaxZs2BiYoJTp04hLS2N++o5f/58PHjwAJaWlpgxo+LY+K5du1ZaHwC4du0awsPD0bJlS1y/fh2HDx9GeHg4NDU1cfXqVfj4+CA6OrpaOQqKiopw48YN2NraKp1v0aIFXF1dERgYCLFYDCKqVSR80aJF2LlzJ44fP84N7QgMDOTS7969i/79+1fIZ2trq/S1ojKkUimuXbsGdXV1vHnzBlKpFAMHDoS+vj4A4NGjRwgLC0NhYSHMzMxw8+ZNDBo0CB4eHpXeA6C0XT148ABWVlZVli+Xy3Hs2DGuj1myZAmCgoKgo6ODxMREODo64unTp2jdujXOnDmDli1bQi6XY9iwYThx4oTSPJ7Lly9j5syZOHDggFI/rCA2NhYPHjyAWCzGDz/8AAcHB2zfvh1EhMmTJ2PTpk1cRF4mkyEoKAhNmzYFUNquFPcgLi4Oc+fOxZ9//onOnTtj+fLlmDJlCs6dOwegtI1FRUWhZ8+eXNk6Ojq4desWjh07hmHDhuHo0aNYsWIFfv31VyxfvhwHDx6sto6RkZG4ffs22rZtC29vb2zduhULFizAnj17EBgYiOvXr0NHRweZmZnQ1NSssn0DgJ2dHWbPnl1tG6kK5oRXQ9m1an19fdGxY0dERUWBx+Nh+PDhSsMdPhQiUumUKM5paWlx59TV1av8/Fdejiq55Zk9ezY8PDwQEhICX19fCAQClY7KmjVrkJGRgcjISGhpaeH7779XqndlelIVO95VZ8uydq9KTnkqq3dJSQnU1NRw+/ZtrlMoy/jx4zFp0iRMmTIFZ8+exbp16ypcc/DgQVy5cgXh4eHQ1tbGhg0bEB4eXq1OH3qfK/v81rx58woTPSvLW12ZqtZgVshSrGda9lhxT2trg1OnTim90JZlzZo1iI+P517afHx88OOPPwIofRBeuHABT548gYmJCVe3wMBAGBkZVZDl6uqKc+fO4c2bN9iwYQMePnyIkydP4smTJ7CwsFBZfmO3YUFBQZ3PLWD8Pf7JwZkPCZCcPn0a0dHRsLGx4fKnp6ejsLCwWjkKXr9+DQ0NDTRv3ryCLhMmTMDixYshFotVBkbqClW/cyKq0e8uIyMDEydOxF9//QUej4fXr18jPj6ec8K9vb2hrq6OZs2acUPqyr9wqKKqwAUAbvgjUPoisHr1ahw9ehSJiYkVXuqfP38OPT29KoNply5dwvnz53HhwgUYGBgo5R85ciS0tLTQvHlz7Ny5Ez179qwyGAIAY8aMUfmsBUoDa0OGDEHnzp0BANOnT8eyZcu4Z729vb2SAw4oBziaNGnCLVBgYWGBkydPAqh9wPDs2bOYNm0at3pJmzZtAFTdvjU0NKCnp4fU1FQUFRVVWtfqYGPCP4DMzEx06dIFPB4PCQkJuHTpEpfm4eGBlStXcmtBpqenAyh9e3v79q1Keba2tpDJZLh//z6A0jFvXbp0qRDtrA5XV1duHHh6ejpOnToFFxeXavMlJCTAyMgIkydPhq+vLyIiIlTqnJmZCT09PWhpaSE1NRXHjh2rkV79+vXDvn37UFxcjPz8fBw8eFBJZmW2VCXn6NGjePfuHeRyOXbv3q3yOh0dHQgEAm4i5a1bt7gfmba2Nvr27Qt/f3/u+uTkZLx4UbosZu/evVFSUoIff/wR/fv3536oZcnMzES7du2gra2NnJwcJT2qus8mJiZ4//49Ll++DAC4ceMG0tLSKowvrQ6RSIQHDx7U6Nq6alvlqcoGlVFYWIibN2/CyclJZfqDBw/A5/MxY8YMTJs2jWuHAODm5oYdO3ZgyJAhkMlkAEp/a/7+/pxTm5mZicTERAClv4WQkBAkJSWhV69ecHV1xZIlS+Do6PjBa083BhumpqaCx+OhU6dOf6tMRt1SWXAmOjoaTk5On3Rwprr6lA+QKOY/Kf6Sk5OhoaFRY7s0b968UnvZ2dnh5cuXOHz4MLy9vZXSeDyeUsCitjaXSqW4ceNGhfMRERHc/LCqypo6dSocHR0RGxsLmUyGXr161ShIVR2nTp3C559/Xmm6Yky4TCbDzp070b59exARBg4cqHQ/Xr58CSMjI6VgWkxMDEaNGqWkZ8+ePUFEKsepK8aE37hxg3NoFcEQRTkJCQlKQbyqNtsq38bLt8vqAhyKr6KKY4VNq6vjh96Lqto3UNoOFF/6a0ujcsKHu3Spl7+6YuHChdixYwesrKywcOFCJUd37dq1yMvLg0AggEQiga+vLwDA09MTd+7c4SZmlkVXVxf79u2Dj48PxGIxtmzZUmGFjZqg+GQqEong7OwMPz+/aqMdALBx40bw+XyYm5tj4cKF3BvtmDFjcPDgQW5i5qxZs3Djxg1IJBJMmDABrq6uNdJrypQp+Oyzz2BmZobBgwejb9++XFpVtizPkCFDMGTIEIjFYri4uEAkElV67d69e7Fp0yZIpVJs27ZN6Q32wIEDuH//PoRCIYRCIUaMGKE0S3v8+PHYunVrpRGXr7/+Grm5uTAzM8Pw4cOV6tOvXz+8e/cOYrGYm5ipQENDAydOnICfnx9EIhFmz56NY8eOoUWLFpUbTwX29vZ49uwZN/G1KuqqbZWnKhtUxuXLl2Fvb19pR7VgwQKuHe7btw+LFy9WSndwcMChQ4cwYsQI3Lx5E+vWrQOPx4NEIoFIJIKrqyv3Ob9jx47o2LEjF2lydHREcnJyjdtsWRqDDYOCguDp6Vkj54nRMPzTgjMfEiAZOnQo9u7dy01gLCkp4RYQqKmcVq1aQV9fnxuuVZ4NGzZg1apV0NbWVjqvp6eH4uJiJCQkACjt+yujKntPmzYNoaGhSkGiQ4cO4d69e5gyZQoAoHv37tzCAhkZGUqTtzMzM2FoaAg1NTWEh4dzXx+qoyqdqgtcVMaAAQMQFBTEDV0CwDnV1QXTunbtipCQECxZsqRKWyqoKhhSHf369cP58+fx6tUrAKWTn/v16/e3+7naBgw9PDywZcsWZGdnAwCysrIgl8urbN9A6ao6VfkjNULlSPF6prKJmQwGo2p++eUX+vXXXxtajQ/im2++4SbyMj4Me3t7un//foXz/6aJmQ2Nqgn7iknjRER3794lPp9PlpaWNHLkSPLw8OAm8b19+5YmTJhApqamJBaLadKkSURE9PjxY5JIJNVOzBSJRBUmZpadlB4bG0uGhoYq9VZMzBQKhcTn8ykgIIBLq2rCflX1UZXvwIEDnK6mpqbc5OSq5JTn+++/p02bNlWrX/n679y5k7p160Z9+/aln3/+WeXETCKi7du3U8+ePVVOzCQiiouLI3d3dzIyMqIOHTqQsbExZWVlcekZGRnk7OxMAoGAhg0bRt7e3tzEzIsXL1KPHj3IxsaGxo8fT7a2tlx7Kd92RowYQbt27SKi0snbZmZmKidm/vHHH+Tj46PSVkRVT4K/cOEC2djYkEgkIhMTE27i4NOnT8nKyorEYjENGjSIxo0bx9mgrD3S09PJwsKCm1RefmKmguzsbJo2bRrx+XwSCoUklUrp0qVLKuutSs6ePXtIKBSqnJhZ2SIHRBUn5JZtEzWtIxHRxo0baezYsdyxv78/dz9sbW25Pqiy9k1EtHjx4konHNd0Yqba/yr4UcnLy8P9+/dhamqqchwYg8FQzfv37/Hf//4X06dPb2hVGPVMamoqQkJCMGrUqAppij706dOnePToEYRCIdzd3RtAy7qFPRv+nSQlJWHkyJGIjIz84GFjdc3z588xbNgwDBkypMLX64/F1KlT4erqWieb2zHqh8LCQlhZWSEkJERp2WkFir7s4cOHSE5OxpgxY1TOCWETMxmMTwhNTU3mgP9L6Nixo0oHnMH4p2FoaIh58+YhOTkZXbrU3RDS2mBgYIC7d+82qA4BAQENWj6jep48eYIVK1aodMA/hAZ1whXj5BgMBoNRc1jfyfinwaK+jE8JY2NjbnO4v0ODOOGamppQU1NDSkoK9PX1G/zzE4PBYHwqlJSUIDk5GUSEoqKihlanXmAvGQwG41NG0YdVN+K7QZxwdXV19OjRA4mJidxsVAaDwWDUDCLCixcv/nHOKgvQMBiMT52ygZLCwsIqr22w4Sg6OjoQi8W4efMmoqKioK+vz5bhYjAYjBpQVFSEkpISLhpedv3bTxkWoGEwGP8EFIGSyjbfU9CgY8LV1dVhaGiImJgY/PXXX2xXOAaDwfgA8vPz0bp1a3Tv3r2hVakzdHR0wOfzcfz4cWRnZ0NXV5cFaBgMxidFUVER5HI53rx5g5YtW1a6L0iDLFFYnqSkJMTGxiInJ6ehVWEwGIxPBm1tbfD5fHTr1q2hValzkpKScP78+Uo3NGEwGIzGDBFBW1sbAwYMQK9evVRe0yiccAaDwWAwypOeno6MjIxqP+kyGAxGY0NdXR2tW7eucrdb5oQzGAwGg8FgMBgfGTb1nMFgMBgMBoPB+MgwJ5zBYDAYDAaDwfjIMCecwWAwGAwGg8H4yDAnnMFgMBgMBoPB+MgwJ5zBYDAYDAaDwfjI/B9EauFqDrYl1QAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(figsize=(7.5, 4.8))\n", + "ax1= fig.add_subplot(2,1,1)\n", + "ax2= fig.add_subplot(2,2,3)\n", + "ax3= fig.add_subplot(2,2,4)\n", + "# moving the lower graphs down to create space for the legend\n", + "pos1 = ax1.get_position()\n", + "ax1.set_position([pos1.x0, pos1.y0-0.08, pos1.width, pos1.height])\n", + "\n", + "# add a - c to the graphs\n", + "ax1.text(-0.62, 1.2, \"A)\", fontsize=14, fontweight=\"bold\")\n", + "ax2.text(-0.5, 1.2, \"B)\", fontsize=14, fontweight=\"bold\")\n", + "ax3.text(-0.9, 1.2, \"C)\", fontsize=14, fontweight=\"bold\")\n", + "\n", + "# graph 1\n", + "params = ['mean', 'std', 'area', 'height']\n", + "x = np.arange(len(params)) # the label locations\n", + "width = 0.12 # the width of the bars\n", + "colors = [\n", + " (\"#023d6b\"), # Jülich dark blue\n", + " (\"#adbde3\"), # Jülich light blue\n", + " (\"#af82b9\"), # Jülich hyancith violet\n", + " (\"#eb5f73\"), # Jülich raspberry\n", + " (\"#fab45a\"), # Jülich apricot\n", + " (\"#faeb5a\"), # Jülich lemon\n", + "]\n", + "multiplier = 0\n", + "for metric, result in all_data.items():\n", + " offset = width * multiplier\n", + " if metric in [\"skew normal data, normal model\", \"normal data, skew normal model\"]:\n", + " rects = ax1.bar(x + offset, result[0], width, label=metric, yerr=result[1], color=colors[multiplier], hatch=\"///\")\n", + " else:\n", + " rects = ax1.bar(x + offset, result[0], width, label=metric, yerr=result[1], color=colors[multiplier])\n", + " multiplier += 1 \n", + "ax1.set_ylabel(r\"$\\bf{F_{y / \\^y}}}$ (-)\", fontsize=9, fontweight=\"bold\")\n", + "box = ax1.get_position()\n", + "# legend below\n", + "h, l = ax1.get_legend_handles_labels()\n", + "ph = [ax1.plot([],marker=\"\", ls=\"\")[0]]*2\n", + "handles = ph[:1] + h[:3] + ph[1:] + h[3:]\n", + "# labels = [\"normal model:\"] + l[:3] + [\"skew normal model:\"] + l[3:]\n", + "labels = [\"normal model:\", \"normal data\", \"normal data (higher noise)\", \"skew normal data\",\"skew normal model:\", \"skew normal data\", \"skew normal data (higher noise)\", \"normal data\"]\n", + "order = []\n", + "leg = plt.legend(handles, labels, ncol=2)\n", + "for vpack in leg._legend_handle_box.get_children():\n", + " for hpack in vpack.get_children()[:1]:\n", + " hpack.get_children()[0].set_width(0)\n", + "\n", + "ax1.set_position([box.x0, box.y0 + box.height * 0.15, box.width, box.height * 0.85])\n", + "ax1.legend(handles, labels, loc=\"upper center\", bbox_to_anchor=(0.5, -0.32), fancybox=True, shadow=True, ncol=2, fontsize=8)\n", + "# legend on the right\n", + "# ax.legend(loc=\"center left\", bbox_to_anchor=(1, 0.5), fancybox=True, shadow=True, ncol=1)\n", + "params = [x if not x == \"baseline_intercept\" else \"baseline\\nintercept\" for x in params]\n", + "params = [x if not x == \"baseline_slope\" else \"baseline\\nslope\" for x in params]\n", + "params = [x if not x == \"std\" else \"standard\\ndeviation\" for x in params]\n", + "ax1.set_xticks(x + 2.5 * width, params, fontsize=9, fontweight=\"bold\")\n", + "\n", + "# graph 2\n", + "params = [\"area\", \"height\"]\n", + "x = np.arange(len(params)) # the label locations\n", + "width = 0.2 # the width of the bars\n", + "multiplier = 0\n", + "for metric, result in comparison_dict.items():\n", + " offset = width * multiplier\n", + " rects = ax2.bar(x + offset, result[0], width, label=metric, yerr=result[1], color=colors[multiplier])\n", + " multiplier += 1 \n", + "ax2.set_ylabel(r\"$\\bf{F_{n / sn}}}$ (-)\", fontsize=9, fontweight=\"bold\")\n", + "box = ax2.get_position()\n", + "# legend below\n", + "# ax2.set_position([box.x0, box.y0 + box.height * 0.15, box.width, box.height * 0.85])\n", + "ax2.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2), fancybox=True, shadow=True, ncol=1, fontsize=8)\n", + "# legend on the right\n", + "# ax2.legend(loc=\"center left\", bbox_to_anchor=(1, 0.5), fancybox=True, shadow=True, ncol=1)\n", + "params = [x if not x == \"baseline_intercept\" else \"baseline\\nintercept\" for x in params]\n", + "params = [x if not x == \"baseline_slope\" else \"baseline\\nslope\" for x in params]\n", + "params = [x if not x == \"std\" else \"standard\\ndeviation\" for x in params]\n", + "ax2.set_xticks(x + 0.5 * width, params, fontsize=9, fontweight=\"bold\")\n", + "ax2.set_ylim(0.0, 1.264356410280364)\n", + "\n", + "# graph 3\n", + "categories = [\"overall\", \"single\\npeaks\", \"double\\npeaks\"]\n", + "means = [\n", + " df_comparison_total[\"area MQ / PP\"].mean(), \n", + " df_comparison_single[\"area MQ / PP\"].mean(),\n", + " df_comparison_double[\"area MQ / PP\"].mean(),\n", + "]\n", + "sds = [\n", + " df_comparison_total[\"area MQ / PP\"].std(),\n", + " df_comparison_single[\"area MQ / PP\"].std(),\n", + " df_comparison_double[\"area MQ / PP\"].std(),\n", + "]\n", + "ax3.bar(\n", + " x=categories,\n", + " width=0.5,\n", + " height=means,\n", + " yerr=sds,\n", + " color=(\"#023d6b\"), # Jülich dark blue\n", + " label=\"fraction of area (MultiQuant / PeakPerformance)\",\n", + ")\n", + "ax3.set_ylabel(r\"$\\bf{F_{MQ / PP}}}$ (-)\", fontsize=9, fontweight=\"bold\")\n", + "ax3.set_xticks(np.arange(len(categories)), categories, fontsize=9, fontweight=\"bold\")\n", + "ax3.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.35), fancybox=True, shadow=True, ncol=1, fontsize=8)\n", + "\n", + "fig.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last updated: 2024-10-11T17:56:34.996952+02:00\n", + "\n" + ] + } + ], + "source": [ + "%load_ext watermark\n", + "%watermark -idu" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "nutpie_env", + "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.10.12" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/notebooks/Normal model_normal data_noise level 0.6.xlsx b/docs/source/notebooks/Normal model_normal data_noise level 0.6.xlsx new file mode 100644 index 0000000..47b00fa Binary files /dev/null and b/docs/source/notebooks/Normal model_normal data_noise level 0.6.xlsx differ diff --git a/docs/source/notebooks/Normal model_normal data_noise level 1.2.xlsx b/docs/source/notebooks/Normal model_normal data_noise level 1.2.xlsx new file mode 100644 index 0000000..9b8ef63 Binary files /dev/null and b/docs/source/notebooks/Normal model_normal data_noise level 1.2.xlsx differ diff --git a/docs/source/notebooks/Normal model_skew normal data_noise level 0.6.xlsx b/docs/source/notebooks/Normal model_skew normal data_noise level 0.6.xlsx new file mode 100644 index 0000000..f25f740 Binary files /dev/null and b/docs/source/notebooks/Normal model_skew normal data_noise level 0.6.xlsx differ diff --git a/docs/source/notebooks/Processing_test_1_raw_data.ipynb b/docs/source/notebooks/Processing_test_1_raw_data.ipynb new file mode 100644 index 0000000..5518ff1 --- /dev/null +++ b/docs/source/notebooks/Processing_test_1_raw_data.ipynb @@ -0,0 +1,359 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "import arviz as az\n", + "import json\n", + "import numpy as np\n", + "import pandas\n", + "import pymc as pm\n", + "from matplotlib import pyplot as plt\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "raw_data_files = [\n", + " \"Normal model_normal data_noise level 0.6.xlsx\",\n", + " \"Normal model_normal data_noise level 1.2.xlsx\",\n", + " \"Normal model_skew normal data_noise level 0.6.xlsx\",\n", + " \"Skew normal model_skew normal data_noise level 0.6.xlsx\",\n", + " \"Skew normal model_skew normal data_noise level 1.2.xlsx\",\n", + " \"Skew normal model_normal data_noise level 0.6.xlsx\",\n", + "]\n", + "\n", + "parameters = [\"mean\", \"std\", \"area\", \"height\", \"alpha\", \"baseline_intercept\", \"baseline_slope\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Prepare data in df_results" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'mean'}\n", + "{'std'}\n", + "{'area'}\n", + "{'height'}\n", + "alpha skipalpha\n", + "{'baseline_intercept'}\n", + "{'baseline_slope'}\n", + "{'mean'}\n", + "{'std'}\n", + "{'area'}\n", + "{'height'}\n", + "alpha skipalpha\n", + "{'baseline_intercept'}\n", + "{'baseline_slope'}\n", + "{'mean'}\n", + "{'std'}\n", + "{'area'}\n", + "{'height'}\n", + "alpha skipalpha\n", + "{'baseline_intercept'}\n", + "{'baseline_slope'}\n", + "{'mean'}\n", + "{'std'}\n", + "{'area'}\n", + "{'height'}\n", + "{'alpha'}\n", + "{'baseline_intercept'}\n", + "{'baseline_slope'}\n", + "{'mean'}\n", + "{'std'}\n", + "{'area'}\n", + "{'height'}\n", + "{'alpha'}\n", + "{'baseline_intercept'}\n", + "{'baseline_slope'}\n", + "{'mean'}\n", + "{'std'}\n", + "{'area'}\n", + "{'height'}\n", + "{'alpha'}\n", + "{'baseline_intercept'}\n", + "{'baseline_slope'}\n" + ] + } + ], + "source": [ + "df_results = pandas.DataFrame()\n", + "\n", + "for path in raw_data_files:\n", + " for param in parameters:\n", + " # print(path, param)\n", + " # normal distribution does not have the alpha parameter so skip that when necessary\n", + " if path in [raw_data_files[0], raw_data_files[1], raw_data_files[2]] and param == \"alpha\":\n", + " print(\"alpha skip\" + param)\n", + " continue\n", + " # summary laden\n", + " summary = pandas.read_excel(path, index_col=0)\n", + " # sort summary and calculate differences between true and simulated values\n", + " df = summary.loc[param, [\"mean\", \"sd\", \"true_values\"]]\n", + " print(set(df.index))\n", + " df[\"ratio_mean_to_truth\"] = np.abs(df.loc[:, \"mean\"] / df.loc[:, \"true_values\"])\n", + " df[\"absolute_difference\"] = df.loc[:, \"mean\"] - df.loc[:, \"true_values\"]\n", + " df[\"ratio_std_to_mean\"] = df.loc[:, \"sd\"] / df.loc[:, \"mean\"]\n", + " df[\"within_range_of_1_std\"] = [True if df.iloc[x, 0] - df.iloc[x, 1] <= df.iloc[x, 2] <= df.iloc[x, 0] + df.iloc[x, 1] else False for x in range(len(df))]\n", + " df[\"within_range_of_3_stds\"] = [True if df.iloc[x, 0] - 3 * df.iloc[x, 1] <= df.iloc[x, 2] <= df.iloc[x, 0] + 3 * df.iloc[x, 1] else False for x in range(len(df))]\n", + " df[\"noise_level\"] = len(df) * [list(set(summary.loc[:,\"noise_scale\"]))[0]]\n", + " df[\"draws\"] = len(df) * [list(set(summary.loc[:,\"draws\"]))[0]]\n", + " df[\"tuning\"] = len(df) * [list(set(summary.loc[:,\"tuning_samples\"]))[0]]\n", + " # calculate mean and std of differences\n", + " df2 = pandas.DataFrame()\n", + " df2[\"path\"] = [path]\n", + " df2[\"parameter\"] = [\"\".join(set(df.index))]\n", + " df2[\"ratio_mean_to_truth\"] = [(np.mean(df.loc[:, \"ratio_mean_to_truth\"]), np.std(df.loc[:, \"ratio_mean_to_truth\"]))]\n", + " df2[\"absolute_difference\"] = [(np.mean(df.loc[:, \"absolute_difference\"]), np.std(df.loc[:, \"absolute_difference\"]))]\n", + " df2[\"within_range_of_3_stds\"] = np.count_nonzero(df.loc[:, \"within_range_of_3_stds\"]) / len(df)\n", + " df2[\"within_range_of_1_std\"] = np.count_nonzero(df.loc[:, \"within_range_of_1_std\"]) / len(df)\n", + " df2[\"noise_level\"] = list(set(df[\"noise_level\"]))[0]\n", + " df2[\"tuning samples\"] = list(set(df[\"tuning\"]))[0]\n", + " df2[\"draws\"] = list(set(df[\"draws\"]))[0] \n", + " if path in [raw_data_files[0], raw_data_files[1]]:\n", + " df2[\"data_distribution\"] = [\"normal\"]\n", + " df2[\"model_distribution\"] = [\"normal\"]\n", + " elif path == raw_data_files[2]:\n", + " df2[\"data_distribution\"] = [\"skew normal\"]\n", + " df2[\"model_distribution\"] = [\"normal\"]\n", + " elif path in [raw_data_files[3], raw_data_files[4]]:\n", + " df2[\"data_distribution\"] = [\"skew normal\"]\n", + " df2[\"model_distribution\"] = [\"skew normal\"]\n", + " elif path == raw_data_files[5]:\n", + " df2[\"data_distribution\"] = [\"normal\"]\n", + " df2[\"model_distribution\"] = [\"skew normal\"] \n", + " # save results in one DataFrame for subsequent plotting\n", + " df_results = pandas.concat([df_results, df2])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "model: normal, data: normal, noise level: 0.6\n", + "model: normal, data: normal, noise level: 1.2\n", + "model: normal, data: skew normal, noise level: 0.6\n", + "model: skew normal, data: normal, noise level: 0.6\n", + "model: skew normal, data: skew normal, noise level: 0.6\n", + "model: skew normal, data: skew normal, noise level: 1.2\n" + ] + } + ], + "source": [ + "for model in set(df_results.loc[:, \"model_distribution\"]):\n", + " dfdf = df_results[df_results.loc[:, \"model_distribution\"] == model]\n", + " for data in set(dfdf.loc[:, \"data_distribution\"]):\n", + " dfdf2 = dfdf[dfdf.loc[:, \"data_distribution\"] == data]\n", + " for noise_level in set(dfdf2.loc[:, \"noise_level\"]):\n", + " dfdf3 = dfdf2[dfdf2.loc[:, \"noise_level\"] == noise_level]\n", + " model = list(dfdf3.loc[:,\"model_distribution\"])[0]\n", + " data = list(dfdf3.loc[:,\"data_distribution\"])[0]\n", + " noise = list(dfdf3.loc[:,\"noise_level\"])[0]\n", + " print(f\"model: {model}, data: {data}, noise level: {noise}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "dfdf = df_results[df_results.loc[:, \"model_distribution\"] == \"skew normal\"]\n", + "dfdf2 = dfdf[dfdf.loc[:, \"data_distribution\"] == \"skew normal\"]\n", + "dfdf3 = dfdf2[dfdf2.loc[:, \"noise_level\"] == 0.6]" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'normal data, normal model': [[0.9999001336005343,\n", + " 1.0025702609677298,\n", + " 1.0017147282600856,\n", + " 1.000123572878139],\n", + " [0.002286555631402294,\n", + " 0.028900726978078068,\n", + " 0.02958680525019264,\n", + " 0.022445046960539197]],\n", + " 'normal data (higher noise), normal model': [[0.9997316666666668,\n", + " 1.0059567381829964,\n", + " 1.001356598861276,\n", + " 0.9977187067316658],\n", + " [0.004410296979166418,\n", + " 0.05488690135089093,\n", + " 0.055093378982298734,\n", + " 0.04168657187789078]],\n", + " 'skew normal data, normal model': [[0.9990176666666667,\n", + " 0.7598253910963016,\n", + " 0.9869124703934096,\n", + " 0.9889579711666672],\n", + " [0.04540922653553522,\n", + " 0.1425229338854569,\n", + " 0.029251994462966387,\n", + " 0.02178598822049324]],\n", + " 'normal data, skew normal model': [[0.9993873333333333,\n", + " 1.145324094260921,\n", + " 1.0038603930164334,\n", + " 1.0021702322498285],\n", + " [0.025492314214288193,\n", + " 0.06460165579288266,\n", + " 0.0295645094605588,\n", + " 0.022277250178015084]],\n", + " 'skew normal data, skew normal model': [[1.0003276666666665,\n", + " 1.0178059537564914,\n", + " 0.9995769654521169,\n", + " 0.9994046368514812],\n", + " [0.022164664598810824,\n", + " 0.08144664654979102,\n", + " 0.02553221429137138,\n", + " 0.019596288333603468]],\n", + " 'skew normal data (higher noise), skew normal model': [[0.9975454545454545,\n", + " 1.062975971807339,\n", + " 1.0078594345558298,\n", + " 1.0013061414928683],\n", + " [0.029588612507556917,\n", + " 0.13828870506270582,\n", + " 0.050852728197426554,\n", + " 0.03782158437972263]]}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "all_data = {}\n", + "for model in set(df_results.loc[:, \"model_distribution\"]):\n", + " dfdf = df_results[df_results.loc[:, \"model_distribution\"] == model]\n", + " for data in set(dfdf.loc[:, \"data_distribution\"]):\n", + " dfdf2 = dfdf[dfdf.loc[:, \"data_distribution\"] == data]\n", + " for noise_level in set(dfdf2.loc[:, \"noise_level\"]):\n", + " dfdf3 = dfdf2[dfdf2.loc[:, \"noise_level\"] == noise_level]\n", + " model = list(dfdf3.loc[:,\"model_distribution\"])[0]\n", + " data = list(dfdf3.loc[:,\"data_distribution\"])[0]\n", + " noise = list(dfdf3.loc[:,\"noise_level\"])[0]\n", + " # print(f\"model: {model}, data: {data}, noise level: {noise}\")\n", + " # print(noise)\n", + " dfdf4 = dfdf3[~dfdf3.loc[:, \"parameter\"].isin([\"alpha\", \"baseline_intercept\", \"baseline_slope\"])]\n", + " if noise == 1.2:\n", + " all_data[f\"{data} data (higher noise), {model} model\"] = [[x[0] for x in list(dfdf4.loc[:,\"ratio_mean_to_truth\"])], [x[1] for x in list(dfdf4.loc[:,\"ratio_mean_to_truth\"])]]\n", + " else:\n", + " all_data[f\"{data} data, {model} model\"] = [[x[0] for x in list(dfdf4.loc[:,\"ratio_mean_to_truth\"])], [x[1] for x in list(dfdf4.loc[:,\"ratio_mean_to_truth\"])]]\n", + "all_data" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "dict_keys(['normal data, normal model', 'normal data (higher noise), normal model', 'skew normal data, normal model', 'skew normal data, skew normal model', 'skew normal data (higher noise), skew normal model', 'normal data, skew normal model'])" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rearrange = ['normal data, normal model', 'normal data (higher noise), normal model', 'skew normal data, normal model', 'skew normal data, skew normal model', 'skew normal data (higher noise), skew normal model','normal data, skew normal model']\n", + "reordered_dict = {k: all_data[k] for k in rearrange}\n", + "reordered_dict.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# save processed data in file\n", + "\n", + "# with open('all_data.txt', 'w') as file:\n", + "# file.write(json.dumps(reordered_dict)) # use `json.loads` to do the reverse" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Last updated: 2024-10-11T18:34:57.629742+02:00\n", + "\n" + ] + } + ], + "source": [ + "%load_ext watermark\n", + "%watermark -idu" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "nutpie_env", + "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.10.12" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/source/notebooks/Skew normal model_normal data_noise level 0.6.xlsx b/docs/source/notebooks/Skew normal model_normal data_noise level 0.6.xlsx new file mode 100644 index 0000000..6686fe5 Binary files /dev/null and b/docs/source/notebooks/Skew normal model_normal data_noise level 0.6.xlsx differ diff --git a/docs/source/notebooks/Skew normal model_skew normal data_noise level 0.6.xlsx b/docs/source/notebooks/Skew normal model_skew normal data_noise level 0.6.xlsx new file mode 100644 index 0000000..451f204 Binary files /dev/null and b/docs/source/notebooks/Skew normal model_skew normal data_noise level 0.6.xlsx differ diff --git a/docs/source/notebooks/Skew normal model_skew normal data_noise level 1.2.xlsx b/docs/source/notebooks/Skew normal model_skew normal data_noise level 1.2.xlsx new file mode 100644 index 0000000..1e20add Binary files /dev/null and b/docs/source/notebooks/Skew normal model_skew normal data_noise level 1.2.xlsx differ diff --git a/docs/source/notebooks/test1_all_data.txt b/docs/source/notebooks/test1_all_data.txt new file mode 100644 index 0000000..02e8674 --- /dev/null +++ b/docs/source/notebooks/test1_all_data.txt @@ -0,0 +1 @@ +{"normal data, normal model": [[0.9999001336005343, 1.0025702609677298, 1.0017147282600856, 1.000123572878139], [0.002286555631402294, 0.028900726978078068, 0.02958680525019264, 0.022445046960539197]], "normal data (higher noise), normal model": [[0.9997316666666668, 1.0059567381829964, 1.001356598861276, 0.9977187067316658], [0.004410296979166418, 0.05488690135089093, 0.055093378982298734, 0.04168657187789078]], "skew normal data, normal model": [[0.9990176666666667, 0.7598253910963016, 0.9869124703934096, 0.9889579711666672], [0.04540922653553522, 0.1425229338854569, 0.029251994462966387, 0.02178598822049324]], "skew normal data, skew normal model": [[1.0003276666666665, 1.0178059537564914, 0.9995769654521169, 0.9994046368514812], [0.022164664598810824, 0.08144664654979102, 0.02553221429137138, 0.019596288333603468]], "skew normal data (higher noise), skew normal model": [[0.9975454545454545, 1.062975971807339, 1.0078594345558298, 1.0013061414928683], [0.029588612507556917, 0.13828870506270582, 0.050852728197426554, 0.03782158437972263]], "normal data, skew normal model": [[0.9993873333333333, 1.145324094260921, 1.0038603930164334, 1.0021702322498285], [0.025492314214288193, 0.06460165579288266, 0.0295645094605588, 0.022277250178015084]]} diff --git a/docs/source/notebooks/test2_summary.xlsx b/docs/source/notebooks/test2_summary.xlsx new file mode 100644 index 0000000..5a05df3 Binary files /dev/null and b/docs/source/notebooks/test2_summary.xlsx differ diff --git a/docs/source/notebooks/test3_df_comparison.xlsx b/docs/source/notebooks/test3_df_comparison.xlsx new file mode 100644 index 0000000..7862a78 Binary files /dev/null and b/docs/source/notebooks/test3_df_comparison.xlsx differ diff --git a/paper/Fig1_model_single_peak.png b/paper/Fig1_model_single_peak.png new file mode 100644 index 0000000..82b9a1b Binary files /dev/null and b/paper/Fig1_model_single_peak.png differ diff --git a/paper/Fig2_model_double_peak.png b/paper/Fig2_model_double_peak.png new file mode 100644 index 0000000..de6507e Binary files /dev/null and b/paper/Fig2_model_double_peak.png differ diff --git a/paper/Fig3_PP-standalone.png b/paper/Fig3_PP-standalone.png new file mode 100644 index 0000000..349ce13 Binary files /dev/null and b/paper/Fig3_PP-standalone.png differ diff --git a/paper/Fig4_peak_results.png b/paper/Fig4_peak_results.png new file mode 100644 index 0000000..2f82bff Binary files /dev/null and b/paper/Fig4_peak_results.png differ diff --git a/paper/Fig5_ppc.png b/paper/Fig5_ppc.png new file mode 100644 index 0000000..3669790 Binary files /dev/null and b/paper/Fig5_ppc.png differ diff --git a/paper/Fig6_PP-validation.png b/paper/Fig6_PP-validation.png new file mode 100644 index 0000000..a367f08 Binary files /dev/null and b/paper/Fig6_PP-validation.png differ diff --git a/paper/literature.bib b/paper/literature.bib new file mode 100644 index 0000000..e45eadd --- /dev/null +++ b/paper/literature.bib @@ -0,0 +1,213 @@ +@softmisc{nutpie, + author = {Seyboldt, Adrian and {PyMC Developers}}, + keywords = {Software}, + license = {MIT}, + title = {{nutpie}}, + url = {https://github.com/pymc-devs/nutpie} +} + +@article{scipy, + author = {Virtanen, Pauli and Gommers, Ralf and Oliphant, Travis E. and + Haberland, Matt and Reddy, Tyler and Cournapeau, David and + Burovski, Evgeni and Peterson, Pearu and Weckesser, Warren and + Bright, Jonathan and {van der Walt}, St{\'e}fan J. and + Brett, Matthew and Wilson, Joshua and Millman, K. Jarrod and + Mayorov, Nikolay and Nelson, Andrew R. J. and Jones, Eric and + Kern, Robert and Larson, Eric and Carey, C J and + Polat, {\.I}lhan and Feng, Yu and Moore, Eric W. and + {VanderPlas}, Jake and Laxalde, Denis and Perktold, Josef and + Cimrman, Robert and Henriksen, Ian and Quintero, E. A. and + Harris, Charles R. and Archibald, Anne M. and + Ribeiro, Ant{\^o}nio H. and Pedregosa, Fabian and + {van Mulbregt}, Paul and {SciPy 1.0 Contributors}}, + title = {{{SciPy} 1.0: {F}undamental Algorithms for Scientific Computing in {P}ython}}, + journal = {Nature Methods}, + year = {2020}, + volume = {17}, + pages = {261--272}, + adsurl = {https://rdcu.be/b08Wh}, + doi = {10.1038/s41592-019-0686-2} +} + +@article{matplotlib, + author = {Hunter, J. D.}, + title = {Matplotlib: A 2D graphics environment}, + journal = {Computing in Science \& Engineering}, + volume = {9}, + number = {3}, + pages = {90--95}, + abstract = {Matplotlib is a 2D graphics package used for Python for + application development, interactive scripting, and publication-quality + image generation across user interfaces and operating systems.}, + publisher = {IEEE COMPUTER SOC}, + doi = {10.1109/MCSE.2007.55}, + year = 2007 +} + +@softmisc{matplotlibzenodo, + author = {{The Matplotlib Development Team}}, + title = {Matplotlib: Visualization with Python}, + keywords = {software}, + month = may, + year = 2024, + publisher = {Zenodo}, + version = {v3.9.0}, + doi = {10.5281/zenodo.11201097}, + url = {https://doi.org/10.5281/zenodo.11201097} +} + +@article{RN173, + author = {Hoffmann, Matthew D. and Gelman, Andrew}, + title = {The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo}, + journal = {Journal of Machine Learning Research}, + volume = {15}, + year = {2014}, + type = {Journal Article} +} + +@article{RN150, + author = {Abril-Pla, O. and Andreani, V. and Carroll, C. and Dong, L. and Fonnesbeck, C. J. and Kochurov, M. and Kumar, R. and Lao, J. and Luhmann, C. C. and Martin, O. A. and Osthege, M. and Vieira, R. and Wiecki, T. and Zinkov, R.}, + title = {{PyMC}: a modern, and comprehensive probabilistic programming framework in Python}, + journal = {PeerJ Comput Sci}, + volume = {9}, + pages = {e1516}, + issn = {2376-5992 (Electronic) + 2376-5992 (Linking)}, + doi = {10.7717/peerj-cs.1516}, + url = {https://www.ncbi.nlm.nih.gov/pubmed/37705656}, + year = {2023}, + type = {Journal Article} +} + +@book{RN162, + author = {Kruschke, John K.}, + title = {Doing Bayesian Data Analysis}, + edition = {1st Edition}, + isbn = {9780123814852}, + year = {2010}, + type = {Book} +} + +@article{RN144, + author = {Azzalini, A.}, + title = {A class of distributions which includes the normal ones}, + journal = {Scand. J. Statist.}, + volume = {12}, + pages = {171-178}, + year = {1985}, + type = {Journal Article} +} + + +@article{RN152, + author = {Gelman, Andrew and Rubin, Donald B.}, + title = {Inference from Iterative Simulation Using Multiple Sequences}, + journal = {Statistical Science}, + volume = {7}, + number = {4}, + year = {1992}, + type = {Journal Article} +} + +@article{RN153, + author = {Grushka, E.}, + title = {Characterization of exponentially modified Gaussian peaks in chromatography}, + journal = {Anal Chem}, + volume = {44}, + number = {11}, + pages = {1733-8}, + issn = {0003-2700 (Print) + 0003-2700 (Linking)}, + doi = {10.1021/ac60319a011}, + url = {https://www.ncbi.nlm.nih.gov/pubmed/22324584}, + year = {1972}, + type = {Journal Article} +} + +@article{RN149, + author = {Hemmerich, J. and Noack, S. and Wiechert, W. and Oldiges, M.}, + title = {Microbioreactor Systems for Accelerated Bioprocess Development}, + journal = {Biotechnol J}, + volume = {13}, + number = {4}, + pages = {e1700141}, + issn = {1860-7314 (Electronic) + 1860-6768 (Linking)}, + doi = {10.1002/biot.201700141}, + url = {https://www.ncbi.nlm.nih.gov/pubmed/29283217}, + year = {2018}, + type = {Journal Article} +} + +@article{RN148, + author = {Kostov, Y. and Harms, P. and Randers-Eichhorn, L. and Rao, G.}, + title = {Low-cost microbioreactor for high-throughput bioprocessing}, + journal = {Biotechnol Bioeng}, + volume = {72}, + number = {3}, + pages = {346-52}, + issn = {0006-3592 (Print) + 0006-3592 (Linking)}, + doi = {10.1002/1097-0290(20010205)72:3<346::aid-bit12>3.0.co;2-x}, + url = {https://www.ncbi.nlm.nih.gov/pubmed/11135205}, + year = {2001}, + type = {Journal Article} +} + +@article{RN145, + author = {Vehtari, Aki and Gelman, Andrew and Gabry, Jonah}, + title = {Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC}, + journal = {Statistics and Computing}, + volume = {27}, + number = {5}, + pages = {1413-1432}, + issn = {0960-3174 + 1573-1375}, + doi = {10.1007/s11222-016-9696-4}, + year = {2016}, + type = {Journal Article} +} + +@article{RN146, + author = {Watanabe, Sumio}, + title = {Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory}, + journal = {Journal of machine learning research}, + volume = {11}, + pages = {3571-3594}, + year = {2010}, + type = {Journal Article} +} + +@article{RN147, + author = {Kumar, Ravin and Carroll, Colin and Hartikainen, Ari and Martin, Osvaldo}, + title = {ArviZ a unified library for exploratory analysis of Bayesian models in Python}, + journal = {Journal of Open Source Software}, + volume = {4}, + number = {33}, + issn = {2475-9066}, + doi = {10.21105/joss.01143}, + year = {2019}, + type = {Journal Article} +} + +@article{harris2020array, + title = {Array programming with {NumPy}}, + author = {Harris, C. R. and Millman, K. J. and + {van der Walt}, S. J. and Gommers, R. and Virtanen, P. and + Cournapeau, D. and Wieser, E. and Taylor, J. and + Berg, S. and Smith, N. J. and Kern, R. and Picus, M. + and Hoyer, S. and {van Kerkwijk}, M. H. and + Brett, M. and Haldane, M. and del R{\'{i}}o, J. F. and Wiebe, M. and Peterson, P. and + G{\'{e}}rard-Marchant, P. and Sheppard, K. and Reddy, T. and + Weckesser, W. and Abbasi, H. and Gohlke, C. and + Oliphant, T. E.}, + year = {2020}, + month = sep, + journal = {Nature}, + volume = {585}, + number = {7825}, + pages = {357--362}, + doi = {10.1038/s41586-020-2649-2}, + publisher = {Springer Science and Business Media {LLC}}, + url = {https://doi.org/10.1038/s41586-020-2649-2} +} diff --git a/paper/paper.md b/paper/paper.md new file mode 100644 index 0000000..88a8a69 --- /dev/null +++ b/paper/paper.md @@ -0,0 +1,149 @@ +--- +title: 'PeakPerformance - A tool for Bayesian inference-based fitting of LC-MS/MS peaks' +tags: + - Peak fitting + - Bayesian inference + - Chromatography + - LC-MS/MS + - Python + - Uncertainty quantification + +authors: + - name: Jochen Nießer + orcid: 0000-0001-5397-0682 + affiliation: 1, 2 + - name: Michael Osthege + orcid: 0000-0002-2734-7624 + affiliation: 1 + - name: Eric von Lieres + orcid: 0000-0002-0309-8408 + affiliation: "1, 3" + - name: Wolfgang Wiechert + orcid: 0000-0001-8501-0694 + affiliation: "1, 3" + - name: Stephan Noack + orcid: 0000-0001-9784-3626 + affiliation: 1 + +affiliations: + - name: Institute for Bio- and Geosciences (IBG-1), Forschungszentrum Jülich GmbH, Jülich, Germany + index: 1 + - name: RWTH Aachen University, Aachen, Germany + index: 2 + - name: Computational Systems Biotechnology, RWTH Aachen University, Aachen, Germany + index: 3 + +date: 12th of August 2024 +bibliography: + - literature.bib +--- + +# Summary + +A major bottleneck of chromatography-based analytics has been the elusive fully automated identification and integration of peak data without the need of extensive human supervision. +The presented Python package $\texttt{PeakPerformance}$ applies Bayesian inference to chromatographic peak fitting, and provides an automated approach featuring model selection and uncertainty quantification. +Currently, its application is focused on data from targeted liquid chromatography tandem mass spectrometry (LC-MS/MS), but its design allows for an expansion to other chromatographic techniques. +$\texttt{PeakPerformance}$ is implemented in Python and the source code is available on [GitHub](https://github.com/JuBiotech/peak-performance). +It is unit-tested on Linux and Windows and accompanied by documentation as well as example notebooks. + +# Statement of need + +In biotechnological research and industrial applications, chromatographic techniques are ubiquitously used to analyze samples from fermentations, e.g. to determine the concentration of substrates and products. +Over the course of a regular lab-scale bioreactor fermentation, hundreds of samples and subsequently thousands of chromatographic peaks may accrue. +This is exacerbated by the spread of microbioreactors causing a further increase in the amount of samples per time [@RN149; @RN148]. +While the recognition and integration of peaks by vendor software is – in theory – automated, it typically requires visual inspection and occasional manual re-integration by the user due to a large number of false positives, false negatives or incorrectly determined baselines, ultimately downgrading it to a semi-automated process. +Since this is a time-consuming, not to mention tedious, procedure and introduces the problem of comparability between purely manual and algorithm-based integration as well as user-specific differences, we instead propose a peak fitting solution based on Bayesian inference. +The advantage of this approach is the complete integration of all relevant parameters – i.e. baseline, peak area and height, mean, signal-to-noise ratio etc. – into one single model through which all parameters are estimated simultaneously. +Furthermore, Bayesian inference comes with uncertainty quantification for all peak model parameters, and thus does not merely yield a point estimate as would commonly be the case. +It also grants access to novel metrics for avoiding false positives and negatives by rejecting signals where a) a convergence criterion of the peak fitting procedure was not fulfilled or b) the uncertainty of the estimated parameters exceeded a user-defined threshold. + +# Materials and Methods +## Implementation +$\texttt{PeakPerformance}$ is an open source Python package compatible with Windows and Linux/Unix platforms. +At the time of manuscript submission, it features three modules: `pipeline`, `models`, and `plotting`. +Due to its modular design, $\texttt{PeakPerformance}$ can easily be expanded by adding e.g. additional models for deviating peak shapes or different plots. +Currently, the featured peak models describe peaks in the shape of normal or skew normal distributions, as well as double peaks of normal or skewed normal shape. +The normal distribution is regarded as the ideal peak shape and common phenomena like tailing and fronting can be expressed by the skew normal distribution [@RN144].\\ +Bayesian inference is conducted utilizing the PyMC package [@RN150] with the external sampler $\texttt{nutpie}$ for improved performance [@nutpie]. +Both model selection and analysis of inference data objects are realized with the ArviZ package [@RN147]. +Since the inference data is stored alongside graphs and report sheets, users may employ the ArviZ package or others for further analysis of the results if necessary. + + +# Results and Discussion + +## PeakPerformance workflow +$\texttt{PeakPerformance}$ accommodates the use of a pre-manufactured data pipeline for standard applications (Fig. 1) as well as the creation of custom data pipelines using only its core functions. +The provided data analysis pipeline was designed in a user-friendly way and is covered by multiple example notebooks. + +![](./Fig3_PP-standalone.png) +__Figure 1:__ Overview of the pre-manufactured data analysis pipeline featured in $\texttt{PeakPerformance}$. + +Before using $\texttt{PeakPerformance}$, the user has to supply raw data files containing a NumPy array with time in the first and intensity in the second dimension for each peak according as described in detail in the documentation. +Using the $\texttt{prepare\_model\_selection()}$ method, an Excel template file ("Template.xlsx") for inputting user information is prepared and stored in the raw data directory. + +Since targeted LC-MS/MS analyses essentially cycle through a list of mass traces for every sample, a model type has to be assigned to each mass trace. +If this is not done by the user, an optional automated model selection step will be performed, where one exemplary peak per mass trace is analyzed with all models to identify the most appropriate one. +The automated model selection can be started using the $\texttt{model\_selection()}$ function from the pipeline module and will be performed successively for each mass trace. +The results for each model are ranked with the $\texttt{compare()}$ function of the ArviZ package based on Pareto-smoothed importance sampling leave-one-out cross-validation (LOO-PIT) [@RN146; @RN145]. + +Subsequently, the peak analysis pipeline can be started with the function $\texttt{pipeline()}$ from the $\texttt{pipeline}$ module. +Depending on whether the "pre-filtering" option was selected, an optional filtering step will be executed to reject signals where clearly no peak is present before sampling, thus saving computation time. +This filtering step combines the $\texttt{find\_peaks()}$ function from the SciPy package [@scipy] with a user-defined minimum signal-to-noise threshold and may reject a great many signals before sampling, e.g. in the case of isotopic labeling experiments where every theoretically achievable mass isotopomer needs to be investigated, yet depending on the input labeling mixture, the majority of them might not be present in actuality. +Upon passing the first filter, a Markov chain Monte Carlo (MCMC) simulation is conducted using a No-U-Turn Sampler (NUTS) [@RN173], preferably - if installed in the Python environment - the nutpie sampler [@nutpie] due to its highly increased performance compared to the default sampler of PyMC. +Before sampling from the posterior distribution, a prior predictive check is performed. +When a posterior distribution has been obtained, the main filtering step is next in line which checks the convergence of the Markov chains via the potential scale reduction factor [@RN152] or $\hat{R}$ statistic and based on the uncertainty of the determined peak parameters. +If a signal was accepted as a peak, a posterior predictive check is conducted and added to the inference data object resulting from the model simulation. + + +## Peak fitting results and diagnostic plots +The most complete report created after completing a cycle of the data pipeline is found in an Excel file called "peak_data_summary.xlsx". +Here, each analyzed time series has multiple rows (one per peak parameter) with the columns containing estimation results in the form of mean and standard deviation (sd) of the marginal posterior distribution, highest density interval (HDI), and the $\hat{R}$ statistic among other metrics. +The second Excel file created is denominated as "area_summary.xlsx" and is a more handy version of "peak_data_summary.xlsx" with a reduced degree of detail since subsequent data analyses will most likely rely on the peak area. +The most valuable result, however, are the inference data objects saved to disk for each signal for which a peak function was successfully fitted. +Conveniently, the inference data objects saved as $\texttt{*.nc}$ files contain all data and metadata related to the Bayesian parameter estimation, enabling the user to perform diagnostics or create custom visualizations not already provided by $\texttt{PeakPerformance}$. +Regarding data visualization with the matplotlib package [@matplotlib; @matplotlibzenodo], $\texttt{PeakPerformance}$'s $\texttt{plots}$ module offers the generation of two diagram types for each successfully fitted peak. +The posterior plot presents the fit of the intensity function alongside the raw data points. +The first row of Figure 2 presents two such examples where the single peak diagram shows the histidine (His) fragment with a m/z ratio of 110 Da and the double peak diagram the leucine (Leu) and isoleucine (Ile) fragments with a m/z ratio of 86 Da. + +![](./Fig4_peak_results.png) +__Figure 2:__ Results plots for a single His peak and a double Leu and Ile peak depicting the peak fit (first row) and the posterior predictive checks (second row) alongside the raw data. The numerical results are listed in table 2. + +The posterior predictive plots in the second row of Figure 4 are provided for visual posterior predictive checks, namely the comparison of observed and predicted data distribution. +Since a posterior predictive check is based on drawing samples from the likelihood function, the result represents the theoretical range of values encompassed by the model. +Accordingly, this plot enables users to judge whether the selected model can accurately explain the data. +To complete the example, Table 2 shows the results of the fit in the form of mean, standard deviation, and HDI of each parameter's marginal posterior. + +__Table 2:__ Depiction of the results for the most important peak parameters of a single peak fit with the skew normal model and a double peak fit with the double normal model. Mean, area, and height have been highlighted in bold print as they constitute the most relevant parameters for further data evaluation purposes. The results correspond to the fits exhibited in Figure 2. + +![](./summary_joint.svg){width="100%"} + +In this case, the fits were successful and convergence was reached for all parameters. +Most notably and for the first time, the measurement noise was taken into account when determining the peak area as represented by its standard deviation and as can be observed in the posterior predictive plots where the noisy data points fall within the boundary of the 95 % HDI.\\ +In the documentation, there is a study featuring simulated and experimental data to validate $\texttt{PeakPerformance}$'s results against a commercially available vendor software for peak integration showing that comparable results are indeed obtained. + + +# Conclusions +$\texttt{PeakPerformance}$ is a tool for automated LC-MS/MS peak data analysis employing Bayesian inference. +It provides built-in uncertainty quantification by Bayesian parameter estimation and thus for the first time takes the measurement noise of an LC-MS/MS device into account when integrating peaks. +Regarding peak acceptance, it improves on vendor software solutions with more sophisticated, multi-layered metrics for decision making based on convergence of the parameter estimation, as well as the uncertainties of peak parameters. +Finally, it allows the addition of new models to describe peak intensity functions with just a few minor code changes, thus lending itself to expansion to data from other chromatographic techniques. +The design of $\texttt{PeakPerformance}$ accommodates users with little programming experience by supplying convenience functions and relying on Microsoft Excel for data input and reporting. +Its code repository on GitHub features automated unit tests, and an automatically built documentation (https://peak-performance.rtfd.io). + + +### Author contributions +$\texttt{PeakPerformance}$ was conceptualized by JN and MO. +Software implementation was conducted by JN with code review by MO. +The original draft was written by JN with review and editing by MO, SN, EvL, and WW. +The work was supervised by SN and funding was acquired by SN, EvL, and WW. + +### Acknowledgements +The authors thank Tobias Latour for providing experimental LC-MS/MS data for the comparison with the vendor software MultiQuant. Funding was received from the German Federal Ministry of Education and Research (BMBF) (grant number 031B1134A) as part of the innovation lab "AutoBioTech" within the project "Modellregion, BioRevierPLUS: BioökonomieREVIER Innovationscluster Biotechnologie \& Kunststofftechnik". + +### Competing interests +No competing interest is declared. + +### Data availability +The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. + +# Bibliography diff --git a/paper/summary_joint.pdf b/paper/summary_joint.pdf new file mode 100644 index 0000000..89c93c6 Binary files /dev/null and b/paper/summary_joint.pdf differ diff --git a/paper/summary_joint.png b/paper/summary_joint.png new file mode 100644 index 0000000..a9c9497 Binary files /dev/null and b/paper/summary_joint.png differ diff --git a/paper/summary_joint.svg b/paper/summary_joint.svg new file mode 100644 index 0000000..f7bf037 --- /dev/null +++ b/paper/summary_joint.svg @@ -0,0 +1,1902 @@ + + + +meansdhdi_3%hdi_97%meansdhdi_3%hdi_97%intercept-43.947.41-57.88-30.021115.4038.691040.141185.07slope6.660.515.717.63-21.653.09-27.50-15.94noise103.637.5189.50117.26118.638.01103.52133.2911.730.0211.7011.7612.430.0112.4212.45317.1628.84263.23370.56674.3426.34623.47722.88774.9965.28653.50897.881762.6664.041639.621881.660.160.020.130.200.150.010.140.176.560.725.227.8814.931.1412.8217.14alpha2.960.392.273.71----parameter1512.3237.311879.7237.71single skew normal modelmean25.950.0125.9325.970.5318.240.021.37areaheightstdsn1581.371441.2520.7615.690.560.481950.641809.30double normal model