{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "ikUJITDDIp19", "pycharm": { "name": "#%% md\n" } }, "source": [ "# Intro to Python Acceleration and CSV Analysis\n", "\n", "## Introduction\n", "This notebook serves as an introduction to Python for a mechanical engineer looking to plot and analyze some acceleration data in a CSV file. Being a Colab, this tool can freely be used without installing anything.\n", "\n", "For more information on making the switch to Python see [enDAQ's blog, Why and How to Get Started in Python for a MATLAB User](https://blog.endaq.com/why-and-how-to-get-started-in-python-for-a-matlab-user).\n", "\n", "This is part of our webinar series on Python for Mechanical Engineers:\n", "\n", "1. **Get Started with Python**\n", " * [Watch Recording of This](https://info.endaq.com/why-mechanical-engineers-should-use-python-webinar)\n", "2. [Introduction to Numpy & Pandas for Data Analysis](https://colab.research.google.com/drive/1O-VwAdRoSlcrineAk0Jkd_fcw7mFGHa4#scrollTo=ce97q1ZcBiwj)\n", "3. [Introduction to Plotly for Plotting Data](https://colab.research.google.com/drive/1pag2pKQQW5amWgRykAH8uMAPqHA2yUfU)\n", "4. [Introduction of the enDAQ Library](https://colab.research.google.com/drive/1WAtQ8JJC_ny0fki7eUABACMA-isZzKB6)" ] }, { "cell_type": "markdown", "metadata": { "id": "ou9NjNsdJ72B", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Import Data File\n", "We will assume that the first column is time in seconds. Some example files are provided or you can load your own." ] }, { "cell_type": "markdown", "metadata": { "id": "vX2hd0kHKPz4", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Example Files\n", "Here are some example datasets you can use to do some initial testing. If you have uploaded your own data, you'll want to comment this out or not run it!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Ut2acKSEKN1C", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "filenames = ['https://info.endaq.com/hubfs/data/surgical-instrument.csv',\n", " 'https://info.endaq.com/hubfs/data/blushift.csv',\n", " 'https://info.endaq.com/hubfs/Plots/bearing_data.csv', #used in this dataset: https://blog.endaq.com/top-vibration-metrics-to-monitor-how-to-calculate-them\n", " 'https://info.endaq.com/hubfs/data/Motorcycle-Car-Crash.csv', #used in this blog: https://blog.endaq.com/shock-analysis-response-spectrum-srs-pseudo-velocity-severity\n", " 'https://info.endaq.com/hubfs/data/Calibration-Shake.csv',\n", " 'https://info.endaq.com/hubfs/data/Mining-Hammer.csv'] #largest dataset\n", "filename = filenames[4]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "id": "MaNRqOKcbNsl", "outputId": "1b96883c-ba3d-49c3-d3e8-745a1700de73", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "application/vnd.google.colaboratory.intrinsic+json": { "type": "string" }, "text/plain": [ "'https://info.endaq.com/hubfs/data/Calibration-Shake.csv'" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" } ], "source": [ "filenames[4]" ] }, { "cell_type": "markdown", "metadata": { "id": "yyguF5C1Ju8l", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Install & Import Libraries\n", "\n", "First we'll install all libraries we'll need, then import them. " ] }, { "cell_type": "markdown", "metadata": { "id": "XQIa2uQNOR_K", "pycharm": { "name": "#%% md\n" } }, "source": [ "Note that if running this locally you'll only need to install one time, then subsequent runs can just do the import. But colab and anaconda will contain all the libraries we'll need anyways so the install isn't necessary. Here is how the install would be done though:\n", "\n", "```\n", "!pip install pandas\n", "!pip install numpy\n", "!pip install matplotlib\n", "!pip install plotly\n", "!pip install scipy\n", "```\n", "\n", "You can always check which libraries you have installed by doing:\n", "\n", "```\n", "!pip freeze\n", "```\n", "\n", "We do need to upgrade plotly though to work in Colab" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "rKONh6lJYtF6", "outputId": "e0c78623-92a5-43e8-d0b7-22fddf826f09", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: plotly in /usr/local/lib/python3.7/dist-packages (5.3.1)\n", "Requirement already satisfied: six in /usr/local/lib/python3.7/dist-packages (from plotly) (1.15.0)\n", "Requirement already satisfied: tenacity>=6.2.0 in /usr/local/lib/python3.7/dist-packages (from plotly) (8.0.1)\n" ] } ], "source": [ "!pip install --upgrade plotly" ] }, { "cell_type": "markdown", "metadata": { "id": "N9lHwj3pIozY", "pycharm": { "name": "#%% md\n" } }, "source": [ "Now we'll import our libraries we'll use later." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "4s3mCpAZNAZq", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import plotly.express as xp\n", "import plotly.io as pio; pio.renderers.default = \"iframe\"\n", "from scipy import signal" ] }, { "cell_type": "markdown", "metadata": { "id": "4HshMSZdObSY", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Load the CSV, Analyze & Plot\n", "We'll load the data into pandas, display it, do some very basic analysis, and plot the time history in a few ways." ] }, { "cell_type": "markdown", "metadata": { "id": "AFbCjAMxPx0G", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Load the CSV File and Prepare\n", "Remember we are expecting the first column to be time and will set it as the index. This is loading a CSV file, but Pandas supports a lot of other file formats, see: [Pandas Input/Output](https://pandas.pydata.org/docs/reference/io.html).\n", "\n", "If you must/need to use .MAT files, scipy can read these: [scipy.io.loadmat](https://docs.scipy.org/doc/scipy/reference/generated/scipy.io.loadmat.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 447 }, "id": "sXOwidDEOrS0", "outputId": "21e88975-7eac-4088-8c7f-24599ba1012a", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\"X (2000g)\"\"Y (2000g)\"\"Z (2000g)\"
Time
0.004394-0.122072-0.122072-0.061036
0.004594-0.0610360.488289-0.366217
0.0047940.1831080.122072-0.061036
0.0049940.122072-0.122072-0.122072
0.0051940.1220720.122072-0.244144
............
27.691064-0.427253-0.152590-0.671397
27.691264-0.122072-0.335698-0.305180
27.691464-0.183108-0.152590-0.122072
27.691664-0.3051800.030518-0.244144
27.691864-0.305180-0.030518-0.366217
\n", "

138440 rows × 3 columns

\n", "
" ], "text/plain": [ " \"X (2000g)\" \"Y (2000g)\" \"Z (2000g)\"\n", "Time \n", "0.004394 -0.122072 -0.122072 -0.061036\n", "0.004594 -0.061036 0.488289 -0.366217\n", "0.004794 0.183108 0.122072 -0.061036\n", "0.004994 0.122072 -0.122072 -0.122072\n", "0.005194 0.122072 0.122072 -0.244144\n", "... ... ... ...\n", "27.691064 -0.427253 -0.152590 -0.671397\n", "27.691264 -0.122072 -0.335698 -0.305180\n", "27.691464 -0.183108 -0.152590 -0.122072\n", "27.691664 -0.305180 0.030518 -0.244144\n", "27.691864 -0.305180 -0.030518 -0.366217\n", "\n", "[138440 rows x 3 columns]" ] }, "execution_count": 91, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv(filename) #load the data\n", "df = df.set_index(df.columns[0]) #set the first column as the index\n", "df" ] }, { "cell_type": "markdown", "metadata": { "id": "p_tOgtT5P4Y5", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Basic Analysis\n", "Once in a Pandas dataframe, doing some basic analysis is SUPER easy as shown. Here's a [link to the docs](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.max.html) for the .max() function, but notice the many others readily available.\n", "\n", "The peak will be applied after first finding the absolute value to ensure we don't ignore large negative values. \n", "\n", "$$peak=\\max(\\left | a \\right |)$$\n", "\n", "Then the RMS is a simple square root of the mean of all the values squared. \n", "\n", "$$rms = \\sqrt{\\left(\\frac{1}{n}\\right)\\sum_{i=1}^{n}(a_{i})^{2}}$$\n", "\n", "The crest factor is equal to the peak divided by the RMS.\n", "$$crest = \\frac{peak}{rms}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 142 }, "id": "oxFZ6Oj0O2Pr", "outputId": "98837822-36b4-4ef3-c319-fe32e31956bc", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Peak Acceleration (g)RMS (g)Crest Factor
\"X (2000g)\"1.7090100.2491296.859939
\"Y (2000g)\"1.6479740.2793385.899566
\"Z (2000g)\"8.8502332.6875013.293109
\n", "
" ], "text/plain": [ " Peak Acceleration (g) RMS (g) Crest Factor\n", " \"X (2000g)\" 1.709010 0.249129 6.859939\n", " \"Y (2000g)\" 1.647974 0.279338 5.899566\n", " \"Z (2000g)\" 8.850233 2.687501 3.293109" ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs_max = df.abs().max() #maximum of the absolute values, the peak\n", "std = df.std() #standard deviation (equivalent to the AC coupled RMS value)\n", "crest_factor = abs_max/std #crest factor (peak / RMS)\n", "\n", "stats = pd.concat([abs_max, #combine the stats into one table\n", " std,\n", " crest_factor],\n", " axis=1)\n", "stats.columns = ['Peak Acceleration (g)','RMS (g)','Crest Factor'] #set the headers of the table\n", "stats" ] }, { "cell_type": "markdown", "metadata": { "id": "gJTBXREgRm-U", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Plot Full Time Series\n", "You can create a plot very simply once you have a dataframe by just doing:\n", "~~~\n", "df.plot()\n", "~~~\n", "Here we show how to manipulate this plot a bit with axes labels in Matplotlib which has a very similar interface to MATLAB. There are a lot of [pretty well documented examples on matplotlib's docs site](https://matplotlib.org/stable/gallery/index.html) (but their docs are confusing to navigate)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "id": "K1ZUY6PNP7ab", "outputId": "d6aaccfa-76ad-4f05-9346-50ff67638ee3", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd1xW1f/A3x+GggiIogwxB7jD3FpqOUttmOubmqOlqZW/5lfLhg3NbNk3NcsWti0rsxxpSq7MlXskKg4EVJClDIHz++Ne8AEe4GE+gOf9ej2v594zP+fec+/nnvU5opRCo9FoNJrCcLC3ABqNRqOpHGiFodFoNBqb0ApDo9FoNDahFYZGo9FobEIrDI1Go9HYhFYYGo1Go7GJKqUwRCRcRPraW47CEJGVIjLOxrA+IrJBRBJF5O2ylq0AOWaIyJf2yr+0EZFQEXmomHELvBYi0lxEdpv3bErxpQQRUSISVJI0ygtLWUVkoYi8YB73FJEzZZz3vSLye1nmkU++xS6biDQyr5lTactVVlQphZEf1m5qWb8ARaSaiFwQkZq5z5VSA5RSITYmNQG4AHgopZ4qK3mvBXLfkzLkv8B6pZS7Uup/Nsj1rIjMKkmGxX35iIi/5bMhIqNEZIeIJIlIpPlx072o8iilJiqlXi1qPFuwVlal1FdKqVvLKL/WIvK7iMSKSJyI7BSRgWWRV0XnmlAYduJmYLdSKimfc1tpCBxUeoVlaVDce1BUGgIHihD+dmBFGclSGAOBVQAi8iQwF5gF+ADXAQuAQeUpkIg4lmd+NrAcWAP4AvWAKUCCXSWyF0qpKvMDwoGngb1APPAd4AYkA5lAkvkbBaQBV8zzPWb8UOB1YBtGhVgG1Db9XIAvgRggDtgO+BQgyzvAk9bOzXweMo/vAzYBbwEXgRPAANPvc1PGNFPOvkB1jIf6rPmbC1TPRwZP4BMgEogAXgMcC8vX9G8M/AkkYjws84AvLfy/B6LM67wBaG3hVwf4xbyG24BXgU0FXKtBwG4z/DGgv+nub6YTC4QB4y3izDBl+NKUcR/QDHgWOAecBm7N756Y9+BVYLMZ/3fA2/TrCZyxUrf6WuT9A0b9SgR2ATeYfuuADCDFvGfNMF7KB82wEcDTFul6mfJm3ZdnzPt1FngAUECQ6Xc78I95nU4DMyzSOWWGzarjNwKBpjwxGK3Ur4Baucr1IzDErCtJwPAC7lNn4C+M+h9p1olqFv6Wsn4OvGZ5PYHnTDnCgXst4n0OfIChNC9h1POilvU+LOoYcBPGMxpv/t9k4ZfvvbdSZm8zr1r5+GeV7SnzPkYC91v4F1SORmbaTub5UPPaXI/xMT8N43mIAZZgvovykWM8cMgsz0Ggvek+FaPOJQJHgD4Yz1WyZXpAO/PeOBf4ji2vl3l5/MyLvc28ILXNCzgR6y+AGVi8AC0qUoR5w9yApVlhgIcxvjRqAI5AB4xuIswb+2uutA4Dza2dk1dhXDFvuCMwCeNlIbkfPPP8FWArxpdOXWAL8Go+1+Mn4EOzLPXMa/Owjfn+hfGCrY7xZZ5IToXxAODOVQW228LvW7OCu5nXMoJ8FAbGSyge6Gc+JPWBFqbfBowvXBegLXAe6G1x/1KA2wAnYDGG0psOOJvlOpHfPTHvwTGMF7qreT7b8iVgpW5ZKowrwDAzr6fNvJ1z31/zPBLoYR57YT7M5vkI4BvzuD8QzdX69zU5X8I9gWDzOrUxw95t7eVjugWZ17U6Rl3ZAMy18HfGeEm4m3mnW8a3cq86AF3N690I4/l63MK/IIWRztX6dAuGYmhuETYe6GaWzaUYZb0Ps45hPPsXgTGmrCPN8zqF3XsrZRbgKPArcDe5PhItyvaKeT0HApcBr6LcM+B+jI+irOv3fxjPeYB5zT7ErCdWZByO8Yx1MuUNwmjlNsdQUv4W+QWax+vI+QH2JrCw0HdsebzIy+uH8VCPtjifAyykaApjtsV5K4yve0eMF+QWoI0NcgQCYQWch5JTYVj61TArkW/uB888PwYMtDi/DQi3IoMPkAq4WriNxOhbLzBfjK6IdMDNwv/r3NfLwq+WGdfTvFZXMF/6pv8s8lcYHwLvWnFvgPGl7m7h9jrwucX9W2PhdyfG12bWl7o7Fl+G+dyD5y3OJwOrzGNr9SWcnApjq4WfAzmVQvb9Nc9PYXxweFgp5xfAGPP401z1rxkWL2ErcedmXTusvESthL8b+MfivA/wh3l8LxBVxOftceAni/PCFIZlfVoCvGARdnEheRVYVnIqjDHAtlzx/wLuK+ze55N3AEZr6hhGT8UGoKlF2ZJzyXIO6FqEcjyN0SoIsAh3COhjce6H8Vzlub/AauD/rLgHmbL0JVfLAXgIWGceC4Ziubmwe14VxzCiLI4vA0Ud4DxtcXwS46vBG+PBXg18KyJnRWSOiDjnk8ZAYGUB5/nKrJS6bB7mJ7e/KZeljP5WwjU0ZY80B+riMF7O9WzI1x+4qJS6lCsfwOhjFpHZInJMRBIwXqZgXKe6GF9Mua9jfjTAeBBz4w/EKqUSc6VT3+I82uI4GbiglMqwOM8qD1i/ByWpK9nlU0plYnRLWLsPYHQ1DAROisifInIjgIg4YLQAVpnh/CnguolIFxFZLyLnRSQeo/XsnZ+A5gy7b0UkwrxPX+YKP5CrYycxgHdBg+Yi0kxEfhWRKDO9WQXlnwtr9cnyelmWu8hlzUXuZyQrP8u6Y/Xem7O7kszfcwBKqTNKqUeVUoEYz9UljBZtFjFKqfR80rOlHM8A85VSlhNzGgI/WTy7hzA+oHyslNfqM6SUCsNQ6jOAc2ZdyLrmS4EbRcQPowchE9hoJe0cVEWFYQ1loxsYFz+L6zC0+gWl1BWl1MtKqVYY/aN3AGPzScPyQbR2XhLOYlQmSxnPWgl3GqOF4a2UqmX+PJRSrW3IIxLwEhG3XPlkMQpj3KEvRquikekuGN1G6eS9jvlxGuPrPzdngdoi4p4rnQgb5LdGUe7BJYwWF5A9CFs3V5gGFv4OGF+h1u4DSqntSqlBGMr6Z4yvazC6EE4qpc6b55EUfN2+xhjTaaCU8sRoPUtWNlaynmW6ByulPIDRFuEh5zX5C6O+3G2tDCYfYHTrNTXTey5XegVhrT5ZXq/c8he1rJbkfkay8iu07ihjdldN85dn5ppS6jQwH6Pb0BYKKkcWtwLPi8hQC7fTGGOKtSx+Lkopa2XI7xlCKfW1Uqo7xvVQwBum+0WMsZt7MJ7nb5XZ3CiIa0VhRAN1RMQzl1sj82G3ZLSItBKRGhj9kj8opTJEpJeIBJsvjwQMRZKZOyMzXmdgvbXzUuAbjMpVV0S8gRcxvhxzoJSKxKgQb4uIh4g4iEigiNxSWAZKqZPADuBlcypqd4wunyzcMV4uMRgv1lkWcTMwBlJniEgNEWkFjCsgu0+A+0WkjyljfRFpYT6YW4DXRcRFRNoAD1ora2EU4x78C7iIyO1mK/J5jH5kSzqIyBDzi/xxjOux1Ure1cw1Ap5KqSsYdSer3gwEfrMIvgS4z6L+vZQrOXeMVleKiHTGeNCzOG+m2yRX+CQgXkTqY3zJZsnVGGOyxCEApVQ8Rl2aLyJ3m/fOWUQGiMgci/QSgCQRaYEx7lUUsupTD4wPru8LCFvUslqyAmhmThF2EpF7MLqXfy2ivIiIl4i8LCJBZv30xuieznOvi1GOLA5gjCHNF5G7TLeFwEwRaWjKUVdE8put9jHwtIh0EIMgEWkoxnqg3iJSHWO8L2vyTxZfY3z0DjOPC+WaUBhKqcMYL9rjZhPPn6uVNUZEdlkE/wKjTzUKY/Ata+GVL8bMmASM5uGfZlhE5DkRyeru6A38pZRKyee8pLyG8TLfizEzaJfplrV4yXI651igGkb/6EVTfj8b8xkFdMGYofQSOZvgizGa+BFm2rkfnkcxmuRRGNfyM0tPETkgIvcCKKW2YQz4vYsx8PknV78OR2K0Xs5iDOC/pJRaa6P8lhTpHpgvz8kYD2IERosj9+KsZRhfZ1mDq0NMhWCNMUC42Y0zEWO8AHJNp1VKrcTo416HMQC6Llc6k4FXRCQR4+W+xCLuZWAmsNms412Bl4H2GNf1NwxFnkWeqbxKqbeBJzEU5HmML9dHMVpFYPS1j8KYALEIY5aYrURhXKuzGLO1JprPZX4UtayW5YjBUEhPYXzU/Be4Qyl1oQjyZpGGUQfXYjz7+zE+Du6zMX6+5cgl8x5T5kUiMgB4D6Nl8rsZdyvG8wiA2WXWw4z7Pcb1+Brj3vyMMfBfHZiNMbEhCqOF+6xFtr8ATTHGrvbYUhixoRVyzSAioRgDux+XII0FwH6l1AJr59ciInIfxiBwkReAlVL+Fe4eiIgPxnTL+rZ0BZRB/iuAeUope63/0FRCKs2S9ErEbozpt/mda8qfingPPIGn7KEsTEIpvW5SzTWCVhiljFLqo4LONeVPRbwHSql/McZK7JX/nMJDaTQ50V1SGo1Go7GJa2LQW6PRaDQlp0p2SXl7e6tGjRoVK+6lS5dwc3MrPGAlRJetcqLLVjmpbGXbuXPnBaVU7vVGOaiSCqNRo0bs2LGjWHFDQ0Pp2bNn6QpUQdBlq5zoslVOKlvZRKQgiwyA7pLSaDQajY1ohaHRaDQam9AKQ6PRaDQ2USXHMDQaTcXlypUrnDlzhpSUFDw9PTl06JC9RSoTKmrZXFxcCAgIwNk5P2Pb+aMVhkajKVfOnDmDu7s7jRo1IikpCXd398IjVUISExMrXNmUUsTExHDmzBkaN25c5Pi6S0qj0ZQrKSkp1KlTBxFbLaNrSgsRoU6dOqSkFM8WqlYYGo2m3NHKwn6U5NprhaG5JlFKsSxsGSnppWV1XqOp+miFoakSJGQkMP738VxMuWhT+M1nN/P85ueZu2tuGUumqchkWYTI+t+xYwetW7cmLS0NgGPHjtGkSRMSEhLyxI2MjOSOO+4AYM2aNXTo0IHg4GA6dOjAunVXtzLZuXMnwcHBBAUFMWXKlKw9tYmNjaVfv340bdqUfv36cfGiUXeVUkyZMoWgoCDatGnDrl278uSdmxEjRnD06NF8y1VaaIWhqRKEJoSyNXIrS48uLTDcheQLHI49TFJaUvY5QHBIMBN+n1DmcmoqNh07duSWW27hrbfeAuCRRx5h5syZeHh45An7zjvvMH78eAC8vb1Zvnw5+/btIyQkhDFjxmSHmzRpEosWLeLo0aMcPXqUVauMLdxnz55Nnz59OHr0KH369GH27NkArFy5MjvsRx99xKRJhW9sOGnSJObMKXsDxFphaErMxjMb+eXYL/YWwypLjizhvV3v8XOYsWlcryW9GL58eLb/6vDV2cd/Rf7FvvP7yl1Gjf2oW7dujn+AWbNmsWjRIubMmUN6ejojR460Gnfp0qX0798fgHbt2uHv7w9A69atSU5OJjU1lcjISBISEujatSsiwtixY/n5Z6MuLlu2jHHjjN2Lx40bl8N97NixiAhdu3YlLi6OyMhIMjMzmTx5Mi1atKBfv34MHDiQH374AYAePXqwdu1a0tPT8y1XaaCn1WpKzOQ/JgMwfdN0/hnzD04OeauVUopF+xYxpOkQvF29i5T+lHVT2Byxmd+G/Iavmy8AC3YvIORACH/f+3eOsMfjjpOcnkxsSiz1a9bn1a2vZvvdHXR39nFyenL28ZaILdnHo1aM4osBX9C2XtsiyagpHm/8foyjF5ILD1gEWvl78NKdrW0Ku3379hz/ALVq1WLatGlMnjyZgwcPWo134sQJvLy8qF4991bvhiJp37491atXJzw8nICAgGy/gIAAIiIiAIiOjsbPz9gx2dfXl+joaAAiIiJo0KBBnjibN28mPDycgwcPcu7cOVq2bMkDDzwAgIODA0FBQezZs4cOHTpYLVdpoFsYmhJxPP54jvOsQeTYlFi+Pfxttvu+C/t4/5/3uefXewgOCSY8PpwNZzbk+aKfs30O438fn8Nt/en1pGWm0e+HfgBEX4rmgz0fcDn9MsEhwXx3+DvWJKwBYPnx5XT+qjP9l/bnvlX35UjHUjEoru4D8/Dah3OEG7NyDLEpsUW5DJoqxsqVK/Hx8clXYURGRlr9ej9w4ABTp07lww8/LFJ+IlLo7KVNmzYxfPhwHBwc8PX1pVevXjn869Wrx9mzZ4uUb1HRLQxNsbiQfIFeS3rlcb/xmxt54PoH+HT/pwB08OlAU6+m3LviXgDOXT4HQMjBEH7412hO7xq9i+mbp1PDqUb2GMSmiE341/Rn2oZpOdL/O/JvHvr9oRxur/39mlUZd0bvzHFuqRhe2vJSgeXbFrWNJUeWsD3K+EL74c4faF67eYFxNEVn6q2BFW5x26+//kp8fDyrV69m8ODB3HbbbdSoUSNHGFdX1zxrGc6cOcPgwYNZvHgxgYGBJCYmUr9+fc6cOZMjTP369QHw8fEhMjISPz8/IiMjqVevHgD169fn9OnTVuMUREpKCq6ursUuty3YrYUhIs1FZLfFL0FEHs8VpqeIxFuEedFe8mquciXzCuNWjsvXP0tZAAz5ZQi/Hv81T5gsZQHQ/sv2rDyxMseA9aS1kxj08yAOxeY0rZBbWZQVz/z5TLayAHhtq3WlpKlaJCcn8+STTzJ//nyCg4MZNGgQM2fOzBOuWbNmhIeHZ5/HxcVx++23M3v2bLp165bt7ufnh4eHB1u3bkUpxeLFixk0aBAAd911FyEhIQCEhITkcF+8eDFKKbZu3Yqnpyd+fn5069aNpUuXkpmZSXR0NKGhoTlk+vfff7n++utL+YrkxG4KQyl1RCnVVinVFugAXAZ+shJ0Y1Y4pdQr5Sulxhof7vmQU4mnbA7/7MZny1Ca8uFY3DF7i6ApB1599VUGDx5Mq1atAJgxYwbffPNNjimrAG5ubgQGBhIWFgbAvHnzCAsL45VXXqFt27a0bduW8+fPA7BgwQIeeughgoKCCAwMZMCAAQBMmzaNNWvW0LRpU9auXcu0aUZreuDAgTRp0oSgoCDGjx/PggULABg6dCgBAQG0atWK0aNH0759ezw9PQFjPMTV1RVfX98yvT4VpUuqD3BMKVXoBh4a+5PVrXQtkXglkfjUeDyre9pbFE0ZMmvWrBzn7u7uHD9+3GrYRx99lM8//5zXXnuN559/nueffz6Hf2JiImBM1d2/f3+e+HXq1OGPP/7I4y4izJ8/P4+7g4MDb731FjVr1iQmJobOnTsTHBwMwNdff83DDz+cJ05pU1EUxgjgm3z8bhSRPcBZ4Gml1AFrgURkAjABjL7B3M01W0lKSip23IpOaZVt/Zn1JRemEjJs6TCm+08v93yrWp309PTMfplmZGRkH1c2+vbtS0RERL7yl0XZBg4cSHx8PGlpaTzzzDO4ubmRmJiIi4sLQ4YMsTm/lJSUYtUpyVp1aC9EpBqGMmitlIrO5ecBZCqlkkRkIPCeUqppYWl27NhR6S1a81JaZQsOCS65MJWUfePKf51GVauThw4domXLlkDFtOhaWlTkslnegyxEZKdSqmNB8SrCtNoBwK7cygJAKZWglEoyj1cAziJStEn8Go1GoykVKoLCGEk+3VEi4ivm5GQR6Ywhb0w5yqbRaDQaE7uOYYiIG9APeNjCbSKAUmohMAyYJCLpQDIwQtm7D02j0WiuUeyqMJRSl4A6udwWWhzPA+aVt1wajUajyUtF6JLSaDQau2BpBjwlJYUWLVqwb9/ViQ1vvvmm1emqycnJ3HLLLWRkZLB7925uvPFGWrduTZs2bfjuu++yw504cYIuXboQFBTEPffck202PTU1lXvuuYegoCC6dOmSYxHg66+/TlBQEM2bN2f16tW5s87D008/ncOces+ePQkPDy910+agFYZGo9EA4OLiwty5c5k8eTJKKSIiIli4cGG22XFLPv30U4YMGYKjoyM1atRg8eLFHDhwgFWrVvH4448TFxcHwNSpU3niiScICwvDy8uLTz75BIBPPvkELy8vwsLCeOKJJ5g6dSoABw8e5Ntvv81Oa/LkyWRkZBQo92OPPWZVxrJAKwyNRnPNktsMeP/+/fHz82Px4sU88cQTzJgxAy8vrzzxvvrqq2xTHs2aNaNpU2O2v7+/P/Xq1eP8+fMopVi3bh3Dhg0D8powzzJtPmzYMP744w9jF8hlyxgxYgTVq1encePGBAUFsW3bNsBYhd68eXO6d+/OyJEjs/fsaNiwITExMURFRQFQu3ZtHB0dS920OVSchXsajeYapPr6lyDmSOkm6hsMA2z74rZmBnzu3Ll07tyZpk2b5tgIKYu0tDSOHz9utctn27ZtpKWlERgYyKlTp6hVqxZOTsZr1tK0uaUJcycnJzw9PYmJiSEiIoKuXbtmp5cVZ/v27SxdupQ9e/Zw5coV2rdvT4cOHbLDtW/fns2bNzN06FB+/PHHPGUqLbTC0BSJvef32lsEjaZM8ff3p3fv3tnbr+bmwoUL1KpVK497ZGQkY8aMISQkBAeH0u282bx5M4MGDcLFxQUXFxfuvPPOHP7lYdoctMLQFJEsM+XXKtqeVOmS2utlqlXA1dAODg75vvStmTZPSEjg9ttvZ+bMmdkthNq1axMXF0d6ejpOTk45zJRnmTAPCAggPT2d+Ph46tSpk69pc0sT6dYoD9PmoMcwNJoiEX05j0ECzTWGl5cXGRkZ2UojLS2NwYMHM3bs2OzxCjCMCPbq1St7G9XcJsyzTJv/8MMP9O7dGxHhrrvu4ttvvyU1NZUTJ05w9OhROnfuTLdu3Vi+fDkpKSkkJSXx6685twwoD9PmoFsYGo1GU2RuvfVWNm3aRN++fVmyZAkbNmwgJiaGzz//HIDPP/+cwMBA3njjDUaMGMHzzz9Pu3btePDBBwF48MEHGTNmDEFBQdSuXZtvvzV2p2zdujX/+c9/aNWqFU5OTsyfPx9HR0c6derEXXfdRZs2bfDx8SE4ODjbtPmVK1cICwujY8cCzUCVClphaDQaTS6yXvz58cgjj/Duu+/St29fRo8ezejRo/OESUxMpEmTJtmznCxxcXHh+++/t5r29OnTmT49r1Xkp59+mhkzZnD58mVuvvnm7EHvX3/9lWHDhmUPrpclWmFoNBpNEWnfvj29evUiIyMDR0fHcslzwoQJHDx4kJSUFMaNG0f79u0BSE9P56mnnioXGbTC0GiKwKmEUzTzamZvMTQVgAceeKBc8/v666+tug8fPrzcZNCD3hpNEXgi9Al7i6DR2A2tMDQajUZjE1phaDQajcYmtMLQaDQajU1ohaHRaK5ZLM2bA8yfP5+2bdtm/66//npEhEOHDuWJGxkZmW0+ZM2aNXTo0IHg4GA6dOiQw9z4zp07CQ4OJigoiClTppC1B1xsbCz9+vWjadOm9OvXj4sXLwKglGLKlCkEBQXRpk0bdu3aVWg5RowYwdGjR/MtV2mhFYZGo9GYPPLII+zevTv7d9ddd3HvvffSsmXLPGHfeecdxo8fD4C3tzfLly9n3759hISE5DBaOGnSJBYtWsTRo0c5evQoq1atAmD27Nn06dOHo0eP0qdPn2wT5StXrswO+9FHHzFp0qRC5Z40aRJz5swpjUtQIFphaDSaa5bc5s0t2bBhA0uWLGHBggVW4y5dupT+/fsD0K5dO/z9/QFjtXZycjKpqalERkaSkJBA165dERHGjh1r1cR5btPnY8eORUTo2rUrcXFxREZGkpmZyeTJk2nRogX9+vVj4MCB2WZHevTowdq1a0lPTy+0XCVBr8PQaDR2Y+6euRxPOl6qabao3YKpnafaFNaaeXOAuLg47rvvPr744gs8PDzyxDtx4gReXl5Ur149j9/SpUtp37491atXJzw8nICAgGw/SxPn0dHR+Pn5AeDr60t0tGGnzNL0uWWczZs3Ex4ezsGDBzl37hwtW7bMXgvi4OBAUFAQe/bsoUOHDvmWq6ToFobGZq5kXrG3CBpNuTBx4kTGjBlDt27drPpHRkZa/Xo/cOAAU6dO5cMPPyxSfiKCiBQYZtOmTQwfPhwHBwd8fX3p1atXDv/yMHFu9xaGiIQDiUAGkK6U6pjLX4D3gIHAZeA+pVTho0CaUuez/Z/ZW4QKQXpmOk4Odn90qgSP3/A47hXMvHlISAgnT57kyy+/zDeMNRPnZ86cYfDgwSxevJjAwEASExPzmCa3NHHu4+NDZGQkfn5+REZGUq9ePYB8TZwXRnmYOK8oLYxeSqm2uZWFyQCgqfmbAHxQrpJpsrmQfMHeIlQIfg//3d4iaMqI48eP89xzz/HVV18VaMyvWbNmhIeHZ5/HxcVx++23M3v27BytEj8/Pzw8PNi6dStKKRYvXmzVxHlu0+eLFy9GKcXWrVvx9PTEz8+Pbt26sXTpUjIzM4mOjiY0NDSHTOVh4ryiKIyCGAQsVgZbgVoi4mdvoa5Fvjn8jb1FqBDorrmqyxtvvMHly5cZMmRIjum1GzduzBHOzc2NwMBAwsLCAJg3bx5hYWG88sor2XHOnz8PwIIFC3jooYcICgoiMDCQAQMGADBt2jTWrFlD06ZNWbt2LdOmTQNg4MCBNGnShKCgIMaPH5896D506FACAgJo1aoVo0ePpn379tkmzqOjo3F1dcXX17dMr49kzQm2FyJyArgIKOBDpdRHufx/BWYrpTaZ538AU5VSO3KFm4DRAsHHx6dDln35opKUlETNmjWLFbeiU9KyPXbysVKUpnLzfsP3yy2vqlYnPT09CQoKAihXa6+lzfLly/nnn3948cUXrfqXRdmy6kJMTAy9evVizZo1+Pj4MG/ePDw8PBg7dqxN6YSFhREfH5/DrVevXjvz6eXJpiJ0xHZXSkWISD1gjYgcVkptKGoipqL5CKBjx46qZ8+exRImNDSU4sat6JS4bCGlJkqlpzzrSFWrk4cOHcoet0hMTKxwYxi2MmrUKC5fvpyv/GVRtjvvvJO4uDjS0tJ46aWXshWvr68vY8aMsXlPDBcXF9q1a1fk/O2uMJRSEeb/ORH5CegMWCqMCKCBxXmA6abRaDR25aGHHirX/HKPW2Rx//33l0v+dh3DEBE3EXHPOgZuBfbnCvYLMFYMugLxSqnIchZVo9GUIvbuCr+WKcm1t3cLwwf4yZx/7LaLkJEAACAASURBVAR8rZRaJSITAZRSC4EVGFNqwzCm1ZaPKtXk4Pzl8/YWQVNFcHFxISYmhjp16thblGsOpRQxMTG4uLgUK75dFYZS6jhwgxX3hRbHCnikPOXS5KX3973tLYKmihAQEMCZM2c4f/48KSkpxX55VXQqatlcXFxyrD4vCvZuYWg0lZKEtAQ8quU1GaEpHGdnZxo3bgwYffLFGXytDFTFslWGdRgaTYVjy9kt9hZBoyl3tMLQaIrBM38+Y28RNJpyRysMjUaj0diEVhgajUajsQmtMDSaYpKQlmBvETSackUrDE2h6EVW1rl85bK9RajSZKpMktOTyz3PcSvH8efpP8s138qCVhiaQtkeVbq7dlUVIi9pgwOlxeUrl/O8pGdvm03nrzqTnpleYNzDsYcJjw8vUf77zu8jIimClPQUdp3bxTMbCp7UsCt6F9GXokuUZ2VEKwxNoeh9MKwzdqVtlkE1MHHNRO76+S5mbJlh1Tz8a1tf49F1j3L04tFst6X/LgUgLC6M1/9+nZjkGF7Y/AIp6Tk3Lhq+fDh3/nwnC/cspCAyMjMIDglmxpYZBIcEE3YxLNtv1IpR9F/an9DTofnGPxx7mOXHlgMwbtU4Bi8bbDXclYwrdP6qM7/HV719U/TCPU2hHIs/Zm8RNJWczWc3A3Ai/gRLjy7NdnfEkWcPP8veC3sBQzkE1Qri+3+/Jy0zDYDH1j1G1KUovj78NQA/h/3MKze9AkB1x6t7as/fPZ/5u+cD8ObNb3JzwM08u/FZ1p1eB8D7vQ2T9Fn5D/4l7wt/6kZjL/Dk9GSWHFlCeEI4Q4KG5Ah7Z+CdACReSeRUwin2XtjLsxufZeINE1m4ZyGTb5hMcnoyy+OWM4tZJbpuFQ2774dRFnTs2FHt2LGj8IBWqGqmpC0pbtmCQ4JLX5gqwr5x+8o8j6pQJ6tSHdo3bp/N5SmP+lFaiEih+2HoLimNRlOmVLWP0qqk/IqKVhiaAqlqD3tps+bkGnuLUOHJ6krSVH60wtAUyIYzRd788Jpi34XK0+VgL2Zvm21vETSlhFYYmgL5/MDn9hZBo9FUELTC0BTIjujiTR64VohPjbe3CBpNuaEVhkZTAlIzUu0tQoXm430f21sEu7L73G57i1CqaIWh0WjKjPd2vWdvEezKq1tftbcIpUqBC/dExAW4A+gB+APJwH7gN6XUgbIXT6PRaCov/178194ilCr5tjBE5GVgM3Aj8DfwIbAESAdmi8gaEWlTLlJqNBUUbaROUxhHYo/YW4RSo6AWxjal1Ev5+L0jIvWA64qbsYg0ABYDPoACPlJKvZcrTE9gGXDCdPpRKfVKcfPUFI3TCaftLUKFJ+lKkr1F0FRwjlw8QvPaze0tRqmQr8JQSv1WUESl1DngXAnyTgeeUkrtEhF3YKeIrFFKHcwVbqNS6o4S5KMpJs9vft7eIlQKTiacpKFHQ3uLoamgVKXFr4UaHxSR5RgtAEvigR3Ah0qplLyxCkcpFQlEmseJInIIqA/kVhgaO5CpMtl1bpe9xagU3PHTHZXKZpBGU1xssVZ7HKgLfGOe3wMkAs2ARcCYkgohIo2AdhhjJbm5UUT2AGeBp/MbbBeRCcAEAB8fH0JDQ4slS1JSUrHjVnSKUrY9l/eUrTBVjLKsM5W1Th5L0VaOAV7c/CKeZzztLUapUKi1WhHZrpTqZM1NRA4opVqXSACRmsCfwEyl1I+5/DyATKVUkogMBN5TSjUtLE1trdY6RSnbtWxgrTiUZQujstbJ/kv7E5EUYW8xKgSVoQVaWtZqa4pI9uC2eVzTPE0rgXyIiDOwFPgqt7IAUEolKKWSzOMVgLOIeJckT42mLPhwz4f2FqHCoZXFVYJDgjmZcNLeYpQYWxTGU8AmEVkvIqHARuBpEXEDQoqbsYgI8AlwSCn1Tj5hfM1wiEhnU96Y4uap0ZQV83bPs7cImgrOFwe/sLcIJabQMQyl1AoRaQq0MJ2OWAx0zy1B3t0wxj/2iUjW+vnnMKfqKqUWAsOASSKSjrFocISqSlMONJoqyrW433VhfHfkO6Z3mY75DVwpyVdhiEh3pdQmAKVUKrAnl78HcJ1San9xMjbTLvDKKaXmAfrTTaOpZExcO9HeIlRIktOTqeFcw95iFJuCuqSGisgWEXlRRG4Xkc4icrOIPCAiXwC/Aq7lJKdGU+EZ+stQEtIS7C2G3UlOTyYsLszeYlRIKnPrAgpQGEqpJzDsSEUCw4FXgSeBphjrL25WSm0vFyk15crRi0ftLUKl5N+L//Lm9jftLYbdeWL9E/YWocJyIfmCvUUoEQWOYSilYjHWWiwqH3E09uZI7BGGLR9mbzEqLWeTztpbBLuz+exme4tQYdkUsYmRLUbaW4xio82ba3KglUXJ2Ba1zd4iaCows/6eZW8RSoQtK701VZy95/dSv2Z9DsRoi/UajSZ/tMKobCRGwYa3oP9scCyd23fvintLJR3NtY1SStsfs4GQAyGMaTUGB6l8HTw2SSwiN4nIKBEZm/Ura8E01jn/zUTYvoi4fatKnNa6U+u0CZAyIDgkmG2R117X1I9Hf+S+VffZW4wKz1s73qq0i/gKVRjmFNq3gO5AJ/NXoL2RaxWlFB+EHiPucoksphRIZKyx/8LJ2MslTuu7I9+VOA2Ndd7ZadV4QZVlw5kNzPhrhr3FqDTsu1DxbUtZw5YWRkegm1JqslLqMfM3pawFszcZieeITUoF4FJqOslpGdl+yWkZXE5LzxNnx4EjjFh/Mx9881O+6cYkpaKU4lTkOVIuJxqOv0zhgyV3MXntZBskMxa6l2Q+d0xyDKQm4nB8Y7HT0BTMtTYe9Mgfj9hbhErF6vDV9hahWNiiMPYDvmUtSIUiah8Obzdl1uxpxF5Ko/VLq2n5otkFlJ7KC69Mp9WLebuE6u96By9J4pZY61/uJy5cYvbrz7Pis5lc92FTkt8wra3sCmFB8gk2Rlh5gV9JgUxDWWVkZpCmMgFITbduIWXfmXg2/xvF6j3WDZ1tjdxKzyU9Wf/PIrZVq3x9qJqKR1pG2bWoqzIp6cXaSsiu2PLG8AYOishqEfkl61fWgtmVc4f5xNOD1U3X8sXWf+jhsJfXnD7hp3/O8Pdnz/CW0wfc6rCDK5lXGLp0Ak1e+oSUKxn4H89SFDm//sPOJfL8xz/h/3Eb3nT+iNtPGYu7vCSJ5f+E5wj725fvsWx3BD/uOsOPSz6DmT6w9CEAntv0HPf7Gwt/Fm3KGS+LO+dtIva721m8tR93fNyTL/8JJTgkmNBPjW1L9oevB+CffxaR5lC5V51WdLZGbi1WvJT0FB774zFOJZyyLcK612CPlY+U9a/D7m/yulthQWgYX/9tY365mPePtt5THDp91YnQ06FW/fZHxBN+4VL5CmQDtkyzmVHWQlQUoi5F8djJxxhwtBsra9cCYPnG1ayt9jYAvZasYqHzH+AA9zuu5lDMEP5N+ov6/gfYcuwOarpUp0NKKplKsfGNIfRI/gNmxNP3nQ184fwW1R2vrvL8wsOdAZcu8dzeOzlo5gXQ99iLdI/ZwqC0cJo7nOaDWh54nfqdqO/vZsXlqxvSuF2JRSnF6gNR+Hm68vyOMSSmJQLPsLRuAntdXIAY3tj7GACPOe7mfeDXf06CG3xWQ0+QK2vG/z4+7z4I5w6DyoSMNPBvazXe1sithJ4J5VJcFJ8OWVJgHsEhwXRJTuHjqHPw0wTUC7EsOfgnvRvdRN0/ZxuB2ha+UGzOqiMAjOpyXSEh8/LZgc+KHEdj8Ni6x9g1Zhe7TyUSXN8TF2dHAO5/fzleksTvr+e0yRV2MYyfwn7i6Y5PIyJkpCaxY/nDdBnwP3CrU+by2mKt9k8R8cEY7AbYZu7nXeU4HnccgJXVrq5UvVEO0fW6AC45ONDz339p6BjBG161ePTiIaabX/neDhd4ZeeDxPj5ALDp5Ho8M42uo92n4wBQZqtDAU/V82aNWw3m1PEC4EtPj+z8Rvj7klLtIMb3okUFuJxz97Jx7p/QZrGx7brTsQmkBxqy1ONivuVTSnHa4V+br4em5ASHBDOyxUimdprK4b/exmP9G7hnZlIrMxNmxFuN4xB3GgCXc//A/9qS2Pxp3vt+MI8N+hKHam55wv/t6pJ9/P2/S3ht1yym/zKSLGtOX2w9yZiu5p7jG9+B5gOgXktIiYekc0Q6OeFb93uizuddtJmekcmc1Ud4+OYm1KlZHYCYY7u46dNDUGczmx5+pQRXRwPw0YZjLDw5nCuJrXA8dz/+tVzZXH0K1SQDsFAYe79n4vYXiXZyYmyDfvj4tuXz9f9l7qW9fLDqMbrf+F/wDQYHxzKT1ZY9vf8DvAmEYvS1vC8izyilfigzqezEm6sPQ65r7VLzEJccjJ67nddtoZNLA8B4yV858y3OHhBWrRqkXx0z6N4wgM8io+mYksqFA38y2fE3PGocIdj/Oj6IOscat/ytVYZVq1aonMGNr6PfpauzpNIDP8o+nlJvGm+4elmN9++esaR51bLqpyk7vjn8Dcm7V/KzSxw08KdWRgYbT0UQEvoXb52cwEsXYhiWeImw0ZupvmosF1OioG4dNtZwZc6VOL6IngNA+3eb0M01AIfYMH6qeTtud7ycnccYPx++iIzmVOxhAJydY/gXZ4YG+MGRO1j1y71cavwZHVJSGfPnTDLHLqPJL0+iLhxhWOPGXPLO4I6UVCLWhFO/36MApFzJoOWrX6DSvDkRe46tGY8ywfc2+uz8nIYNGxBZPZ3eP/Qq/wtaxVh4cjgAzu4Hwf0ZooDJNeowMS4e9/DjNKzfkP0R8QSGziPa03hl9109hm/6buJkqrE90N6z2+n+0S1w83+h9/Qyk9WWLVr3AP2yWhUiUhdYq5S6ocykKiHF3aI18NV3qXHdpzncfNPTiXIqfvfNB1Hn8MzI5HuPmvzkXrPwCJprgoVR53imrjeJjsbHSIfkFHZatBTyY9mZs4zx8yHBMe9XZKfkFE47O9lcXz8/G819/j453ByUoq3jEF7ociszN+1lh+MHNE1LwzvVlb/cM/JJSVOW1M1oTK2MIzwbc5EH/K7eL7eEQBpVP8iB6kbLb98Jcwwqn5ZrYdiyRastCmOfUirY4twB2GPpVtEorsJ4fVYfvq5fJXvbNBpNFWdRZDSdUlJxLEOFYcunyCoRWQ1kTbe4B1hRLIkqOFpZaDSaysp4Px/cMzLZUoZ52DLo/YyIDMXYUhXgI6VU/ivTNBqNRmMXsro4ywqbOjuVUkuBpWUqiUaj0WgqNAXt6b1JKdVdRBLJskdhegFKKeWRT1SNRqPRVEEK2qK1u/nvrpTysPi5l5ayEJH+InJERMJEZJoV/+oi8p3p/7eINCqNfDUajUZTdGy1VluoW1EREUdgPjAAaAWMFJFWuYI9CFxUSgUB7wJvlDRfjUaj0RQPW0ZIWlueiIgT0KEU8u4MhCmljiul0oBvgUG5wgwCQszjH4A+UhIzrRqNRqMpNgWNYTwLPAe4ikhCljOQBnyUX7wiUB84bXF+BuiSXxilVLqIxGPYy7iQKxwiMgGYAODj40NoaGgpiKjRaDSVi7J89+WrMJRSrwOvi8jrSqlny0yCUkIp9RGmIuvYsaPq2bNn0RMJKTyIRqPRVGSK9e6zEVvWYTwrIl5AU8DFwn1DCfOOABpYnAeYbtbCnDG7wjyBmBLmq9FoNJpiYIvxwYeA/8N4oe8GugJ/Ab1LmPd2oKmINMZQDCOAUbnC/AKMM/MbBqxThdky0Wg0Gk2ZYMug9/9hmDY/qZTqBbQD4kqasVIqHXgUWA0cApYopQ6IyCsicpcZ7BOgjoiEAU8CeabeajQajaZ8sGWld4pSKkVEEJHqSqnDItK8NDJXSq0gl10qpdSLFscpwPDSyEuj0Wg0JcOWFsYZEakF/AysEZFlgPUNozUaTYXiliS937am9ChUYSilBiul4pRSM4AXMLqJ7i5rwTSaqsRXZ6PKPc89J06x84r1bWA1FY8uySn2FqFQClQYIuIoIoezzpVSfyqlfjEX2lU5toWf5rPI6ALDzI0+X07SFMzI+ER7i6ApAm1S02iVmmrV77awnvzPrFeZabVz+PW8dJmOuV4k7hmZ/F+sMYzonGsOyCeR0YyLT7BwsX2da7fLyVbdv40of2V3LZIJZKT4UC1Tse/EKR65WOKh4lKnQIWhlMoAjohI0XeGr4QsqDaZjik5H+oH4nJuRlI7w/quY5ZbppY2K09H5HoJwH9jr+7d/fnZq0pu14lT3JBi/cWkKX+6JKfwcWQ0+9y6EpzkajXMpQwvel1O5uGjN3Dp2DPsPXGKm82X9y2JznRKvpInzoPxCew5cYrhR2/K4d45JZUJMak8eLQDL1x5gJ41+9r05eqamckH0edZdiyVvVk7twEvXIildVqV/D6scPyV2ZLLJx5nx0ljPfPDcQk57kUW7kcn8sepCDonp/DnyTNMvBjPCxdiy0VGW8YwvIADIvKHiPyS9StrwezBU1Nf43CtW3K43R+fyE3mw/vq+Rh8cimM96POMz4unnfO5Vl8Xir8fiKWgPQMno69+rXxd/hpnDC25Hw7+jxemVdlcgbez9UKap6qH3h78VHUObqkpNLSz5NeLf+bx18pYWVmZ9Y2fYH/pQ8DBAFuNF/ytYd9z9tJjwPQPsVwc0Aht73Onnavce+Tc7PTmhd1jgwleLxwksdf+5yZM9/l9RF38XHUOXo6Xt3LvXFaXgX065lIBJh65WEE2HTyDG+eu8CtDx3MEc4r0Y+0CD0PpSxo5utG1v0Ho20owA8RkQxIusTf4acZFXUHZ9MbUS8jg0+izlE7M5NH4uL5T2IS805lwMl7y1RGWxTGC8AdwCvA2xa/qoeDA1Ftn2TjPRvpY7YYamVmZt9A74wM/NMzWH/qTHaUnsnJTLlotEKyvuQet/j6Ly57T5xi34lT1H5kI9RqaMiSWgOAaamTAfg06hxr4sfgkZmZHe+C8mBlWk+6O3plu02LvYh7RiZtdMujXJkbfT77AXNycafbgHtpXauZ4XBoOolHXibpyAxAaNp/EuJk7M18LNOPe29+lXXD19E7MJi3O7VhpVdfhkXVB8xOphsn027QY1xXp0Z2fusT/kOz1MVgpgNAzXrw9FHeuWcNK+5axt/hp1liZTylnvkhtF214IG0p/HMzKT/pcvU8vTMEa5fm8F8OepR1g1fx7o+n9BMf4yUiM7JKdS+5A2Ah4sT6566hYiA27P9l2XcRPO0K8w5H0MNpcjAzWo64Zk+jLv0Ju/f80CZymvLSu8/RaQh0FQptVZEagB5d6CvQtRyqcVcixZDVj/xX6odjdRRAjKstCYGvsXsVf/lF3c37o9PZG5tr7xhcnElIRhnj3053C6deASf6mEIiwCoXrcJTNoMaZc4PWsL4pjCb5k1+V9PT2h2G2/63UDXZ4NZcfpx/NIzeDZ9As8+N5PTfz3Npoj1AHRMSaXbsYH06dSILWfe4if3msW9NJoi0ODGJ6DlEDj4M3R6CIDPbv+S+NR4PJ3r8uGGY7Tw9cDD1YmGddw48toAGk37jXtd5rO1Yx/qmunUcnEgoP+71Lg1Fr67Jd/8Rj42i1anrfR716yHM9DAqwnc/SHUbQ5rxgBw+fRYejYNoHmKC14kMveetvx1rAHsf8tqHo28a9KpkTnO4upNdY8ASNVbGxeHx2Mv8mB8Ij8O/4yXdkwik0ya1K0JD30NGem88+LDfJwxkEEv/wpfDoXj6+na83Y+/ulsnrRSqMaU3kH0bF6vTGW2ZaX3eAyjfrWBQAyDgAuBPmUqWQXipZhYmly5QvfbQoiv04DxH3wPfJAzUDU36ng04IE4Y8ZxjROjudz4y3zT/PFMJP2SXsWh2gUcXSJ5qdV0Gvt1poGnP6FHzsOvi64Gru4O1d3ZMnUAN81eRwZAz6nZ3uepRYP0DDKU8H1GT6YpRY8W9/BlxHo6XTa+AL/P6MkAzxq8sjuWvpcu84jv1Yo1MOkSK2pa/3LRFI+pnabSrNVo4+TmZ7LdXZ1ccXUyxjIe79ssT7w/n+lJLddqVtPMMnKQexjb3dmdujXq0tzXnea+7gULdsM9OU5/e+hhGnu70fyfVURRh7vb1Sc4wBP2Q1KdYHJ/Wrg5W9QTEcTdVyuMYuKsgC6TCKpnZXshRyf+lzGEB7o1BgdHGPszAH2B8C7t4NIxeDPQCHvLVMT3Tv6ved76VNrYsnDvEQxT5H8DKKWOikjZqrEKhndGJu9FfkAX9+tofZ03Qwf2p97fr4FUgxtGQexxuH4osv1jMBVGekqDAtP86PL9oJzJTPPG0SWSzs26cZ2nEec/HRvAr3nj+NeyPmj63IAWsB6Uxauks38X0mK7siOmE6H9Agm/pU+2Fcubk1MYnzCWRR6LjfiJGazQjY5SZXSWsigiDevkr7idHZ0BaKpyPrZbRm0pcj4L+izgVOIpWvjmfVkF1q0JE0Kp6dXYcHjyEBP2LeKjsB/wqeGTI2wjz0bsvbCXDnV6sTNmfZHluJaR/q9DqzG0zsxgZIuRjG01Nod/+Ozb84kJuHlfPe71HKWyktoGbFEYqUqptKxtKEwjgFXfnpNrbUiOhfHree7nA3AaMs1SP9SjCTRaBp4B4OF3Nc71wyBiJwAxeNLkUh3Ou121ldjEswlz/G/D9fcXGJ/ZBICdDy9iz/k92coim9vfgXTbxhwevLkpxI1m5FYjTbfqTjg5OLLu/v/x9d8nueXmq9XpibRJvFvtA1q3u4mBZ9ezIvk0nu1GQ1SVnMdQpfCo5sHHbR6nZaOSN+57BPTIcf7Xs72JTrCob/7tLDL2Z+KNzxF83S10q98tR7znuz5P/0b9iUuNY+cmrTBsobFnYxq6N+TuIGM5m6ODI891ea7oCT15CBLydk+VJbYojD9FJGtfjH7AZGB52YpVARj2KWx8G/xuYPTgIA78uJf2Da/ONKFBp7xxuk6C1YYl+Ffvvp79Z+by2uDWsPRB0pvdRvUbRuIgDkQFjuTIm1tof10t3Kq5cVP9m/Km1elB22UVgUHz+WKgMXDp4mwMMdWv5cozt7XIEfSnzB6sSOnC2x6Neb36MF777QloXIvWdVpzIOaA7Xlq7EKXdkWoF0XAz9MVP0/rLVgAZwdnejbomcfd1cmVHgE9WH6s6r8SSovZPWbTqk7uzUWLgYe/8StHbFEY0zC2St0HPIxh++njshSqQhDYy/gBrfw9WPZo98LjmK2wSKcGjOnaEDBmN/GfxThbBPOt48Xbw2/gluZ18yRRGGueuJnTF62v+chSFAUx4HpfVu6PQhAcOozDITMdOtzHvCsJ9FrSq8jyaDSW9GzQk9DTofYWo0JTKsrCTtgySyoTWGT+NIVwfuJ+vNw9Cw03tENAsdJv6uNOU59CBjYL4IHujVm5P4pOjb2MwbQuEwDwdvIuJKZGUzg1nGoUHuga5t6WZbtOoqwpaIvWfRQwVqGUalMmElVy6voWPNhtbzo1ql3wYJqmxDze/nF7i1DuZI1xavInqFYQ0zpX7h0aCmph3FFuUmgqDHcF3sUvx/QAeEmoV+OamkSYB58aPkRfLtgm27VIsQa2KxgF7emdbcI818I914LiaSov+8YZiwhndp9JcEiwnaWpfHTz74ZndU/6Nexnb1HKnb7X9WVD4w082eFJnByc+HT/pyw+uNjeYlUo2tVrV3igCk6hpkHMhXs/AB+aTgEYe2NoNBoLAtwDeOPmN3BxcrG3KOWOi5MLc26eg4+bD3Vc6/BMp2cKj3SN4eRQ+b+zbbEl9QjQDUgAY+EecG23uTUaK9R01qsfNVUbWxRGquX+F9fMwj2NpohMvGGivUXQaMoUWxRG7oV733MtLNzTaIpAQM2Aa7IrSmMbvRv0trcIpYItCmMacJ6cC/eeL0mmIvKmiBwWkb0i8pO5Z7i1cOEisk9EdovIjpLkqdGUJW3q6lnmuflvp7z7f1yrVJX6YYvCcAU+VUoNV0oNAz413UrCGuB6cy3Hv8CzBYTtpZRqq5TqWMI8NUXg50F6XkNRmHHTDHuLUOEY02qMvUWoMNzX+j57i1Aq2KIw/iCngnAF1pYkU6XU70qpdPN0K8bMK00FIrBWIMsGLbO3GJWGLJPlGk1uqjlUw9GhamwhZMs8LxelVFLWiVIqydxEqbR4APguHz8F/C4iCvhQKfVRfomIyASMfTvw8fHJNuVdVJKSkoodt6JT1LJFXcm7M5smL67iWqZ1pirXyWuBwbUGV5n7Z4vCuCQi7ZVSuwBEpAOQXFgkEVkL+Frxmq6UWmaGmQ6kA1/lk0x3pVSEuf/GGhE5rJTaYC2gqUw+AujYsaPq2bNnYSJaJTQ0lOLGregUtWypGanM/HJm2QlURfBx9ynTOlOp62SIvQWwP8Etg+kZ1NPeYpQKtiiMx4HvReQsxmZfvsA9BUcBpVTfgvxF5D4M8yN9VNZWYnnTiDD/z4nITxgbOVlVGJrSp7pj9cIDaTQFcH2d69kfs9/eYmhKCVus1W4XkRaQvanTEaXUlZJkKiL9gf8CtyilrNrqFhE3wEEplWge3wq8UpJ8NZqyQBveyx99baoWtpgGeQRwU0rtV0rtB2qKyOQS5jsPcMfoZtotIgvNvPxFZIUZxgfYJCJ7gG3Ab0qpVSXMV6PRaMoVVYXWOdvSJTVeKTU/60QpddG0L7WguJkqpYLycT8LDDSPjwM3FDcPjUajqQg09GhobxFKDVum1TqKRbtSRByBamUnkkajqSpM7zLd3iLYnWDvqmP52RaFsQr4TkT6iEgf4BvTTXMN4FPDx94iVHi6+XeztwgVltbere0tgt2pClZqs7ClJFMx1jdMMs/XoLdrvWb45e5f6PJ1F3uLUWFZM2wN3q56e1uNdVrWbmlvEUqVQlsYSqlMpdRCpdQw0zTIQeD9shdNUxGo4az3aC4IeScozQAAD/NJREFUXzffKvUFWRbUdqltbxHsxrw+8+wtQqliS5cUItJOROaISDjG1NbDZSqVRqOpMvw2+Dd7i2AX5jSYU+W2681XYYhIMxF5SUQOY7QoTgOilOqllNItDI1GYxM1q9VkWLNhxY7fuo5t4yD5bYH6eo/Xi533uFbjso8X9l1YpLiuDlXPvlhBLYzDQG/gDqVUd1NJZJSPWBqNpirxSNtHih3X29Wb+jXrFxhm/X/Ws3jA1T3EezXoBcC0ztO4o8kd/HjXjzblNbrl6OzjF7q+wMiWIwHwqu5Ft/p5JzcMaTrEpnSrCgUpjCFAJLBeRBaZM6T0sk2NxqSqbIpTHjhIwb3fo1qMYteYXVb9FIp8rAcB8H/t/y/PxAOPah4ANPJoBEBTr6b8NfKvbP+VQ1bi5+aXI041h2pM7TyVm/xvAqBtvbbUr1mfd3q+w4+D8iqcBu4NeLz949nnXXyr/uSQfEfrlFI/Az+bZjkGYdiUqiciHwA/KaV+LycZNZoKiTZ7YTu1XWoz48YZhCeE06tBL8atutrVs2/cvuzjz277DIWik28n/jz9J4+uezRPWu/3fp/zyefpUb8Hq8NXM7LFyDxhpnWeRvPazbNf/mB0jWUR4B7A78N+Z3vUdgCiLkVld2nN7jGb1eGraebVDIB+DftZLdOKIStyyHw26Sx/R/1t8zWpjNhiS+oS8DXwtYh4AcMxptpqhXGNMKDxAFaeWGlvMSocWd0eGtsY2mxo9vHktpNZsDuvsYiOvlf3Sctaw3Fvi3t5+a+XAfig7wd0r989O8y41uOwRs1qNW3awKmTb6c8bl4uXoxoMcJqePdq7iSmJVqVeVnY1f1j2tdrX2jelZEizQdUSl3EMCGe774UmqqHs4OzvUWokNwZeKe9Rai0OErhGwp5u3pntz5a1mnJ2UtnbVrX4O7snq/fN7d/w+5zu20XNBdbRm7h57CfrS5ozWpx3tnkTmb1mFVl9sCwRE8g1xRK7r5ejUFh/fKa/Onq15X3/7F9suWs7rP49+K/1HGtU2C4H+78ocAw13tfz/Xe19ucrzXuDrq7QP+qZGwwN7rGawrl4RsetrcImipGm7ptmN5lOoNqDbIpfA3nGrSt17bQcM1rN7fbynsx5wRphaG5pnF2cNY2pTSlzogWI+jrWeA+a5oKhlYYGpt4tdur9hZBo6nQZO1QWcOp6prT0WMYGpu40f9Ge4ug0VRo+lzXhyntpjCq5Sh7i1JmaIWh0Wg0pYCjgyPj24y3txhliu6S0mg0Go1NaIWh0Wg0Gpuwi8IQkRkiEiEiu83fwHzC9ReRIyISJiLTyltOjSY/XrnpFXuLoNGUO/Ycw3hXKfVWfp7m3uHzgf9v7+5jrKrvPI6/P8wMnWGQp8JMhwflebpBBGQKUgQHqYpQrVBkYc0u7aYdSSS7/aPZVkm7dk0TtbrZaLPdtSsJ9kHadBdrWpVqy6xs//ChrHV8qCshNGVqZRsjdrJN68N3/7hnppfhznBmuHPPvXc+r4TMOb9z7r3fX87lfnJ+99zfuQI4ATwj6eGIeKlUBZoNZPOCzVmXYFZy5TwktQI4GhHHIuKPwH5ykyCamVkGsgyM3ZKel7Q3mdSwvxnkbtrU60TSZmZmGRixISlJTwAfKLBpD/A14DYgkr93A399jq/XAXQANDc3D3vir56enqqcNAzOvW8z6mbQ/XZ38QqqYKV8j/g9WZmqsW8jFhgRkeo3/5K+DvygwKZuYFbe+sykbaDX65tFt62tLdrb21PXmq+zs5PhPrbcnWvf9j+xn+5uBwZQ0veI35OVqRr7ltVVUvnTn24GXiiw2zPAAklzJI0FtgMPl6I+K2zZtML3TB5ttrcWvleCWbXL6iqpOyUtJTckdRy4EUDSdODfImJjRLwjaTdwEKgB9kbEixnVa8CSpiVZl1AWVrZU/604zQrJJDAiouCtsCLi18DGvPVHgEdKVZcN7pKWS7IuoSx45l4brcr5slqzsjR30tysSzDLhAPDzMxScWCYmVkqDgyzIeq9FafZaOPAMBuicXXVe0c1s8E4MMzMLBUHhpmZpeLAMDOzVBwYZmaWigPDbAiWNy/PugSzzDgwzIZgy4ItWZdglhkHhpmZpeLAsCGZO3F0z6NUo5qsSzDLjAPDhmT3st1Zl5Cp1dNXZ12CWWYcGDYk62aty7qETEmeFsRGLweGDUntmKzuuWVmWXNgmJlZKg4MsyEYXzc+6xLMMuPAMBuCmjG+SspGLweGmZmlksk3mJK+A7Qmq5OANyNiaYH9jgO/A94F3omItpIVaWZmp8kkMCLiz3uXJd0NnBpk93UR8duRr8rMzAaT6TWSyl3Uvg24PMs6zMzs7LK+qH4N8HpEvDrA9gB+JCmAf42I+wZ6IkkdQAdAc3MznZ2dwyqop6dn2I8td8Xq24SaCbz17lvnXlAFyuK94fdkZarGvo1YYEh6AvhAgU17IuL7yfIO4MFBnubSiOiW1AQ8LukXEfFkoR2TMLkPoK2tLdrb24dVd2dnJ8N9bLkrVt+mPjSVt06NzsDI4r3h92Rlqsa+jVhgRMRHBtsuqRbYAgx4g4GI6E7+npR0AFgBFAwMK52VLSs5dupY1mWYWYlleVntR4BfRMSJQhslNUo6r3cZuBJ4oYT12QCuuOCKrEswswxkGRjb6TccJWm6pEeS1WbgvyT9HHga+GFEPFbiGq2AhZMXZl1CUWycszHrEswqSmaBERGfiIh/6df264jYmCwfi4glyb9FEfHlbCq1/hpqG85oW3/++pLW0Dq59ew7ncUda+8oQiVmo4d/6W1DNrZmLEf+8giXzbwMgKZxTdy2+jamNUwb8DG3r7m9qDXsu3ofWxduHXD7sqZlp63PmziPey+/t2/9mrnXDPk1185cO+THmFWTrC+rtQpVN6aOr67/6mltP9n2ExbvW1xw/01zN3HoV4c4ePwgAFfNvqpv+Wy+cfU3qBtTx41P3MipP5ziC5d8gca6xoJ3v1t//npuXnEzzY3NHDx+kOXNyzl84jAb5mygobaBgx8/yCtvvMK689ed8Ro3H76ZBzc9yJrvrClYx4JJC1LVa1atfIZhRdXS2ALAriW7AJhSP4Wn/uIpAO667K6+/fKXBzNv4jyWNi1l0dRFfW1XXnAlANtat52x/56Ve2hubAZyoTS1YSqbF2zuG0abPn76GWEBsLRpKY9+/FEm1U+ivqYegAeufoCf7vhp3z7Xzr82Vc1m1cpnGDYirpt/HXVj6ti6cCvj6sb1tTfUNvCpxZ8qymssnLyQQ9sOse6761jRuIL7t95flOeded5Mjr55lMa6RiaMncAtK29hVcsqZk+cXZTnN6tUPsOwovrkhZ8EcmcWHRd1MKV+ymnbn77haTou6uhb37VkFweuPXDaPl9c9cWCz/3pxZ8GoLGusa9tasNUvnfN99j+/u1Fqb+QHR/c4bAww4FhRbbjgzvo2tlV8Eqq/rp2dnHT0puYP3l+X9vFTRdz/cLrC+6/c9FOunZ2UVdTd1p765RW6lRX8DFn86UPf+mM4bHeL9ObxzUP6znNqpUDw8rKPZffA/zpbGKkbVmwhatmX3Va2w1/dgNdO7uY+L6JJanBrFI4MKys9H5I+0d1ZuXHgWFlaeZ5M2ka18RnP/TZrEsxs4SvkrKyVF9bz4+v/3HWZZhZHgeGlYWvrP0KE8ZOyLoMMxuEA8PKwoY5G7IuwczOwt9hmJlZKg4MMzNLxYFhZmapODDMzCwVB4aZmaXiwDAzs1QcGGZmlooDw8zMUlFEZF1D0Un6X+CXw3z4VOC3RSynnLhvlcl9q0yV1rcLImLaYDtUZWCcC0nPRkRb1nWMBPetMrlvlaka++YhKTMzS8WBYWZmqTgwznRf1gWMIPetMrlvlanq+ubvMMzMLBWfYZiZWSoODDMzS8WBkZC0QdIrko5K+nzW9RSTpOOSuiQ9J+nZrOs5F5L2Sjop6YW8timSHpf0avJ3cpY1DtcAfbtVUndy7J6TtDHLGodL0ixJhyS9JOlFSX+btFf8sRukb1Vx7PL5OwxAUg3wP8AVwAngGWBHRLyUaWFFIuk40BYRlfQjooIkrQV6gAci4sKk7U7gjYi4PQn7yRHxuSzrHI4B+nYr0BMRd2VZ27mS1AK0RMQRSecBPwOuAz5BhR+7Qfq2jSo4dvl8hpGzAjgaEcci4o/AfuBjGddkBUTEk8Ab/Zo/BuxLlveR+89acQboW1WIiNci4kiy/DvgZWAGVXDsBulb1XFg5MwAfpW3foLqOuAB/EjSzyR1ZF3MCGiOiNeS5d8AzVkWMwJ2S3o+GbKquCGb/iTNBpYBT1Flx65f36DKjp0DY3S4NCIuBq4GbkqGPqpS5MZYq2mc9WvAPGAp8Bpwd7blnBtJ44F/Bz4TEW/lb6v0Y1egb1V17MCB0asbmJW3PjNpqwoR0Z38PQkcIDcEV01eT8aRe8eTT2ZcT9FExOsR8W5EvAd8nQo+dpLqyH2gfisi/iNpropjV6hv1XTsejkwcp4BFkiaI2kssB14OOOaikJSY/JFHJIagSuBFwZ/VMV5GNiZLO8Evp9hLUXV+2Ga2EyFHjtJAu4HXo6If8zbVPHHbqC+Vcuxy+erpBLJJW//BNQAeyPiyxmXVBSS5pI7qwCoBb5dyX2T9CDQTm7q6NeBvwceAr4LnE9uWvttEVFxXx4P0Ld2ckMaARwHbswb868Yki4FDgNdwHtJ8y3kxvor+tgN0rcdVMGxy+fAMDOzVDwkZWZmqTgwzMwsFQeGmZml4sAwM7NUHBhmZpaKA8NsEJLenzfb6G/yZh/tkfTPI/San5H0V4Ns/6ikfxiJ1zYbjC+rNUupFDPHSqoFjgAXR8Q7A+yjZJ/VEfF/I1WLWX8+wzAbBkntkn6QLN8qaZ+kw5J+KWmLpDuTe5A8lkwbgaTlkv4zmQTyYL9fAve6HDjSGxaS/ia5z8LzkvZD35xLncBHS9JZs4QDw6w45pH7sL8W+CZwKCIWA78HNiWhcS+wNSKWA3uBQr+4X03ufgq9Pg8si4iLgF157c8Ca4reC7NB1GZdgFmVeDQi3pbURW56mceS9i5gNtAKXAg8nhtRoobcDKb9tZC7n0Kv54FvSXqI3BQovU4C04vZAbOzcWCYFccfACLiPUlvx5++HHyP3P8zAS9GxKqzPM/vgfq89U3AWuAaYI+kxclwVX2yr1nJeEjKrDReAaZJWgW56bAlLSqw38vA/GSfMcCsiDgEfA6YCIxP9ltIFcx+apXFgWFWAsmtf7cCd0j6OfAc8OECuz5K7owCcsNW30yGuf4buCci3ky2rQN+OLJVm53Ol9WalRlJB4C/i4hXB9jeTG6a+vWlrcxGOweGWZmR1EruXtdPDrD9Q8DbEfFcaSuz0c6BYWZmqfg7DDMzS8WBYWZmqTgwzMwsFQeGmZml4sAwM7NU/h8v9glw8CWfOQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots() #create an empty plot\n", "\n", "df.plot(ax=ax) #use the dataframe to add plot data, tell it to add to the already created axes\n", "\n", "ax.set(xlabel='Time (s)',\n", " ylabel='Acceleration (g)',\n", " title=filename)\n", "ax.grid() #turn on gridlines\n", "\n", "fig.savefig('full-time-history.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "xme7PeHrWzU8", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Plot with Plotly\n", "Matplotlib may be familiar, but Plotly offers much more interactivity directly in a browser. They also have really [good documentation with a ton of examples online](https://plotly.com/python/).\n", "\n", "The trouble is that plotting too many data points may be sluggish. I've long maintained that plotting 10s of thousands of data points isn't very useful anyways, you just get a shaded mess. \n", "\n", "So here we'll plot:\n", "\n", "- The moving peak value\n", "- The moving RMS\n", "- The time history around the peak" ] }, { "cell_type": "markdown", "metadata": { "id": "-qamxSVYbHLz", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Moving Peak\n", "This takes advantage of Pandas [rolling() function](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.rolling.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 447 }, "id": "03td9T_yeIFn", "outputId": "9d282cb6-2681-4b6a-a12f-91262d836c1f", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\"X (2000g)\"\"Y (2000g)\"\"Z (2000g)\"
Time
0.004394NaNNaNNaN
0.2812030.9765771.1596860.976577
0.5580251.1596861.4038301.068132
0.8348701.4038301.1596861.129168
1.1116891.2207221.2817581.220722
............
26.5768631.1596861.0376130.793469
26.8536550.9765771.0376130.854505
27.1304631.0681321.0681320.946059
27.4072600.9765771.2817580.793469
27.6840641.0986501.1291680.915541
\n", "

101 rows × 3 columns

\n", "
" ], "text/plain": [ " \"X (2000g)\" \"Y (2000g)\" \"Z (2000g)\"\n", "Time \n", "0.004394 NaN NaN NaN\n", "0.281203 0.976577 1.159686 0.976577\n", "0.558025 1.159686 1.403830 1.068132\n", "0.834870 1.403830 1.159686 1.129168\n", "1.111689 1.220722 1.281758 1.220722\n", "... ... ... ...\n", "26.576863 1.159686 1.037613 0.793469\n", "26.853655 0.976577 1.037613 0.854505\n", "27.130463 1.068132 1.068132 0.946059\n", "27.407260 0.976577 1.281758 0.793469\n", "27.684064 1.098650 1.129168 0.915541\n", "\n", "[101 rows x 3 columns]" ] }, "execution_count": 97, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n_steps = 100 #number of points to plot\n", "n = int(df.shape[0]/n_steps) #number of data points to use in windowing\n", "df_rolling_peak = df.abs().rolling(n).max().iloc[::n] #finds the absolute value of every datapoint, then does a rolling maximum of the defined window size, then subsamples every nth point\n", "\n", "df_rolling_peak" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "6uubzmj7X8-V", "outputId": "3050a9a0-88ab-4451-9de6-67d82648bbdd", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = xp.line(df_rolling_peak)\n", "fig.update_layout(\n", " title=\"Rolling Peak\",\n", " xaxis_title=\"Time (s)\",\n", " yaxis_title=\"Acceleration (g)\",\n", ")\n", "fig.show()\n", "fig.write_html('rolling_peak.html',full_html=False,include_plotlyjs='cdn')\n" ] }, { "cell_type": "markdown", "metadata": { "id": "QdJZtyIpaImV", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Moving RMS\n", "Now we'll plot the rolling RMS using the standard deviation. Notice that these rolling value plots make it much easier to compare the datasets than by trying to plot all values which result in a shaded mess.\n", "\n", "Also in this example I'm showing how easy it is to [change the theme of the plotly figure, see their documentation](https://plotly.com/python/templates/) for more examples and information. You can also make custom themes." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "emiCK-7vYJP0", "outputId": "d9aa4fd8-b9ad-4eca-db73-93c962f4ce44", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df_rolling_rms = df.rolling(n).std().iloc[::n] #does a rolling standard deviation of the defined window size, then subsamples every nth point\n", "\n", "fig = xp.line(df_rolling_rms)\n", "fig.update_layout(\n", " title=\"Rolling RMS\",\n", " xaxis_title=\"Time (s)\",\n", " yaxis_title=\"Acceleration (g)\",\n", " template=\"plotly_dark\"\n", ")\n", "fig.show()\n", "fig.write_html('rolling_rms.html',full_html=False,include_plotlyjs='cdn')" ] }, { "cell_type": "markdown", "metadata": { "id": "naRQofdha_3M", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Time History Around Peak\n", "Now let's find the time that had the maximum value and display the time history around that." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "EiNO0_vDKKsy", "outputId": "21220e37-f7b0-4aa7-e0ec-12183c3390cf", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/plain": [ "8.659936" ] }, "execution_count": 100, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.abs().max(axis=1).idxmax()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "vbsRy0N0aY1E", "outputId": "bfdafee7-7c6e-4ee5-cedf-a1fe096580c7", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "peak_time = df.abs().max(axis=1).idxmax() #get the time at which the peak value occurs\n", "d_t = (df.index[-1]-df.index[0])/(len(df.index)-1) #find the average time step\n", "fs = 1/d_t #find the sampling rate\n", "\n", "num = 1000 / 2 #total number of datapoints to plot (divide by 2 because it will be two sided)\n", "df_peak = df[peak_time - num / fs : peak_time + num / fs ] #segment the dataframe to be around that peak value\n", "\n", "fig = xp.line(df_peak)\n", "fig.update_layout(\n", " title=\"Time History around Peak\",\n", " xaxis_title=\"Time (s)\",\n", " yaxis_title=\"Acceleration (g)\",\n", " template=\"plotly_white\"\n", ")\n", "fig.show()\n", "fig.write_html('time_history_peak.html',full_html=False,include_plotlyjs='cdn')\n" ] }, { "cell_type": "markdown", "metadata": { "id": "f7KTFZNijXgI", "pycharm": { "name": "#%% md\n" } }, "source": [ "## PSD\n", "Now using SciPy we can easily compute and plot a PSD using a custom function we'll make to ease the interface to [SciPy's Welch function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.welch.html). This is very similar to [MATLAB's version](https://www.mathworks.com/help/signal/ref/pwelch.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Ihz-UEolbyNT", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def get_psd(df, bin_width=1.0, window=\"hann\"):\n", " d_t = (df.index[-1]-df.index[0])/(len(df.index)-1)\n", " fs = 1/d_t\n", " f, psd = signal.welch(\n", " df.values, fs=fs, nperseg= int(fs / bin_width), window=window, axis=0\n", " )\n", "\n", " df_psd = pd.DataFrame(psd, columns=df.columns)\n", " df_psd[\"Frequency (Hz)\"] = f\n", " df_psd = df_psd.set_index(\"Frequency (Hz)\")\n", " return df_psd\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 447 }, "id": "9iXUTdMfodir", "outputId": "b5e8f5d9-e5a0-4e18-bd01-5ae033d62d76", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\"X (2000g)\"\"Y (2000g)\"\"Z (2000g)\"
Frequency (Hz)
0.0000000.0000480.0000490.000072
4.0000480.0002940.0002760.000332
8.0000950.0002560.0002540.000287
12.0001430.0001890.0002060.000230
16.0001910.0001700.0001560.000193
............
2484.0296060.0000070.0000060.000007
2488.0296540.0000070.0000060.000006
2492.0297020.0000070.0000070.000007
2496.0297490.0000070.0000070.000007
2500.0297970.0000030.0000040.000004
\n", "

626 rows × 3 columns

\n", "
" ], "text/plain": [ " \"X (2000g)\" \"Y (2000g)\" \"Z (2000g)\"\n", "Frequency (Hz) \n", "0.000000 0.000048 0.000049 0.000072\n", "4.000048 0.000294 0.000276 0.000332\n", "8.000095 0.000256 0.000254 0.000287\n", "12.000143 0.000189 0.000206 0.000230\n", "16.000191 0.000170 0.000156 0.000193\n", "... ... ... ...\n", "2484.029606 0.000007 0.000006 0.000007\n", "2488.029654 0.000007 0.000006 0.000006\n", "2492.029702 0.000007 0.000007 0.000007\n", "2496.029749 0.000007 0.000007 0.000007\n", "2500.029797 0.000003 0.000004 0.000004\n", "\n", "[626 rows x 3 columns]" ] }, "execution_count": 121, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df_psd = get_psd(df,bin_width=4) #compute a PSD with a 1 Hz bin width\n", "df_psd.to_csv('psd.csv') #save to a CSV file\n", "df_psd" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "jE2ZLBW-aeZr", "outputId": "2ff24aca-75ef-4434-c45d-37e988e426ce", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = xp.line(df_psd)\n", "fig.update_layout(\n", " title=\"Power Spectral Density (PSD)\",\n", " xaxis_title=\"Frequency (Hz)\",\n", " yaxis_title=\"Acceleration (g^2/Hz)\",\n", " xaxis_type=\"log\",\n", " yaxis_type=\"log\"\n", ")\n", "fig.show()\n", "fig.write_html('psd.html',full_html=False,include_plotlyjs='cdn')" ] }, { "cell_type": "markdown", "metadata": { "id": "xA_ZSrAM48Sp", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Cumulative RMS from PSD\n", "Now that we have the PSD, we can easily compute and plot the overall RMS value. This is partially thanks to the [cumulative sum function](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.cumsum.html) in Pandas. \n", "\n", "The nice thing about a PSD (in addition to the easy control of the bin width) is that the area directly relates to the RMS level in the time domain. The equation is as follows.\n", "\n", "$$ g_{\\text{RMS}}=\\sqrt{\\int \\text{PSD}(f)\\ df} $$\n", "\n", "Let's demonstrate by quickly using the PSD just calculated, integrating, and taking the square root and compare to the values we calculated from the time domain." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "MIikvleC5BIe", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def rms_from_psd(df_psd):\n", " d_f = df_psd.index[1] - df_psd.index[0]\n", " df_rms = df_psd.copy()\n", " df_rms = df_rms*d_f\n", " df_rms = df_rms.cumsum()\n", " return(df_rms**0.5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "5w283bcf5MVd", "outputId": "f4b6f52c-54fd-445f-d135-7bdf634e9360", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df_rms = rms_from_psd(df_psd)\n", "\n", "fig = xp.line(df_rms)\n", "fig.update_layout(\n", " title=\"Cumulative RMS\",\n", " xaxis_title=\"Frequency (Hz)\",\n", " yaxis_title=\"Acceleration (g RMS)\",\n", " xaxis_type=\"log\",\n", " #yaxis_type=\"log\"\n", ")\n", "fig.show()\n", "fig.write_html('cum_rms.html',full_html=False,include_plotlyjs='cdn')" ] }, { "cell_type": "markdown", "metadata": { "id": "tqZQytK3jLju", "pycharm": { "name": "#%% md\n" } }, "source": [ "## FFT" ] }, { "cell_type": "markdown", "metadata": { "id": "PQ6Asg0zjOWn", "pycharm": { "name": "#%% md\n" } }, "source": [ "### Typical FFT (Or Should We Say DFT)\n", "This uses SciPy's discrete [Fourier transform function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fft.html). The trouble here is that this may be very long and therefore plotting a LOT of data." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "rDTAargQjSNz", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "from scipy.fft import fft, fftfreq\n", "\n", "def get_fft(df):\n", " N=len(df)\n", " fs = len(df)/(df.index[-1]-df.index[0])\n", " \n", " x_plot= fftfreq(N, 1/fs)[:N//2]\n", " \n", " df_fft = pd.DataFrame()\n", " df_phase = pd.DataFrame()\n", " for name in df.columns:\n", " yf = fft(df[name].values) \n", " y_plot= 2.0/N * np.abs(yf[0:N//2])\n", " \n", " phase = np.unwrap(2 * np.angle(yf)) / 2 * 180/np.pi\n", " phase = phase[0:N//2]\n", " \n", " df_phase = pd.concat([df_phase,\n", " pd.DataFrame({'Frequency (Hz)':x_plot[1:],\n", " name:phase[1:]}).set_index('Frequency (Hz)')],axis=1)\n", " df_fft = pd.concat([df_fft,\n", " pd.DataFrame({'Frequency (Hz)':x_plot[1:],\n", " name:y_plot[1:]}).set_index('Frequency (Hz)')],axis=1)\n", " \n", " return df_fft, df_phase" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "3qE8BdHX1pGu", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "df_fft, df_phase = get_fft(df)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "id": "sEuX5lKd0cWo", "outputId": "0f55d51f-529f-41c3-ff42-67adafb5217a", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXwV9b3/8dc7YQk7CBLWChJQQahCRHvdcAEVrRaXiq2Iu6LWq629Ymtba9WitV7vvWIVqz+graKVVtGqiEtqXZDFiyKgEhCvCRE0GCBICEk+vz9mEifhnJOTmJOE5PN8PM7jzPKd+X6/M3POZ9bvyMxwzjnn4klr6gI455xr3jxQOOecS8gDhXPOuYQ8UDjnnEvIA4VzzrmEPFA455xLqEUECkkbJJ3Y1OWojaTnJU1NMm2mpNckbZf0+1SXLUE5bpH056bKv6FJypF0aT2nTbgsJB0gaUW4zq6tfylBkknK+ibzaCzRskp6QNIvwu5xkvJSnPcPJb2Yyjzi5FvvukkaFC6zNg1drlRpEYEinlgrM9V/fJLaSfpCUuea/WZ2ipnNSXJWlwNfAF3N7CepKm9rUHOdpNB/AK+aWRcz++8kynWTpDu+SYb1/dOR1C/625D0A0nLJBVLKgh3ao6qa3nM7Eoz+01dp0tGrLqa2V/MbEKK8hsh6UVJWyQVSVouaWIq8mruWnSgaCLHACvMrDhOf7L2A1abPxHZEOq7DupqP2BVHdKfCjyXorLUZiLwAoCkHwP3AncAmcC3gPuBMxqzQJLSGzO/JDwDLAL6AL2Ba4FtTVqipmJme/0H2ADcALwHbAUeBzoBO4EKoDj8/AAoBXaH/e+G0+cAvwWWEGwITwP7hOMygD8DhUARsBTITFCWe4Afx+oP87k07L4QeB24G/gS+Bg4JRw3OyxjaVjOE4H2BD/mjeHnXqB9nDJ0Ax4GCoB84DYgvbZ8w/GDgX8C2wl+JPcBf46M/yvwWbicXwNGRMb1BBaEy3AJ8Bvg9QTL6gxgRZh+HXByOLxfOJ8tQC5wWWSaW8Iy/Dks40pgGHATsBn4FJgQb52E6+A3wBvh9C8CvcJx44C8GNvWiZG8nyTYvrYD7wDfDse9ApQDJeE6G0bwZ7w6TJsP3BCZb4+wvJXr5afh+toIXAwYkBWOOxX433A5fQrcEpnP/4VpK7fx7wBDwvIUEhyV/gXoXqNefwPODLeVYuCcBOtpLPAWwfZfEG4T7SLjo2WdDdwWXZ7Az8JybAB+GJluNvAHgmC5g2A7r2tdLySyjQH/RvAb3Rp+/1tkXNx1H6POvcK8uscZX1m3n4TrsQC4KDI+UT0GhfNuE/afFS6bgwl23qcT/B4KgScI/4vilOMyYE1Yn9XA6HD4jQTb3HbgQ+AEgt/Vzuj8gEPDddM24X9sY/2Zp/ITLuQl4YLYJ1xwVxL7h38LkT++yAaUH66oTsD8yjTAFQR7Fh2BdGAMwekgwhX6bI15fQAcEKufPQPF7nBFpwPTCP4kVPMHF/bfCiwm2LPZF3gT+E2c5fF34MGwLr3DZXNFkvm+RfDH2p5gT3w71QPFxUAXvg5cKyLj5oUbdqdwWeYTJ1AQ/PlsBcaHP47+wIHhuNcI9mgzgEOAz4HjI+uvBDgJaAPMJQh2PwfahvX6ON46CdfBOoI/8g5h/4zojz/GthUNFLuBs8O8bgjzbltz/Yb9BcDRYXcPwh9x2D8ZeCzsPhnYxNfb36NU//MdB4wMl9OoMO33Yv3phMOywuXanmBbeQ24NzK+LcGfQ5cw77Lo9DHW1RjgiHB5DyL4fV0XGZ8oUJTx9fZ0LEFAOCCSditwZFi3jHrU9ULCbYzgt/8lMCUs63lhf8/a1n2MOgtYCzwLfI8aO4eRut0aLs+JwFdAj7qsM+Aigp2hyuX37wS/8wHhMnuQcDuJUcZzCH5jh4XlzSI4qj2AIDj1i+Q3JOx+heo7Xr8DHqj1P7Yx/shT/SH4MZ8f6b8LeIC6BYoZkf7hBHvz6QR/jG8Co5IoxxAgN0F/DtUDRXRcx3Dj6VPzBxf2rwMmRvpPAjbEKEMmsAvoEBl2HsG584T5EpxyKAM6RcY/WnN5RcZ1D6ftFi6r3YR/9uH4O4gfKB4E/jPG8IEEe+ZdIsN+C8yOrL9FkXHfJdi7rNwz70JkTzDOOrg50n8V8ELYHWt72UD1QLE4Mi6N6sGgav2G/f9HsKPRNUY9/wRMCbsfqbH9DSPy5xtj2nsrlx0x/jxjpP8e8L+R/hOAl8PuHwKf1fH3dh3w90h/bYEiuj09AfwiknZuLXklrCvVA8UUYEmN6d8CLqxt3cfJewDB0dM6gjMTrwFDI3XbWaMsm4Ej6lCPGwiOAgZE0q0BToj09yX4Xe2xfoGFwL/HGJ4VluVEahwpAJcCr4TdIggox9S2zlvSNYrPIt1fAXW9cPlppPsTgr2EXgQ/6IXAPEkbJd0lqW2ceUwEnk/QH7fMZvZV2Bmv3P3CckXL2C9Guv3CsheEF+CKCP6UeyeRbz/gSzPbUSMfIDiHLGmGpHWSthH8iUKwnPYl2EOquRzjGUjwA6ypH7DFzLbXmE//SP+mSPdO4AszK4/0V9YHYq+Db7KtVNXPzCoITj/EWg8QnFKYCHwi6Z+SvgMgKY1gj/+FMF0/Eiw3SYdLelXS55K2Ehwt94pXwPCOuXmS8sP19Oca6Sfy9bWRQqBXoovhkoZJelbSZ+H87kiUfw2xtqfo8orWu851raHmb6Qyv+i2E3Pdh3drFYefnwGYWZ6ZXWNmQwh+VzsIjmArFZpZWZz5JVOPnwIzzSx6w81+wN8jv901BDtOmTHqG/M3ZGa5BMH8FmBzuC1ULvP5wHck9SU4Y1AB/CvGvKtpSYEiFktyGAQLvdK3CKL4F2a228x+bWbDCc5/ngZcEGce0R9grP5vYiPBRhQt48YY6T4lOKLoZWbdw09XMxuRRB4FQA9JnWrkU+kHBNcVTiQ4ihgUDhfB6aEy9lyO8XxKsLdf00ZgH0ldaswnP4nyx1KXdbCD4AgLqLq4um+NNAMj49MI9jpjrQfMbKmZnUEQpJ8i2JuG4FTBJ2b2edhfQOLl9ijBNZuBZtaN4GhZldnEyPqOcPhIM+sKnB9JD9WXyVsE28v3YtUh9AeC03dDw/n9rMb8Eom1PUWXV83y17WuUTV/I5X51brtWHC3Vufws8edaGb2KTCT4PRgMhLVo9IE4GZJZ0WGfUpwzbB75JNhZrHqEO83hJk9amZHESwPA+4Mh39JcG3mXILf8zwLDy8SaemBYhPQU1K3GsMGhT/yqPMlDZfUkeC845NmVi7pOEkjwz+NbQQBpKJmRuF0Y4FXY/U3gMcINqp9JfUCfkmwp1iNmRUQbAi/l9RVUpqkIZKOrS0DM/sEWAb8Oryl9CiCUzuVuhD8qRQS/KHeEZm2nOAC6S2SOkoaDkxNkN3DwEWSTgjL2F/SgeEP8k3gt5IyJI0CLolV19rUYx18BGRIOjU8aryZ4Dxx1BhJZ4Z74NcRLI/FMfJuF97j383MdhNsO5XbzUTgH5HkTwAXRra/X9WYXReCo6wSSWMJfuCVPg/nu3+N9MXAVkn9CfZcK8s1mOAmiDUAZraVYFuaKel74bprK+kUSXdF5rcNKJZ0IMF1rbqo3J6OJtjR+muCtHWta9RzwLDwVt82ks4lOI38bB3Li6Qekn4tKSvcPnsRnIbeY13Xox6VVhFcI5op6fRw2APA7ZL2C8uxr6R4d5/9EbhB0hgFsiTtp+B5nuMltSe4nld5U0+lRwl2ds8Ou2vVogOFmX1A8Ae7PjyU68fXG2mhpHciyf9EcM70M4KLapUPTPUhuNNlG8Fh4D/DtEj6maTK0xrHA2+ZWUmc/m/qNoI/8fcI7vR5JxxW+dBR9LbMC4B2BOc/vwzL3zfJfH4AHE5wx9GvqH6oPZfgUD4/nHfNH801BIfenxEsy/8XHSlplaQfApjZEoILef9JcEHzn3y9N3gewdHKRoIL878ys5eSLH9UndZB+Kd5FcEPMJ/gCKPmQ1VPE+yNVV40PTMMBLFMATaEp2uuJLgeADVuizWz5wnOYb9CcGHzlRrzuQq4VdJ2gj/1JyLTfgXcDrwRbuNHAL8GRhMs138QBPBKe9ySa2a/B35MEBg/J9hTvYbgKAiCc+k/ILix4SGCu76S9RnBstpIcPfVleHvMp661jVaj0KCQPQTgp2Z/wBOM7Mv6lDeSqUE2+BLBL/99wl2Ci5Mcvq49ahR5nfDMj8k6RTgvwiORF4Mp11M8HsEIDw1dnQ47V8JlsejBOvmKYIL+u2BGQQ3LHxGcER7UyTbBcBQgmtT7yZTGSVx1NHiScohuGD7x28wj/uB983s/lj9rZGkCwku7tb5wa0Gyr/ZrQNJmQS3TfZP5pA/Bfk/B9xnZk31/IbbC+01j5DvBVYQ3EYbr981vua4DroBP2mKIBHKoeFOh7pWwgNFAzGzWYn6XeNrjuvAzD4iuBbSVPnfVXsq56rzU0/OOecSatEXs51zzn1zLerUU69evWzQoEH1mnbHjh106tSp9oQtiNe55Wtt9QWvc10tX778CzOr+bxQNS0qUAwaNIhly5bVa9qcnBzGjRvXsAVq5rzOLV9rqy94netKUqIWFAA/9eScc64WHiicc84llLJTT5IeIXjicLOZ7dE+iqSf8vWTqm2Ag4B9zWyLpA0ETxqWA2Vmlp2qcjrnnEssldcoZhM00Ts31kgz+x1BW+hI+i5wvZltiSQ5rp6P3jvnmqHdu3eTl5dHSUlDtWoTW7du3VizZk1K82hukqlzRkYGAwYMoG3beI1fx5eyQGFmr0kalGTy8wjaZHLOtVB5eXl06dKFQYMGISXb+Gzdbd++nS5dutSesAWprc5mRmFhIXl5eQwePLjO80/pA3dhoHg21qmnSJqOBA2vZVUeUUj6mKAhMQMeTPSEraTLgcsBMjMzx8ybN69eZS0uLqZz57q+wmLv5nVu+ZpTfbt168aQIUNSGiQAysvLSU9vbq/fTq1k6mxmrFu3jq1bt1Ybftxxxy2v7fR+c7g99rvAGzVOOx1lZvmSegOLJH1gZq/FmjgMIrMAsrOzrb63iPktda1Da6tzc6rvmjVr6Nq1a8rz8SOK+DIyMjj00EPrPP/mcNfTZGqcdqp8SYeZbSZoZnpsE5Qrrnc/f5cPtiRqKdk551qOJg0U4QuFjiVo479yWKfKt5uFb8aaQNAWfLNx/nPnc84z5zR1MZxz9VTZgkPl97JlyxgxYgSlpaUArFu3jv33359t27btMW1BQQGnnXYaAIsWLWLMmDGMHDmSMWPG8MorX79KZPny5YwcOZKsrCyuvfbayndWs2XLFsaPH8/QoUMZP348X375JRCcGrr22mvJyspi1KhRvPPOO3vkXdPkyZNZu3Zt3Ho1lJQFCkmPEbxm8QBJeZIukXSlpCsjySYBL9Z4p24m8Lqkd4ElwD/M7AWccy5FsrOzOfbYY7n77rsBuPrqq7n99ttjniq75557uOyyywDo1asXzzzzDCtXrmTOnDlMmTKlKt20adN46KGHWLt2LWvXruWFF4K/sRkzZnDCCSewdu1aTjjhBGbMmAHA888/X5V21qxZTJtW+4sEp02bxl13pb5B4FTe9XReEmlmE9xGGx22Hvh2akrlnHOw7777VvsGuOOOOzj00ENp06YNZWVlnHde7L+w+fPnc9tttwFUO98/YsQIdu7cya5du9iyZQvbtm3jiCOCl/BdcMEFPPXUU5xyyik8/fTT5OTkADB16lTGjRvHnXfeydNPP80FF1yAJI444giKioooKCggMzOTa665hldeeYWBAwfStm1bLr74Ys4++2yOPvpoLrzwQsrKyuLWqyE0h4vZzrlW5tfPrGL1xj1P63wTw/t15VffHZFU2qVLl1b7BujevTvTp0/nqquuYvXq1TGn+/jjj+nRowft29d8lXoQQEaPHk379u3Jz89nwIABVeMGDBhAfn4+AJs2baJv3+DNxH369GHTpk0A5OfnM3DgwD2meeONN9iwYQOrV69m8+bNHHTQQVx88cUApKWlkZWVxcqVKznmmGNi1qshNIeL2c451yw8//zzZGZmxg0UBQUFMffWV61axY033siDDz5Yp/wk1Xq78Ouvv84555xDWloaffr04bjjjqs2vnfv3nz22Wd1yreu/IjCOdfokt3zb0zPPvssW7duZeHChUyaNImTTjqJjh07VkvToUOHPZ4sz8vLY9KkScydO5chQ4YA0L9/f/Ly8qql6d+/PwCZmZkUFBTQt29fCgoK6N27d9U0n376acxpEikpKSEjI6N+lU6SH1E451q9nTt38uMf/5iZM2cycuRIzjjjDG6//fY90g0bNowNGzZU9RcVFXHqqacyY8YMjjzyyKrhffv2pWvXrixevBgzY+7cuZxxxhkAnH766cyZMweAOXPmVBs+d+5czIzFixfTrVs3+vbty5FHHsn8+fOpqKhg06ZNVdc3Kn300UcMHz68gZdIdR4onHOt3m9+8xsmTZpU9Yd7yy238Nhjj1W79RSgU6dODBkyhNzcXADuu+8+cnNzufXWWznkkEM45JBD2Lx5MwD3338/l156KVlZWQwZMoRTTjkFgOnTp7No0SKGDh3KSy+9xPTp0wGYOHEi+++/P1lZWVx22WXcf//9AJx11lkMGDCA4cOHc/755zN69Gi6desGBNc7OnToQGZmZkqXj596cs61enfccUe1/i5durB+/fqYaa+55hpmz57Nbbfdxs0338zNN98cM112djbvv7/nI2A9e/bk5Zdf3mO4JGbOnLnH8LS0NO6++246d+5MYWEhY8eOZeTIkQA8+uijXHHFFbXW75vyQOGcc3UwadIkCgsLGzXP0047jaKiIkpLS/nFL35Bnz59gOBOrSlTprBz586U5u+Bwjnn6ujSSy9t1PxqXpeodNFFFzVK/n6NwjnnXEIeKJxzziXkgcI551xCHiicc84l5IHCOdfqRJvjLikp4cADD2TlypVV43/3u9/FvO10586dHHvssZSXl7NixQq+853vMGLECEaNGsXjjz9ele7jjz/m8MMPJysri3PPPbeq+fJdu3Zx7rnnkpWVxeGHH17t4b3f/va3ZGVlccABB7Bw4cJa63DDDTdUa9Z83LhxbNiwocGbGAcPFM65Vi4jI4N7772Xq666CjMjPz+fBx54oKr576hHHnmEM888k/T0dDp27MjcuXNZtWoVL7zwAtdddx1FRUUA3HjjjVx//fXk5ubSo0cPHn74YQAefvhhevToQW5uLtdffz033ngjAKtXr2bevHlV87rqqqsoLy9PWO4f/ehHMcuYCh4onHOtTs3muE8++WT69u3L3Llzuf7667nlllvo0aPHHtP95S9/qWpyY9iwYQwdOhSAfv360bt3bz7//HPMjFdeeYWzzz4bCJoSf+qppwB4+umnmTp1KgBnn302L7/8MmbG008/zeTJk2nfvj2DBw8mKyuLJUuWAMFT4wcccABHHXUU5513XtU7M/bbbz8KCwurGgTcZ599SE9Pb/AmxsGfo3DONYXnp8NnK2tPVxd9RsIpye1hx2qO+95772Xs2LEMHTq02guIKpWWlrJ+/fqYp3aWLFlCaWkpQ4YMobCwkO7du9OmTfD3Gm1iPNqUeJs2bejWrRuFhYXk5+dXvbsiOs3SpUuZP38+7777Lrt372b06NGMGTOmKt3o0aN54403mDBhAn/729/2qFND8UDhnHMERwXHH3981WtOa/riiy/o3r37HsMLCgqYMmUKc+bMIS2tYU/SvPHGG5xxxhlkZGSQkZHBd7/73Wrje/fuzcaNGxs0z1g8UDjnGl+Se/6NLS0tLe6ffawmxrdt28app57K7bffXnVE0LNnT4qKiigrK6NNmzbVmguvbEp8wIABlJWVsXXrVnr27Bm3ifFoU+WxlJSU0KFDh29S5aT4NQrnnEtCjx49KC8vrwoWpaWlTJo0iQsuuKDqegQEjfsdd9xxPPnkk8CeTYlXNjH+5JNPcvzxxyOJ008/nXnz5rFr1y4+/vhj1q5dy9ixYznyyCN55plnKCkpobi4mGeffbZamT766CMOPvjglNfdA4VzziVpwoQJvP766wA88cQTvPbaa8yePbuqifEVK1YAcOedd3LPPfeQlZVFYWEhl1xyCQCXXHIJhYWFZGVlcc8991TdtTRixAi+//3vM3z4cE4++WRmzpxJeno6hx12GKeffjqjRo3ilFNOYeTIkVVNjO/evZvc3Fyys7NTX3EzS8kHeATYDLwfZ/w4YCuwIvz8MjLuZOBDIBeYnmyeY8aMsfp69dVXk0578OyD7eDZB9c7r+aiLnVuKVpbnZtTfVevXt0o+Wzbti1l816+fLmdf/75KZt/LNu3bzczsx07dtiYMWNs+fLlZmb2t7/9zW6++WYzS77OsdYBsMxq+W9N5TWK2cB9wNwEaf5lZtWuHElKB2YC44E8YKmkBWYW+yW2zjnXSEaPHs1xxx1HeXk56enpjZLn5ZdfzurVqykpKWHq1KmMHj0agLKyMn7yk580ShlSFijM7DVJg+ox6Vgg18zWA0iaB5wBeKBwzjW5iy++uFHze/TRR2MOP+eccxqtDE1919N3JL0LbARuMLNVQH/g00iaPODweDOQdDlwOQQvLY/XbnttiouL6zxtffNqLupT571da6tzc6pvt27d2L59e8rzKS8vb5R8mpNk61xSUlKv7aEpA8U7wH5mVixpIvAUMLSuMzGzWcAsgOzsbBs3bly9CpOTk0PS0wY3LSSfvpmqU51biNZW5+ZU3zVr1tClS5eU57N9+/ZGyac5SbbOGRkZHHrooXWef5Pd9WRm28ysOOx+DmgrqReQDwyMJB0QDnPOOdcEmixQSOojSWH32LAshcBSYKikwZLaAZOBBU1VTueca+1SFigkPQa8BRwgKU/SJZKulHRlmORs4P3wGsV/A5PDu7XKgGuAhcAa4Inw2oVzzjWIaDPjADNnzqx6FuKQQw7h4IMPRhJr1qzZY9qCgoKqZj4WLVrEmDFjGDlyJGPGjKnW7Pfy5csZOXIkWVlZXHvttZW3/rNlyxbGjx/P0KFDGT9+PF9++SUQPKpw7bXXkpWVxahRo3jnnXdqrcfkyZNZu3Zt3Ho1lJQFCjM7z8z6mllbMxtgZg+b2QNm9kA4/j4zG2Fm3zazI8zszci0z5nZMDMbYma3p6qMzjkHcPXVV7NixYqqz+mnn84Pf/hDDjrooD3S3nPPPVx22WUA9OrVi2eeeYaVK1cyZ86cao0JTps2jYceeoi1a9eydu1aXnjhBQBmzJjBCSecwNq1aznhhBOqHrp7/vnnq9LOmjWLadOm1VruadOmcddddzXEIkjIn8x2zrU6NZsZj3rttdd44oknuP/++2NOO3/+fE4++WQADj30UPr16wcET1fv3LmTXbt2UVBQwLZt2zjiiCOQxAUXXBCzqfGaTZBfcMEFSOKII46gqKiIgoICKioquOqqqzjwwAMZP348EydOrGoe5Oijj+all16irKys1np9E019e6xzrhW6c8mdfLDlgwad54H7HMiNY29MKm2sZsYBioqKuPDCC/nTn/5E165d95ju448/pkePHrRv336PcfPnz2f06NG0b9+e/Px8BgwYUDUu2tT4pk2b6Nu3LwB9+vRh06ZNQPUmyKPTvPHGG2zYsIHVq1ezefNmDjrooKpnOdLS0sjKymLlypUcc8wxcev1TXmgcM650JVXXsmUKVM48sgjY44vKCiIube+atUqbrzxRl588cU65SeJ8J6euF5//XXOOecc0tLS6NOnD8cdd1y18b179656eVGqeKBwzjW6ZPf8G9OcOXP45JNP+POf/xw3TaymxvPy8pg0aRJz585lyJAhAHs0ER5tajwzM5OCggL69u1LQUEBvXv3rpomVlPjtSkpKSEjIyP5itaDX6NwzrV669ev52c/+xl/+ctfqt5MF8uwYcPYsGFDVX9RURGnnnoqM2bMqHYU0rdvX7p27crixYsxM+bOnRuzqfGaTZDPnTsXM2Px4sV069aNvn37cuSRRzJ//nwqKirYtGnTHk9Wf/TRRwwfPryBlkRsHiicc63enXfeyVdffcWZZ55Z7TbZf/3rX9XSderUiSFDhpCbmwvAfffdR25uLrfeemvVNJs3bwbg/vvv59JLLyUrK4shQ4ZwyimnADB9+nQWLVrE0KFDeemll5g+fToAEydOZP/99ycrK4vLLrus6mL6WWedxYABAxg+fDjnn38+o0ePrmpqfNOmTXTo0IHMzMyULh8/9eSca/UefPBBHnzwwaTSXnPNNcyePZvbbruNm2++mZtvvjlmuuzsbN5///09hvfs2ZOXX355j+GSmDlz5h7D09LSuPvuu+ncuTOFhYWMHTuWkSNHAkGDgVdccUVS5f4mPFA451wdTJo0icLCwkbN87TTTqOoqIjS0lJ+8Ytf0KdPHwC6d+/OlClT2LlzZ0rz90DhnHN1dOmllzZqfvFafL3ooosaJX+/RuGcazSVzVi4xvdNlr0HCudco8jIyKCwsNCDRRMwMwoLC+t9G62fenLONYoBAwaQl5fH559/ntJ8GuO5guYmmTpnZGRUe1q8LjxQOOcaRdu2bRk8eHDK88nJyanXy3n2Zqmus596cs45l5AHCueccwl5oHDOOZeQBwrnnHMJeaBwzjmXkAcK55xzCaUsUEh6RNJmSXu2ihWM/6Gk9yStlPSmpG9Hxm0Ih6+QtCxVZXTOOVe7VB5RzAZOTjD+Y+BYMxsJ/AaYVWP8cWZ2iJllp6h8zjnnkpCyB+7M7DVJgxKMfzPSuxio3yODzjnnUqq5PJl9CfB8pN+AFyUZ8KCZ1TzaqCLpcuByCF4xGK+VxdoUFxfXedr65tVc1KfOe7vWVufWVl/wOqeEmaXsAwwC3q8lzXHAGqBnZFj/8Ls38C5wTDL5jRkzxurr1VdfTTrtwbMPtoNnH1zvvJqLutS5pWhtdW5t9TXzOtcVsMxq+W9t0rueJI0C/gicYWZVbwIxs/zwezPwd2Bs05TQOedckwUKSd8C/gZMMbOPIsM7SRtZKUUAABo1SURBVOpS2Q1MAGLeOeWccy71UnaNQtJjwDigl6Q84FdAWwAzewD4JdATuF8SQJkFdzhlAn8Ph7UBHjWzF1JVTuecc4ml8q6n82oZfymwx/sEzWw98O09p3DOOdcU/Mls55xzCXmgcM45l1DCU0+SMoDTgKOBfsBOggvL/zCzVakvnnPOuaYWN1BI+jVBkMgB3gY2AxnAMGBGGER+YmbvNUI5nXPONZFERxRLzOxXccbdI6k38K0UlMk551wzEjdQmNk/Ek0YPgy3ucFL5Jxzrlmp9fZYSc8QtL0UtRVYRtAOU0kqCuacc655SOaup/VAMfBQ+NkGbCe4VvFQ6ormnHOuOUjmgbt/M7PDIv3PSFpqZodJ8jufnHOuhUvmiKJz2C4TUNVGU+ewtzQlpXLOOddsJHNE8RPgdUnrAAGDgavCBvvmpLJwzjnnml6tgcLMnpM0FDgwHPRh5AL2vSkrmXPOuWYh7qknSUdVdpvZLjN7N/yUhOO7Sjq4MQrpnHOu6SQ6ojhL0l3AC8By4HOCJ7OzCN5Ktx/BaSnnnHMtWKIH7q6XtA9wFnAO0Jegrac1BM9PvN44RXTOOdeUEl6jMLMtfP38hHPOuVbImxl3zjmXkAcK55xzCXmgcM45l1BS78yW9G/AoGh6M5ubojI555xrRmo9opD0J+Bu4CjgsPCTnczMJT0iabOk9+OMl6T/lpQr6T1JoyPjpkpaG36mJlUb55xzDS6ZI4psYLiZ1WxqPBmzgfuAeEcfpwBDw8/hwB+Aw8Pbcn8V5m3AckkLzOzLepTBOefcN5DMNYr3gT71mbmZvQZsSZDkDGCuBRYD3SX1BU4CFpnZljA4LAJOrk8ZnHPOfTPJHFH0AlZLWgLsqhxoZqc3QP79gU8j/XnhsHjD9yDpcuBygMzMTHJycupVkOLi4jpPW9+8mov61Hlv19rq3NrqC17nVEgmUNySstwbgJnNAmYBZGdn27hx4+o1n5ycHJKeNmwzt755NRd1qnML0drq3NrqC17nVKj11JOZ/RP4AOgSftaEwxpCPjAw0j8gHBZvuHPOuUaWzF1P3weWELT39H3gbUlnN1D+C4ALwrufjgC2mlkBsBCYIKmHpB7AhHCYc865RpbMqaefA4eZ2WYASfsCLwFP1jahpMeAcUAvSXkEdzK1BTCzB4DngIlALvAVcFE4bouk3wBLw1ndGrY75ZxzrpElEyjSKoNEqJAkn+g2s/NqGW/A1XHGPQI8kkw+zjnnUieZQPGCpIXAY2H/uQRHAs4551qBZF6F+lNJZwFHhoNmmdnfU1ss55xzzUVSbT2Z2XxgforL4pxzrhmKGygkvW5mR0naTtCMRtUogssLXVNeOuecc00u0atQjwq/uzRecZxzzjU3ybYeW+sw55xzLVMyt7mOiPZIagOMSU1xnHPONTdxA4Wkm8LrE6MkbQs/24FNwNONVkLnnHNNKm6gMLPfhtcnfmdmXcNPFzPraWY3NWIZnXPONaFknqO4KWxvaSiQERn+WioL5pxzrnmoNVBIuhT4d4IWXFcARwBvAcentmjOOeeag2QuZv87wXuyPzGz44BDgaKUlso551yzkUygKDGzEgBJ7c3sA+CA1BbLOedcc5FMEx55kroDTwGLJH0JfJLaYjnnnGsukrmYPSnsvEXSq0A34IWUlso551yzkTBQSEoHVpnZgVD1WlTnnHOtSMJrFGZWDnwo6VuNVB7nnHPNTDLXKHoAqyQtAXZUDjSz01NWKuecc81GMoHiFykvhXPOuWYrmYvZ/5S0HzDUzF6S1BFIT33RnHPONQfJNDN+GfAk8GA4qD/BrbK1knSypA8l5UqaHmP8f0paEX4+klQUGVceGbcgueo455xraMmceroaGAu8DWBmayX1rm2i8I6pmcB4IA9YKmmBma2uTGNm10fS/4jgqe9KO83skKRq4ZxzLmWSeTJ7l5mVVvaE76OwBOkrjQVyzWx9OP084IwE6c8DHktivs455xpRMkcU/5T0M6CDpPHAVcAzSUzXH/g00p8HHB4rYXgNZDDwSmRwhqRlQBkww8xinu6SdDlwOUBmZiY5OTlJFG1PxcXFdZ62vnk1F/Wp896utdW5tdUXvM6pkEygmA5cAqwErgCeA/7YwOWYDDwZPrdRaT8zy5e0P/CKpJVmtq7mhGY2C5gFkJ2dbePGjatXAXJyckh62jnBV33zai7qVOcWorXVubXVF7zOqZDMXU8VwEPhpy7ygYGR/gHhsFgmE1wLieabH36vl5RDcP1ij0DhnHMuteIGCkkrSXAtwsxG1TLvpcBQSYMJAsRk4Acx8jmQ4KG+tyLDegBfmdkuSb2AI4G7asmvyS3ftJwdu3dwzIBjmroozjnXYBIdUZz2TWZsZmWSrgEWEjx38YiZrZJ0K7DMzCpveZ0MzDOzaFA6CHhQUgXBBfcZ0bulmqsLX7gQgJVTVzZtQZxzrgHFDRRmVtWUeI0H7jokmq7GPJ4juKYRHfbLGv23xJjuTWBkMnk455xLrfo8cDeAJB+4c845t/dL5jmKqwmuEWyD4IE7oNYH7pxzzrUMqXzgzjnnXAuQTKCo+cDdX0nugTvnnHMtQDKBYjrwOdUfuLs5lYVyzjnXfCRz91IHgltbH4Kqxv46AF+lsmDOOeeah2SOKF4mCAyVOgAvpaY4zjnnmptkAkWGmRVX9oTdHVNXJOecc81JMoFih6TRlT2SxgA7U1ck55xzzUky1yiuA/4qaSMgoA9wbkpL5ZxzrtlIpvXYpWHDfQeEgz40s92pLZZzzrnmIpkmPK4GOpnZ+2b2PtBZ0lWpL5pzzrnmIJlrFJeZWVFlj5l9CVyWuiI555xrTpIJFOmSVNkTPkfRLnVFcs4515wkczH7BeBxSZWtx14RDnPOOdcKJBMobgQuB6aF/Yuo+2tRnXPO7aVqPfVkZhVm9oCZnW1mZwOrgf9JfdGcc841B0m9qU7SocB5wPeBj4G/pbJQzjnnmo+4gULSMILgcB7wBfA4IDM7rpHK5pxzrhlIdOrpA+B44DQzO8rM/gcor8vMJZ0s6UNJuZKmxxh/oaTPJa0IP5dGxk2VtDb8TK1Lvs455xpOolNPZwKTgVclvQDMI2jCIynhbbQzgfFAHrBU0gIzW10j6eNmdk2NafcBfgVkE7xNb3k47ZfJ5u+cc65hxD2iMLOnzGwycCDwKkGbT70l/UHShCTmPRbINbP14atU5wFnJFmuk4BFZrYlDA6LgJOTnNY551wDSuaupx1m9qiZfRcYAPwvwS2ztekPfBrpzwuH1XSWpPckPSlpYB2ndc45l2JJ3fVUKdy7nxV+GsIzwGNmtkvSFcAcgusiSZN0OcFzHmRmZpKTk1OvghQXF9d52njp61uGxlafOu/tWludW1t9weucCnUKFHWUDwyM9A8Ih1Uxs8JI7x+BuyLTjqsxbU6sTMysKnBlZ2fbuHHjYiWrVU5ODklPOyf42iN9vOHNVJ3q3EK0tjq3tvqC1zkVkmnrqb6WAkMlDZbUjuDC+IJoAkl9I72nA2vC7oXABEk9JPUAJoTDnHPONbKUHVGYWZmkawj+4NOBR8xslaRbgWVmtgC4VtLpQBmwBbgwnHaLpN8QBBuAW81sS6rK6pxzLr5UnnrCzJ4Dnqsx7JeR7puAm+JM+wjwSCrL55xzrnapPPXknHOuBfBA4ZxzLiEPFM455xLyQOGccy4hDxTOOecS8kDhnHMuIQ8UzjnnEvJA4ZxzLiEPFM455xLyQOGccy4hDxTOOecS8kDhnHMuIQ8UzjnnEvJA4ZxzLiEPFM455xLyQOGccy4hDxTOOecS8kDhnHMuIQ8UzjnnEvJA4ZxzLqGUBgpJJ0v6UFKupOkxxv9Y0mpJ70l6WdJ+kXHlklaEnwWpLKdzzrn42qRqxpLSgZnAeCAPWCppgZmtjiT7XyDbzL6SNA24Czg3HLfTzA5JVfmcc84lJ5VHFGOBXDNbb2alwDzgjGgCM3vVzL4KexcDA1JYHuecc/WQsiMKoD/waaQ/Dzg8QfpLgOcj/RmSlgFlwAwzeyrWRJIuBy4HyMzMJCcnp16FLS4urvO08dLXtwyNrT513tu1tjq3tvqC1zkVUhkokibpfCAbODYyeD8zy5e0P/CKpJVmtq7mtGY2C5gFkJ2dbePGjatXGXJyckh62jnB1x7p4w1vpupU5xaitdW5tdUXvM6pkMpTT/nAwEj/gHBYNZJOBH4OnG5muyqHm1l++L0eyAEOTWFZnXPOxZHKQLEUGCppsKR2wGSg2t1Lkg4FHiQIEpsjw3tIah929wKOBKIXwZ1zzjWSlJ16MrMySdcAC4F04BEzWyXpVmCZmS0Afgd0Bv4qCeD/zOx04CDgQUkVBMFsRo27pZxzzjWSlF6jMLPngOdqDPtlpPvEONO9CYxMZdmcc84lx5/Mds45l5AHCueccwl5oHDOOZeQBwrnnHMJeaBwzjmXkAcK55xzCXmgcM45l5AHCueccwl5oHDOOZeQBwrnnHMJeaBwzjmXkAcK55xzCXmgcM45l5AHCueccwl5oHDOOZeQBwrnnHMJeaBwzjmXkAeKGk77f9/nf566DoA31z3PB4+eCTuLWFe0jg1bNzRt4Zxzrgmk9FWoe4vZ78+mrKSMcYzjk7Q1zNq6hivKKrji9f8AYOU7c/le7kNB97BpTVlU55xrdH5EAfx++e/5r03/hZlVDdu8vaSqe+vO0q8TL7zp6+7Sr5Ka/49e/hEj5/grwJ1ze6eUBgpJJ0v6UFKupOkxxreX9Hg4/m1JgyLjbgqHfyjppFSWs1JFxdeBYseu8qruaACp5o6+bNqxidwv18JHC+PONycvp6q7ZGcRbxe8XdVfXFpcmTmUlwFQXlHOFzu/qEcNkjf/o/lsL9+e0jyccy1Dyk49SUoHZgLjgTxgqaQFZrY6kuwS4Eszy5I0GbgTOFfScGAyMALoB7wkaZiZlZNCh/z521XdZy/8t6ruozfOreoeOfhb1bufPHGP+USPHgaWpENGOJ+HR1DUJojN3XeLorZBAOpb0p7itjvZnr5n3D5gpyhNL6O9pbMzvTfdrZTP03fRYXtnhqVvJz+9HMjgvQ7bGfLlAEo7t6OszWcMsK6UlZbSJ703FWYsbPch+1cMpM9X7Uhv04l/ZbwHwMJH5tK9Ux8og50V2+ioLrTLKKdb+kGUlW7kwx1dKEvPYzlvMo5DqPhyO6/2WMuU9t8lzcrYP72cVdpKl5J9sd79KN32f+zfuScf7v6c8lIxqlMXiuhKQfFHdG+bQZdO2XxSupZeHQew4bPlfLD9bS4c8hNyi5aRlpFBp686Ud5xH0Zm7sfOnVtYt3MLUEav9gNJJ51uuzewo11f1hVv5J/Fr3BSz0vpXFZMehuxtegjhn3rQNZ98Rlj+45gK7uhoi3atYtuX63lgy6j0BfreXt1OvlFn5HRKYPBHfuw4atC2rbvRknxJg7aN4u31rzMxi/eY0TWWPr1+Dad2nfF2E1ndvJFwceU9RhG5/Ri2pbu5v+2ltGjR0faqQNFtoV9Swoppx27duyg46Ax5Bd8QLtdX5DV5wB27zLaduyBFX7Elvad6GglpHftR1pxETsy+tCjvSiyNvRKM3a370L6to20a9seddqHL7YUk/fBS3TsM5wDBx9AefsupLfJoLx0O1a6E7XJIL1tO9IFpeVl2NY8KjK6s277CkZtHUpa2S46d9wXlZewpayE7VuK6dJzX7qmVZDetgM7ykXHtu34rHgjvTtmkt6mHWkVpZSpLWlWBmlp7C4z0tPTSU8XVlFOhRmYUbx9E1269iZNbSirKEdqQ5os2PFKqwBrQxsZuysqaJPWBiREBZKoqDCUlkZFRRkind07Pqdtp14gkNLADMMwMySRpjTKKoy0tDQqystJTxMVZaWkt20PQLlV8FXZV1RUlAMKhlVUYIR5W5CfWQUgDKOiopy0tHTS09KDncLwY0CFQXqaqvKvMEhLE+LrHUgLczIgTaIiTGsYMigtL6WN0iAtnTSlYRUVQRkqykHVf/Pl5RVBfhgKFgKSwIwKjLSwTlUktpYUU16R0r9GFHdv+ZvOWPoOcIuZnRT23wRgZr+NpFkYpnlLUhvgM2BfYHo0bTRdojyzs7Nt2bJldSqnmTFq7qg6TeOcc83Nkh8so0MYMOtC0nIzy06UJpUXs/sDn0b684DD46UxszJJW4Ge4fDFNabtHysTSZcDlwNkZmaSk5NTp0KWV6QmUDrnXGN6+/U39jhCaSh7/V1PZjYLmAXBEcW4cePqPI+VrCQnJ4f6TLs38zq3fK2tvuB1ToVUXszOBwZG+geEw2KmCU89dQMKk5zWOedcI0hloFgKDJU0WFI7govTC2qkWQBMDbvPBl6x4KLJAmByeFfUYGAosCSFZXXOORdHyk49hdccrgEWAunAI2a2StKtwDIzWwA8DPxJUi6whSCYEKZ7AlgNlAFXp/qOJ+ecc7Gl9BqFmT0HPFdj2C8j3SXAOXGmvR24PZXlc845Vzt/Mts551xCHiicc84l5IHCOedcQh4onHPOJZSyJjyagqTPgU/qOXkvILUt8TU/XueWr7XVF7zOdbWfme2bKEGLChTfhKRltbV30tJ4nVu+1lZf8Dqngp96cs45l5AHCueccwl5oPjarKYuQBPwOrd8ra2+4HVucH6NwjnnXEJ+ROGccy4hDxTOOecSavWBQtLJkj6UlCtpelOXpyFJ2iBppaQVkpaFw/aRtEjS2vC7Rzhckv47XA7vSRrdtKVPjqRHJG2W9H5kWJ3rKGlqmH6tpKmx8mou4tT5Fkn54bpeIWliZNxNYZ0/lHRSZPhese1LGijpVUmrJa2S9O/h8Ba7nhPUuWnWs5m12g9B8+frgP2BdsC7wPCmLlcD1m8D0KvGsLuA6WH3dODOsHsi8DzBe+KPAN5u6vInWcdjgNHA+/WtI7APsD787hF292jqutWxzrcAN8RIOzzcrtsDg8PtPX1v2vaBvsDosLsL8FFYrxa7nhPUuUnWc2s/ohgL5JrZejMrBeYBZzRxmVLtDGBO2D0H+F5k+FwLLAa6S+rbFAWsCzN7jeBdJlF1reNJwCIz22JmXwKLgJNTX/r6iVPneM4A5pnZLjP7GMgl2O73mm3fzArM7J2wezuwBuhPC17PCeocT0rXc2sPFP2BTyP9eSReGXsbA16UtFzS5eGwTDMrCLs/AzLD7pa0LOpax5ZS92vCUy2PVJ6GoYXVWdIg4FDgbVrJeq5RZ2iC9dzaA0VLd5SZjQZOAa6WdEx0pAXHrC36/ujWUMfQH4AhwCFAAfD7pi1Ow5PUGZgPXGdm26LjWup6jlHnJlnPrT1Q5AMDI/0DwmEtgpnlh9+bgb8THIZuqjylFH5vDpO3pGVR1zru9XU3s01mVm5mFcBDBOsaWkidJbUl+MP8i5n9LRzcotdzrDo31Xpu7YFiKTBU0mBJ7Qje2b2gicvUICR1ktSlshuYALxPUL/Kuz2mAk+H3QuAC8I7Ro4AtkYO6/c2da3jQmCCpB7hofyEcNheo8b1pEkE6xqCOk+W1F7SYGAosIS9aNuXJOBhYI2Z3RMZ1WLXc7w6N9l6buqr+039IbhD4iOCOwN+3tTlacB67U9wh8O7wKrKugE9gZeBtcBLwD7hcAEzw+WwEshu6jokWc/HCA7BdxOcf72kPnUELia4AJgLXNTU9apHnf8U1um98I+gbyT9z8M6fwicEhm+V2z7wFEEp5XeA1aEn4kteT0nqHOTrGdvwsM551xCrf3Uk3POuVp4oHDOOZeQBwrnnHMJeaBwzjmXkAcK55xzCXmgcC2OpPJI65orwiYQWgRJh0p6OOy+UNJ9NcbnSMpOMP08SUNTXU7XsrRp6gI4lwI7zeyQWCPCB5lkwZOte6OfAbd9g+n/APwHcFnDFMe1Bn5E4Vo8SYPC9vjnEjzJOlDSTyUtDRtX+3Uk7c8lfSTpdUmPSbohHF61py6pl6QNYXe6pN9F5nVFOHxcOM2Tkj6Q9JcwSCHpMElvSnpX0hJJXSS9JumQSDlel/TtGvXoAowys3eTqPPpkSOqDyV9HI76F3CiJN9JdEnzjcW1RB0krQi7PwauJ2jSYKqZLZY0IewfS/AU74KwwcQdBE0cHELw23gHWF5LXpcQNBFxmKT2wBuSXgzHHQqMADYCbwBHSloCPA6ca2ZLJXUFdhI013AhcJ2kYUBGjICQzddNNlQ6V9JRkf4sADNbQNhUg6QngH+Gwysk5QLfTqJuzgEeKFzLVO3UU3iN4hML3k0AQRs/E4D/Dfs7EwSOLsDfzeyrcLpk2sSZAIySdHbY3y2cVymwxMzywnmtAAYBW4ECM1sKYGErqJL+CvxC0k8JmpmYHSOvvsDnNYY9bmbXROqaEx0p6T8IlsfMyODNQD88ULgkeaBwrcWOSLeA35rZg9EEkq5LMH0ZX5+qzagxrx+ZWbXG5SSNA3ZFBpWT4PdmZl9JWkTwUpnvA2NiJNtZI++EJJ0InEPwRryojHBeziXFr1G41mghcLGCtv6R1F9Sb+A14HuSOoTXA74bmWYDX/95n11jXtMUNAmNpGFha73xfAj0lXRYmL5L5HrBH4H/BpZa8Aa2mtYQnlqqjaT9CBrGO8fMagaFYex5Csu5uPyIwrU6ZvaipIOAt8Lry8XA+Wb2jqTHCVrc3UzQRHOlu4EnFLwp8B+R4X8kOKX0Tnix+nO+fiVnrLxLJZ0L/I+kDgR79icCxWa2XNI24P/FmfYDSd0kdbHg9ZiJXEjQuupTYR03mtlESZkEp6I+q2V656p467HOxSHpFoI/8LsbKb9+QA5wYLzbdyVdD2w3sz/WM4/rgW1m9nC9C+paHT/15FwzIOkCgnci/7yWZzz+QPVrH3VVBMz5BtO7VsiPKJxzziXkRxTOOecS8kDhnHMuIQ8UzjnnEvJA4ZxzLiEPFM455xL6/5osO8ZjIYM1AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots() #create an empty plot\n", "\n", "df_fft.plot(ax=ax) #use the dataframe to add plot data, tell it to add to the already created axes\n", "\n", "ax.set(xlabel='Frequency (Hz)',\n", " ylabel='Acceleration (g)',\n", " title=filename)\n", "ax.grid() #turn on gridlines\n", "\n", "fig.savefig('fft.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "883JemLPjSk9", "pycharm": { "name": "#%% md\n" } }, "source": [ "### FFT from PSD\n", "Here we can use the output of a PSD and convert it to a typical DFT. This has the benefit of allowing you to explicitly define the frequency bin width." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 447 }, "id": "s1EOA2jaOnMF", "outputId": "0cef351f-f361-43ac-92d3-2d285a43bf73", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\"X (2000g)\"\"Y (2000g)\"\"Z (2000g)\"
Time
0.004394-0.122072-0.122072-0.061036
0.004594-0.0610360.488289-0.366217
0.0047940.1831080.122072-0.061036
0.0049940.122072-0.122072-0.122072
0.0051940.1220720.122072-0.244144
............
27.691064-0.427253-0.152590-0.671397
27.691264-0.122072-0.335698-0.305180
27.691464-0.183108-0.152590-0.122072
27.691664-0.3051800.030518-0.244144
27.691864-0.305180-0.030518-0.366217
\n", "

138440 rows × 3 columns

\n", "
" ], "text/plain": [ " \"X (2000g)\" \"Y (2000g)\" \"Z (2000g)\"\n", "Time \n", "0.004394 -0.122072 -0.122072 -0.061036\n", "0.004594 -0.061036 0.488289 -0.366217\n", "0.004794 0.183108 0.122072 -0.061036\n", "0.004994 0.122072 -0.122072 -0.122072\n", "0.005194 0.122072 0.122072 -0.244144\n", "... ... ... ...\n", "27.691064 -0.427253 -0.152590 -0.671397\n", "27.691264 -0.122072 -0.335698 -0.305180\n", "27.691464 -0.183108 -0.152590 -0.122072\n", "27.691664 -0.305180 0.030518 -0.244144\n", "27.691864 -0.305180 -0.030518 -0.366217\n", "\n", "[138440 rows x 3 columns]" ] }, "execution_count": 117, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "RD-ya2m7jUfb", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def get_fft_from_psd(df,bin_width):\n", " fs = len(df)/(df.index[-1]-df.index[0])\n", " f, psd = signal.welch(df.to_numpy(), \n", " fs=fs, \n", " nperseg=fs/bin_width,\n", " window='hanning',\n", " axis=0,\n", " scaling = 'spectrum'\n", " )\n", "\n", " df_psd = pd.DataFrame(psd**0.5,columns=df.columns)\n", " df_psd.columns\n", " df_psd['Frequency (Hz)'] = f\n", " return df_psd.set_index('Frequency (Hz)')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "axwfbM7T1Gvs", "outputId": "ea167946-0426-49e5-8e02-228616f8c3ea", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df_fft_from_psd = get_fft_from_psd(df,.25)\n", "\n", "fig = xp.line(df_fft_from_psd)\n", "fig.update_layout(\n", " title=\"FFT from PSD\",\n", " xaxis_title=\"Frequency (Hz)\",\n", " yaxis_title=\"Acceleration (g)\",\n", " #xaxis_type=\"log\",\n", " #yaxis_type=\"log\"\n", ")\n", "fig.show()\n", "fig.write_html('fft_from_psd.html',full_html=False,include_plotlyjs='cdn')" ] }, { "cell_type": "markdown", "metadata": { "id": "yBbcLxtXlrO6", "pycharm": { "name": "#%% md\n" } }, "source": [ "## Moving Peak Frequency\n", "Now we'll introduce a for loop to go through all the columns and find the moving peak frequency." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "kLcuacqMnlQA", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "def moving_frequency(df,steps):\n", " #df.index = np.linspace(df.index[0],df.index[-1],num=df.shape[0]) #resample to have uniform sampling\n", " d_t = (df.index[-1]-df.index[0])/(len(df.index)-1)\n", " fs = 1/d_t\n", " length = df.shape[0] #full length of time series data\n", " nperseg = length/5 #average 5 FFTs per segment\n", " n = int(length/steps) #the number of data points per step, we will not compute\n", " df_peak_freqs = pd.DataFrame(columns=df.columns)\n", " for i in range(0,length-n+1,n):\n", " df_sub = df.iloc[i:i+n]\n", " df_psd = get_psd(df_sub, nperseg/fs)\n", " time = (df_sub.index[-1]-df_sub.index[0])/2+df_sub.index[0]\n", " df_peak_freqs.loc[time] = df_psd.idxmax()\n", " return df_peak_freqs\n", " \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 542 }, "id": "ZDARGyY_pGc3", "outputId": "d516952c-3d85-4061-dc3f-b11e2e901466", "pycharm": { "name": "#%%\n" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "
\n", "\n", "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "df_peak_freqs = moving_frequency(df,40) #get rolling peak frequency for 20 different steps\n", "\n", "fig = xp.line(df_peak_freqs)\n", "fig.update_layout(\n", " title=\"Rolling Peak Frequency\",\n", " xaxis_title=\"Time (s)\",\n", " yaxis_title=\"Peak Frequency (Hz)\",\n", ")\n", "fig.show()\n", "fig.write_html('rolling_peak_frequency.html',full_html=False,include_plotlyjs='cdn')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "FkRATc7ImdFo", "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [] } ], "metadata": { "colab": { "name": "Webinar-Intro-Python-Acceleration-CSV-Analysis.ipynb", "provenance": [], "toc_visible": true }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 4 }