SSSI_notebook.ipynb

{
 "cells": [
  {
   "attachments": {
    "Steps.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0oAAAB3CAIAAACytCdUAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAGInSURBVHhe7b2FdxXJ+rb9/UHv75w5M4PrwODu7u7u7jpogOBuwQlBQ0IEQjzE3SHEE+IBhu/q/TRNsyPEIFbX6rVXVXXt6uqSp+5qqf7/vigUCoVCoVAomhBK3ikUCoVCoVA0KZS8UygUCoVCoWhSKHmnUCgUCoVC0aRQ8k6hUCgUCoWiSaHknUKhUCgUCkWTQsk7hUKhUCgUiiaFkncKhUKhUCgUTQol7xQKhUKhUCiaFEreKRQKhUKhUDQplLxTKBQKhUKhaFIoeadQKBQKhULRpFDyTqFQKBQKhaJJoeSdQqFQKBQKRZNCyTuFQqFQKBSKJoWSdwqFQqFQKBRNCiXvFAqFQqFQKJoUSt4pFAqFohFQWFh4/MTpZ8+cdH81SUpKfvjoid2N28XFxXpQpYSFhYeHR+L4/Pnf9PQMCVQoGgtK3ikUCoWiofPvv/9u3rrj5KlzU6fPefXqtR5aZRITk9q265KS8v5/f7T19vHTQyumpKR067ZdpaWlt27fGzx09Nx5i/UdCkUjQck7hUKhUDR0Pn786Ob2EsenT58O29hKYNW5dOnqocNHcbi/9MjNzZXASkBNnj5zvkevAYOGjHR4+Pjz58/6DoWikaDknUKhUCgaE2fOXvz8+V/dYyE+PgFBpnvKw5B3kJmZuXL1+oKCAvGWCyLy6LGT4RHazVmDvLy83NwPukehaNgoeadQKBSKRsPnz59nz1mI/NL9litt+w/YlJaW6v4y3Lx1Z8CgEb37Dp4zd9HjJ89QaecvXM7KytJ3l8fTp44uru665ythYRGeXj66R6Fo2Ch5p1AoFIpGA1rtP7+1SktL0/0W4uLiq3717ocUFRXPmrOwrF5EF+bk/PjGrkLREFDyTqFQKBSNg+zsnE5/9Rg3furNW3f1IMv1vGnT51T+PqxZ3iUlvz1sY1vJy7CIRXf3V7rHhK+v/737DrpHoWjYKHmnUCgUisbB3bv2v/3eJjAoZPrMuXqQRd6NnzitqMryDv7q0isqKlr3lCEwKPjjx4+6xwTyziwrFYqGjJJ3CoVCoWgEILmGDh+7Y+de3KPGTCosLJTwGsi7I0ePVyLvAgICy31VVsk7RSNCyTuFQqFQNALu3XvQ5e/eclP17NkLiC0czs4uyD5D3uXnFxQU6LLPTEXyjv9mZWVLoEFwcOgDh0e6x4SSd4pGhJJ3CoVCoWjolJZ+7Nt/6OYtO+S6WmZm1sLFK+7cvb//gM2nT59E3hUVFQ8ZOnrXrn3yFzPIuylTZ+sek7xbt37z//tPCwk0yM7Jadeh6/yFSzdt2b5z1z9sl69cRwgqeadoRCh5p1AoFIqGTkFh4b79h5PfvtP9X76cu3B56LAxoaFhuH18/RB5Iu/u3LWXCAYJiYmOz51lw1tcXOzk7BIaGo4beTd5ykxLrO9ISkpesmxV/4HD2fr1H/pnyw7e3r4FBQUkpcdQKBo2St4pFAqFoikQFBTSvWf/wsIi3f8jsrOzW7bqWJVvlPn5Bzg8fKJ7FIrGgJJ3CoVCoWgK3Ld3iIyM0j3lERYeobssxMUn3Lv/QPeUx7///isvcDxzdH4TGCSBCkWjQMk7hUKhUDR9Ll66On/BUn51fxX48CFv3ISp48ZPtT1xutx3aRWKBouSdwqFQqFo4qDqjhw7gYPfaik8haKRouSdQqFQKJoyhrYTlMJTNAeUvFMoFApFk8VK2wlK4SmaPEreKRQKhaJpUq62E5TCUzRtlLxTKBQKRRPkw4cPkybPXLN2ExsOPfTLl/kLlhqBxNFDFYqmhZJ3CoVCoWjiIOZ01/duhaKpouSdQqFQKJo4St4pmhtK3ikUCoWiiaPknaK5oeSdQqFQKJo4St4pmhtK3ikUCoWiiaPknaK5oeSdQqFQKJo4St4pmhtK3ikUCoWiiaPknaK5oeSdQqFQKJo4St4pmhtNU955eftOmjxz/oKlamvC285d/+j1XSlh4REXL11VW/Pc9EagaPYoeadobjRZeacse5MHhae7KgUV+NzpBU1Cbc1tYxRH3OvtoFKMbxuorSlt3Xv2Nz5KgVccYHarqm+GWxUvDTR2lLxTNFaqLu/evn2nexTNiSPHTlRF3qEAmom5b24wkFcu79hrDlc0E5pJpSt5p2isKHmnqBwl75o5St4pykXJu0aMknfNgV8p7+Li4n18/Ky2vLw8fXej4t9//w0LCyf/yclvCwsLf81ZJCYmJSQk6p5fRT3KOwp53vwl/+8/Ldju2z/UQy0cOHiEwI6du2/YuFUPqpjDNrZ/tGg/ZNiYktJSPagCzEfEwTZx0vQ/WrS7eeuOHuN7Pn/+nJWVLe5Pnz6PGTd55ap1JSUlElIDtASzjQQ/jR47aeXq9bVJsPbUu7wz2w3cemjF5OTkSuSioiI9qKZQ8mFhEZJaZmYWKVeUZkhIGHGMxlAVsBuScumPmqVAZvz8AgoLa3tSdYWSd40YJe+aA79S3mVnZ6/bsIWB8+GjJwkolcSkHj0HvPb01nc3HjDHvfoMmjFrPqeAdOjbb4iXt4++76fx8eOn1m3/at22M4JSD/ol1PvVuz9bdkDljJ8wTfdbaNm6Y7sOXafPmKv7KyUnN3f23EVVkXeCyDtxl5Z+3PvPwYrkHcP5nLmLxI00lD8iAiSkBpDgvAVLxG0kmJtb8wRrT73LO+yGzZHjlMOOnXvT0zP00Io5d/7S8BHjiB8TE6sH1QhmU9179key083ZFi5a/lfXXkzq9N0m/v33y4yZ8zji48fP9KAqEB0TO2DgcP5V0UnRAHx9/XXPly9u7q/+778tjxw9ofvrGyXvGjFK3jUHfvHN2T17D2DODJsVEhq2a9c+cTci4uIT/vu/1t7evuK9cPHKL5B3sGz5miVLV/3i6531Lu9at+lMK/2jRbuoqBgJOXPuoqeXd4dO3aoo72DBwmU1k3eMqYyyZ89eFK8VmzZvnzhphu758iU8PKKgoED31Ij1G7dOnjpL91jeWK9lgrWn3uUdXL5ynRqhKer+H/H4yTPi10belX78OGToaGZxKSnvJQQD2H/g8HLlHSQlJXPEask7QInyr4rkHT29a7c+usciN7v16Ofi6q776xsl7xoxSt41B+pR3kVERPK7d99By57GBCaes3jg8IiBH+/TZ8+9fXSp9/PgWMUW5KC/jHqXd23adYmNjaPADT29/8DhzKwsK3lnb//wjxbt2ZycXCTk06dPCQmJf7bs0L1H/2nT5xjyjkCJuWPX3nLFE8diw5FfULBl604JhHcp7xlf+SNiMSw8PDMzi2gjR0+Ii4vnWO7ur+SPjNZ379rLIdLS0hmhBw8dVVr68fad+9pBd+yV1Pz83/ToNYCQq1dvSIgkOHrsJEnw+XNnSTAjI5O9np7ekub585fxEkG8N27cXrx0JY6qq59q0bjknfQRZBbxayPv/APekIK59mH/AZtfKe8iI6PYa/R3ObWPHz+Kt95R8q4Ro+Rdc6C+5N3nz//+/mc7CczMzGzfsVvL1h0PHjo6dtyUPn0HFxUVMcItWbpqx869c+YtXrhouTx+ZHfjNqOpzRHb1Ws3Dhk2OiAgcO8/B//3R9tWbTq9fu3VuUtPEre7cYuYPj5+i5estDlyfPLUWes3bi0tLZ01ewF7t23fvW//YUbc+QuXMkAS08vLZ9GSlUeOHp8ydfabN0Fp6ekzZ2kxHz56yt4NG7f6+PrhMPPuXQoR2JYuX716zcYXL9wkPPfDh63bdq1dt3nylFmHbWwtj3y1Gzh4pK3tKTQKcpA4d+7aczqME8NHjrt2/Qa2Ozn57YqV6/gjdcHJYsGvXb/52+9t+vYfivWnQHr1GUQgc/f//q81B42PTyADLVt3wu38wnXVmg0oGF8/TTF7eftMmjxz7bpNSJlFi1d0+qtHUXGxJWs1pyHIO345WQpTQk6cPGsl7w7ZHBs5akJ2djZyivZw8tRZAmkJUlzU/qjRE0XeMWR26do7OiaWmEgiil1SMMO/2Gg/HTt3Nw/wc+ctio6ORfNRp6GhDPQRRBswcPiTp45UEIP0qtXrCZHRevyEabjv2z+MiorBMXLU+NCw8OCQUNwurm6pqWk4zpy9QDa69+w/YeI0UpAESVwSJA4NgxDknf2DR337DaUNFBYWDRw0Yuu23Z8+f3746Al7aTYcYsas+T17D6QrSVbrkIYm74zG/9zpxcrV61u06hgUHCLRVq5aTy3fuXt/wsTpRBB5FxgYTFXuP2gzbcZcOgvyaNfuffQv+iCzMgr/7LmLBQWFkoIBTYgUdM9X+G9u7gd+qRcSWbd+y6zZC5ndscss75DyGJnlK9bOmbeITsGsgDbM3rS0NHL4f/9tOWnKTEt63+QdE5gx4ybTsHv0HLB+w9bCwsK3X40MjSo6RrtujR3Au3vPftzMDejpxBdTQ4hRLE7OLqtWb2jRqgPmUTvGT0PJu0aMknfNgXqRd/0HDkPD4dBDv3xxd/fAO2rMxEuXr/XtP+RDXt6YsZN79x3MLsY2LPjDh09CQ8P/81vL/Px8AkWrif2aMGk68g4HOSSQQT0vL/+vrr02bNxG4KPHT//zWyv+i/Rh74qVawnE8qIFGSyJyRDOYEkgBvGBg/b8PiMuQmrN2o24T5+5wG9Znjk6McqSIBvp5+drF4GQAp3+6okDDTp67KSSEk1Tojbi4hNmzpp/zPZkREQk44pchULJ9eg1gDNFgP7drS8hHz9+atWm83W7WwUFBQMGDWcvOZFrQs8sQ4iLixtu9Aruf/YdxC1DC4616zfjWLpsNeYeBwKawPBw7fpoLWkg8k6kTEJC4ps3gfftHczyDonMLiYDluj6U1CfP3/u02/wuK9P7Bk3Zzds2ErtREXHsIkaK3s1lEA2HFSryDsp86HDx86ZuygwKDgyMlqLZ4lpvjl79uxFQkTerVylKTMJx3Hw0BEczFtwn7+gXX5DcTKdIBvU9YBBI6Rhs9d8c/b0mfOE0AW69eiP0JfAHTv2Eogi8fT0RiugZgik6mnzScnJEqcOaWjyDpBWeJF3uHFs3KT1dPeXr4wZo3FztqiouHefQfPma48zenr50PuY+OGmYQweOhqFXtFTleXKO4GZJ8WOOsSNuKddUa1meXf6zLn2Hf/GQZvEiNGKmMixF3lHIJksK++OHjtBmh6vPQkkhFmrONgsETVoIa3b/oW8I0Emw/YP9JeNOCl6Bw6xro6OTrhxrN+wxbL/Z6HkXSNGybvmQL3IOy8vn+zsbHSPHvpV3ukey40nvH37DWE0nTFrPqJn1659N2/dRTzJm2vR0doVkYrkXVp6Oo5hlsF4/MRpyC/H587IO0Tby5cexDTkncR88tSRwIiIKDQBDrnYI3JB7iCXCzkJDgklZVKQgYcxeJ1FZr1/n8p0HAfybvjIcVpsC0T+/Y+2KDAyNmjIKP6YmqZdxWFeTsjsuYvatu/COC3y7rDNMf5SibwbNWaSlqjFlIu8W7FqHeeLQ+RdVJSuQmpDA5F3KN2u3fogztq0/Ytqqoq8ozrKlXeDBo987uTCRnU8c3SWCGb4O5u4g4JCSH/33gOIwOTkt5wju2hR8gon7hrLO/S9q9vLwMDgfgOG1Vje0YmCgoIJbG7ybux4vTQIF3l3+Igt/UgCDXn34UMejsFDRtG/EFUUF2ZE4tAe6C8VPclarrxDLMqLTShFUibNlq07lZV34yZMnb9Qs6tUnNzM/aG8Q7GR5pvAINIkpHJ5J1bLkHe0H2ofB9aVWbEEEkHJuzpByTtFY6W+bs7i3rR5uwRCufLuwEEbRnHZGMyqK+9u3Lxj/L2kpBR5xzSX0ZSY5co7M9t37iG8y9/fnms2Y/X8TZe/e8sLIoa8Mygr7zp06h4YFGxkTM6U4jVCMOIi727cvM1fKpF3c+cv1hK1mHKRd8HBoUOGjT5sY4v0mTJttuytJQ1E3n38+HHw0NG4/fwC8FJQZnnX9e8+hrwbMnR09579CZw2Y25ZeXfM9iQRJJC2JE3RCsqTTfdY0pdreGvWbiy2PPx04uQZUQlEE3kn79ZUXd5RleSfOQAJDhw8sqy8kwQNeTdp8kxD3q1YuY7/kpQu7yy3JpubvFu0eLm4CRd5d+jwsYrk3YWLV4z+VVysrzLDFO7//tvylYd2wawsMZbHPdes++68iBwdE5uamor1kD+We/XOkHcGP5R3N2/eodHm5eltwEre+VrafFXknWHPidBU5R398V3KezYc8tAOowPezMwsiWAg0YgD4mb7aHkmp+ooeac9wITFoUn92bLDX116sTHmRURGUQF6jEYF2cbIHrM99flzzfNPInHxCfQ9SuPosRMvX3qItmAUnDt/yceP1WtkZfH01C5Htevwd23eratHeWfGSt59/vx5wsRpQ4eNES+G28PDE1XEYIb4IaRyeYcdpPlt37FH+/OXL5OnzIyLSyhX3hETB7kiEO/Zc/oLkpgDkpK7P2VB3g0bPlb3aKqu3317BxyUUo9eAyTw4sWrGH0reUe2Me4iUGDe/MUYphkz5/XrP1RCGPg5aI3l3YaN2lM7jEzJyW/rquvVr7zjLJB0ci7y5gEnaHFnYnAQcLLrxQu31m07JyQkUr90CswRgUFBIcT38PDKycn9z2+tEL5UK/H7DRi6WeTauk1lH6wkAv9iwwHUNeOuzRHt8abuPfqfOKk91bd4ycqQ0DAcg4aMGjNu8q3bd6Ojtce8zpy9wB/T09NxVy7v3rwJ+rNFexpwUFAwNrNtuy5oOCIMHDRy/IRpN2/dkQRPndbkHcN/bGw80TgvonX6q8fdu/bsVfIOCBd55+PjW1beoZ779h+6bLmmfoA5QLilMVMFd+7aL1qyYvachbLLCqp+7z8H/mjRPtx0/X7wkFEk6PDwcddufd6+eydXlMvKu/MXrmB/5C83btym+f1Q3nXvOWCtRUqKJbSSd7K4oyHvMBpYPEPeaXcnLDdkm4O8o/tv37FX2jm9YM+eAzJGM1gPHDTCEkXH1dUdgy8PWAPGQQqzulQm70pLS5OS3pbdjKM2WKp79U6eCzaa1NWrdnjDI348MDRA5JoK/bA2i1CkpLzH+D59qg3MkZHR02fMZezBarRq05nEExKSJFqNycrKXrZ8TSOSdwgXtA7nzlCEkfr49RUwzNaixSsIx6gZyi8uLh7Lu2jJSgKN+ykbNm7DmGLrZ86aT3yRd3Y3blPORFu3fgv2nSE/OzuHYRuBvmXrzn/2HXr67Dl6cemyVdjBFSvXOTu79B8wrEWrDustRvO+/cM58xZzggcPHWVybzmOxtTpc8iY7vkehvx16zcjHzmLCROnL1uxBoNLOGZ6+cq1HIKDMsA4Obmg9hiMyVturj5A3r5zf8SoCYTs239YBubkt28XLFw2b/4SAq9dv0nIhYtX2nf8mwbz8NETKRn2vn+fSmXhXrtuM4UzctQETtbmyHG5wDNw8Eg/P//rdrc4nLG5uLrT3iyHrTn1KO/IPDJINrxUYmKS1msYVo3w6K9vR9LdJCQ1VRtEBVQgIYFBIQzSRmT6i8TMKPMigvmI5k1eWQgNC6fq8b57lyLxGbnxMnJ8+vT5w4c8I352drY4IiKiqGVxc3RxsBUVF9OhcCQlv5X8SN7MCRqJsLELYSduuT9ozqokJZuWrTql3uUddmPylFk0cpQ086jEpGTMSJ9+Q44eOym3UAcPHcXcj5jbduxGOdErR4waT/iatRtpMwjoadPn0GH37juIMiMacpCB39XtJXMnJg8LFy03txkDxm70+sZN20mQzebocenF71NTJ06esWr1BgwRbf6PFu0O29hifDji7LmLmCUiARkHkZLYAWwXecjJzR06fMzmLTu2bd+9ZdvOv7v3tbd/SMyJk7RXQDZt3v7A4XHP3gOPHD3B/GHKtNkdOnaLjIrmxNFzO3btlad1OQqGjnN//Zq5qh9nisUYP3EaR2EvxcJZ9+47mETEJgwaMlKK5SdRL/IuKjpG2r+A0RZ55+jo/He3vuZVY54/d0bw/Vx5FxERuXrtxmO2Jw/bHENxHzl6HDcmm6p9l5LCvL/sFcWfBwaC1pxqmUP8kFrKO+wpXmNgbnS4uLjJdYIaQxf8/Y+2xo28Gzdvy6Uj5nnyeG/twbg0InnHgERnM9BDLehBnz6Z5QhmUQJ1vwnz1TsjWSO+JPLVqz1OBxa3hhENZJeE8CteCWe0Fm+5cAirwwmSGUnKyBjIXkFCjMOBEVOSMrzGIcCcmjnczFDTNcWIyKievQY26jdnSRMrpLZqbYYUqyvqXd6ZWz5YeQ3MMa36puG1pPetA1pFK5dy40ggv2UzI9G+hn+7+2P8RRxg/i8RxEGgEcEcaLgFI0H4/ojWsOsnUS/yDlvExFuefwCmT3L6iHUfX/89/2j3YYAZfnR07E+Xd8HBoeIoLCr6q2svJgTijY2N+/jxI4cs+YUfnOFU4+MTmFvo/krBUtRe3t29Zx+fkDhx0oy58xYHBYVMmjzz1m1d8D16/IzwhYtXME2REIPjJ04z31q8dOX2HXtwkAiBTLCIP3vOQuTR5KmzmFdJvdo/eEQ4W2jotyeicEugLErJiTPRwUuCoqcpeZIi5JDNMZmL01jxbtqynRnPq1eve/cZ1H/gcCzXAwc9ffLPJFLugsGHvDzmjmQP4f7M0QmHhBvQ2iiBrt36ZGVlk4G8vHzqmqOTMps8mo3UpsT4r/2Dh0j/a9dvMv+QwxFh5uwFy1esKSrS6isvP//osROEz52/hGm95QjfyTtZYZWzoOJkb1X4lfKO5lRX2/4DNpQtubIKr/22b/+h//3ehnayeMlKq10Nf8PCrFq9XtzrN27t3KWnPHdfG+pX3hmnprYqboYUqyswKUaa5hHdcLPXHF63qDZQm61aY0F1+XmVXgmopt9+b9O+498MrK9fe4kGAOQd4qrr16el8wsKCgsLf7q8S0rSLyRaybusrCzUd1FRkWUC8C9ik7EfBYCDDEFOTq48aCmgydhV7vUkwtmMyCSLl79L+h8+5PFLCEfhWJYj6pcN0A2WNLVn1ckYbiSLkU7N5N3osZMuXb7GxskuWboqLT0jLS19/sJlf7ZsP3PW/Hnzl/z3f6051o5de4cMHU3eTp0+93//bWl1bZwqHD5yXOu2nUPDwt++S+nyd++t23YhjJYsW/1Hi3aDh45et37zf35rxbnYHj89avTE/Px8NNmfLTvICmTIwVZtOiFkZ81eMHHSdAzEBMsV7KLi4slTZvYfOIw/IjSTk98yCej0V4+Y2Lhjtqfu33eg2BcsXObm/ioqOmbM2Mm//9GWciTBfgOGcWq379zz9NIW1/D186cM5WFnS1Iz2rbvsmzFGkvev0Fhjp+orX0l29LlWgSyunzlWryRkVHUJufOqSFYW7XuNH7CVA5BnCHDxhDh+vWbcrgXL1ypu7XrNv3+Z7vc3Ny16zYTKA9gGfKOP8rreBs2bSv7okAl/Ep5R4uqkw31v3jxStk8PDyt9tZy83jtuWz56nPnL1uFN4rt5SsPm6PHpWQOHDzi8PBR7Wutfp+9U9Q79SvvwKqRq63qW1V6bo2pF3kHaKFDh48xw/+7e9892ivt+tU7fhnl7957gEOWwvnp8s7ASt6Bn59/2/Zdw8MjyN+OnXvPnb+UkJgUERF54NCRGzdvI7YYqiXrHh5ezMXz8wtOnTonT/kYOLu4+vkHsGvV6g3FxSVotWnT5+JFDD1/7lxa+nHegiVnz128ffsewi0iIgqJEBDwhj9SHByUmDZHjpOm7YnTaKyYmNjZc/VHTWkcNZB3y1eszcjIZCP/+o4v2oLy7OLcyQ91g6pD1thZlp+Ni4vv0Kkb2ZOYgsi7nr0HihehgyhMz8h4/vwF6aCDiUA6QDq7d2vLPMKo0ZOIhgPlt22HtpLZyVNnjx476e3thyiUy6gODo/5S2JiEjOALdt2IbBIhEIm21OnzcGL5BX5e/DQUUPeLVy0nONKOA7kBY1m8JBRsvQRv5279DQ/tmVAtLCwiD9atKfY+WPAG+1moixagbzLzMwkw/JcBSnI20+AkjMfTh6tQMgyZSGr79+nEigPaRny7siR46PGTET8meV7VfiV8k7RGFHyrplT7/JO0TCpl0q3el6WsV4EnMg7FAuygQFRPtNcn/IO5sxdFGp5Nuvhw8fy7r3IHfQHbmSZzRFb/tL17z4WbVSK7EOyaP/8yu49++Msy7Lfv+/AL4oNWYYjL097JZAIy5avlis9woxZ85CDROj0V09+4dLla6Q5bsJUVBdlse/AYYlZM3ln3Jw1I/JO92jmQHtZff7CZWjKY8dPHzl6IiQ0bOXq9bLJbWuzvLt7z574aenpIu8kEIqKi/Ea8m7M2Mmyl1/zSGN5ZbL1lq27tMPZnrI5ejwnJ9fV1V1WGELnZWRkUPfytXuaS0yMtmKZlbwzEiQO8g4JtW//4UlTZqKPJ02eQVbRVRLBINZy+1XglFu27iR1Z8g7/jJn7uLtO/bExye2at1pt+X9TUDeyVtUQEyRd6hAZiqvPDxPnDxDoJW8S3mfunXbLsQi+a/WjE3JO0XlKHnXzFHyTlEu9VLp8QmJusvCzNkL5HKGyDv0DGO98SmaepZ3R44eN+TdLcsVLJF3lp26vEPEDB0+pqCgULbS0o+Oz52HDh/LhqO4uLhj5x7tO/6N5vj333+nTpsdEhJmRCYRDmG+kIYXeUeaI0aNN6JJmh07d+fQyV9H8Z8n7/Ly8n7/s+3Fi1fEi25A6KDeZKM+rOSd3Y1bqJb0dP3qnQRCcXHJHy3aG/JuyNDR7dp3xdG6bWdj+bHMzCw/vwDtmqXlyhkkJSVziI2bttEsSkpK9u07xLnv3XsANwV45859WcPMSt7t3qN/+Z4MUNQ4kKfnzl9q37HbocNHy71ghqyMjf2m8Dp36SXlacg7DteiZQdUJon4+fnjlZhoOyP/xETekeEWrTrK8BkTE0uglby7fOW6PKL39Nlz4ypgVVDyTlE5St41c5S8U5RLvVR6dHTM+o1bE5OSEhITY+Pir9tp9wDz8vN37tonUgodZf/gYUlJaWJiUp++g4mDRkLkeLz2ZNwMC4tAXVhSqio/V97xl3kLtBvJAiJUd1l45fEabVFYWIhi5TT2HTgsy+VDUJC27EK58q5smo+fPOO0MzIyjSV5fp68Q8fs3ru/d59B71JSOPSChcusnikUedeqTaeQkNDk5LeU287d+5A4VvKOdE6eOjds+Nj8/AIfX3+k3mtPb8KvXruBniNlhB3ni3YcO37KwsUr2PXkqePmrTtFLTk81D7ksmjxitDQcCSRXCS7cPGKfOClcnn3+fO/U6fP2bV7n2yyGpYVyLs27brIFTu7G7d/+70NyhK36epd8R8t2hmJeLz2svyvHHlHFXfv2W/mrPmETJmqrREgKyEZ8u7Q4WMbNm7jvGjZdpYWX0XqRd5xOrm5H6w2q8ufHz99IpAz0v0NA8wEjU331CkN8GSFpifvtmzZ2eXv3mxW9i08PELC2czPPVdEVFS0RNb9VeOFixt/6dFrgLz11fD5efJOOj59Cne5M+TKETNiNaSWRaLl5eUZ82eBKiacWb3urxgGVt1VMXIuVU+wNp0djYJOELccVAZQzlRO9sMH65PNz9dOlqIWuyqbttROYM1X0qkXeccJcmr5BQWckVGnCIa09HS8uClYipdoGGpOkGi4icBe4sg7CZaUqsoP5B1JYwhcXN3btu/i7eOL7JA8USWLl6589PgpGWJ43r1nPyFOzi5Ei49PJGd79h7YtGU7g/fpM+evXb8ZERH19KljcIj+Kq5w2OaYm/tLdhGHho6AQA3gRTVeunQtOzt7ydKVXl4+yFsic9qomXv3HnCGx2xP3rh5h5hkgDRRdcjByMjoHTv1xWCrJe8Qajdv3TU2cwnS1Izw96mpeuiXL/7+bwh5+vR52f4g8q57j/4HDtpMnjJTVpctKi420jFrXD8/fwmM+aprAY1LiKvrS6PdU7CEyINuwHEfODwiBAWJFynm6vYSr4dlLXJmBrhlM9ykcPvOPXG/T02bM3cReZNt8NBR68roWhpTUvLboOAQ+QtuAs1nkZ2TgzgzEsHiu7i+NA5HXSOLxU0g7UECMzIzqXHcVLHsdXzuHB0TS85xyzXqqlNfV++odOQ4VYADaIfdevQj3MfHF6tBfdFHWrbuWPsXPw1Ic8bMebX8Thdddcy4yT8cUSri2vUbV65c1z0m3r9PbdGqg6/fd6s9Hz12Qj4nWr80yat3mzZvnzZjbsdO3XW/hW3b93Tu0ss8gSwLU7LVazbonkq/TFoR7i891q7b/L8/2pqX72rI0B9/hrxjpn3u/CX6PoPd+o1bUlNTUR7jJ0w1hEtVmDVnITZE91QMFnLS5BlIPN1vWTQUe8sIy1ijB1XAP/8cnGGZV1cOJ9KydacTJ8/8MEGP156t2nQy7ibVgIePnsgIi+PS5WscOiYmdu36zbJcCPMHRhNZW1tAyfUfMGzL1p3kDbuampaG7eVfuIODQ22PazejakANKr0x8gN5RyFSlKg0hE5hYZExNlBDhMhFC4sU1WQmXlphSUkpNYHW5l9SkaSAOOPX8tdvkJrsKi3VWxV/xAv8EeQQkggYRwRzmgTKv4xZRbXkHWN/HW6ot0FDRnbr3s8qvOFs2KCRo8YbXl9f/3YduhreKm5ubi9R24b3/IXL69ZvMby12fRaqQL1eHOWcfGh5QKqsP+g9p4KZpfZiMx/jhw9Lsu0mhtwbfD28TM/h1ozyFWN5V1CYpIhWHGLQ1i+cq3VaohRUTEy96hfmqq8u3/fgRa4a7d+VR6re/HyteEjteVwJaRcmEG1aNVR99RI3sHlK9ebubyjg+/YuVf3fPmycvV6ZjgEvvb0llscVeS+vcO9+9qbkpXj7OyCdT1z9oLutzyis2Dhst179Ad7Kocur7sq5b//a4280z2VQoL+/vqna6rLiRNn/P21lyMprr2Wj1sI2joVaZqkY4Z/3e6WreUJIoGyZaAxPgKZnZ1tbrSHDh3NyNCXkasWSt41Yqol74iJca+rbd6Cpa1ad2KGMX/hUqtdDWTDNg0bPnb4yHGTJs9kGzBo+NjxU6zi/HDbtGVHvwHDxo6bIon07D1o6bLVVnFqthnm+Ic0BHkny+LIdUdslqHkRN4Rwnw0IyOTiQdu2VUzap8C1EbeafMty9kxGfv9z3bkxzjZY7Ynvb19zSGyBqq465GmKu9oUd169Fv19VJcbm7u3Xv2VvIuKiqaEfTAwSMRkVF4qfcRI8e3aNmBSbJUjcg7GtXBw0cZ7QxpkpiYxB/ZXr16LSHADHrZijVs6Awl77r37J+enmH0R7nyVN0eWnV5R/qt2nQ2Eg8KClm1an3l8o7Ikp+6knccXDp4jeUdE93xE6eKm4z16jMoPeNbGcqSrhhSJskdO/cwwgMDg9dv2FKRvHvl8Vq+sFddqlvpjZRGL+/KHeOr++xdHUIfKCouZqvxOPoLoPOQPcmnXPTWd1QH/lVSUmIkoof+Qqol7xj/aBV6UK3BxGCdKcN7lpWKhLv3HmzesoMGgFvkHdbw0eOnSEAKSgr54qVrCYlJbHZ2tz5//vfa9ZtkzNfP/+pVO0saOvn5+ffuOxANe+frqy1VOGXabBGUeK9du+Hl7WN7/NSTJ8/y8vIeODz+55+DMbFxxN+ydZc0vJSUlICAQEK279hjPB4q8s7P/42Xl098fMKu3fusRJjzC1esfE5O7sNHT0+dPk9IaGg4MqKkpJThBBlHCKkh7zgjo4UTfvHSFULOX7gsCa5ctf7kqXM4cnM/bNi4NTQsnJy4uLrLQGiGZOxu3GbvpcvXZAZfh9RY3vEv/qt7Ghgy1E2cNOP3P9syFuJmpKRSzPKOxkkdxccnXre79X//bUmx05DY+5/fWk2dPltmOyLvaBK0NBxyfSg6OrZzl54vXrjRAhExVA3tFkVIhJOWh3T79h/atOVdVaqe1tK2fZdpM+bu3rNPTB99YcrU2c+fO9P+9/xzgN5Nk2ajkKXjOzx8/MzR6emz50ePnTh1+hxax5B371Le022J7Ob2UhbFMIO8I/F27bvKcz7g4+u3bv1mQ95RffSvqKgY+6+3eqmd+/YP6Y8OD5/IF4eBkKDgEI5y4ICNVa8HQ97FxsZt2747MjKKmBcuXpGVZYGj5OcXhIdHHrY5JvIuODiU/k60K1fsiotL6OlHj51ct157zueVh2dERKRYQoPMrKxhI759mQZT2a6D9vXCffsPGzFppZTn3937YjYlxNPTe+u2XRXJu+jomAOWOyfVpU7knbl1NUwavbwbOmzMpMn6GxUG9SjvFL+Mask7folfVwoPE2Nv/xCFtGfvd3PoOXMXibkXeYfDw8MTk2TZqbkxythWwNIVFBS+9vRq3baz7DWD0fT29hUr/Nry2sq585dOn9H01szZC7QYlq9PiiMuPqHzX9pCQtDl797yJRUHh8frN2xlaJkxc54xMIi8Q29hlUgcg25lglPT0jDQ/AtJKq/IZGRkOjm7WBwZxrtHf7RoJw4BeXfdTnsbesvWneg8HHFxCZu37rDs/LJ2nf6qDQMVckHcBomJSUOGjeGIDBWTv74aVVfUWN5RYliVhqnwZKijVdAI5ZpHuw7aS/eGvKNO27TrIuWf++FD7z6DjhzVToS9ZW/OSgo4pDeNGj1x9JhJlv3a83zEz8vPP37idK/e+gfvm/zNWUwEIZVXPc11wcJlbJTVb7+3Kba8vrD/gM0VyyTNze2VvOsGg4eOZgKDo9NfPSSE0hOHIe+mTpttMQmf/QPeoNplrwG9hl/tATVLP5L30gx5R2dHhOEAppqurtptBLqhPKoPkyZrCWJPyKocZer0OQmJ373jCIa841xGjJoggVgY95faFVzkvtNXG8JsAXmHlevZe6AkePmKtjwZuzAULVt3wnBdunRNIptJS083fwLAKMORoyZYylB7zkpug8TGxi+zLKS/Y9deolUi75KSkv/q2kv3VAerSq8ZcueqISu8pnD1jsHbSuEpedccqK68E1NeJwoPEyPX0rDREiKgnyqRdxhBmiWTYDbsMpYRebfn65KBZhBhiJ6+/YcOHzlOLpK98ngt8o65tXx8etES7X1qQN7NmKk/QG3IO7KxeOnKS5evvX+fZiXv3r9PZX7cr//QO3e1L+ZZYXPkONKBY9kcPY5tJf8Snp+fz9/FXVbeybN3hrxDFBqR9+3Xp9flyju4devuzFnzKZMGJe/u3rOnshqgwjOGOqQDFUGjkldlzfIOx6Aho5YuX7102eoFi5bfsnxRhkAreWcoOXZJb+rWo78RuH+/9vW8Dx/ydu3a13/gcAlsDvJOrvRXUvU+vvqTpiUlJW3bd5H3ihyfOxvyjn5n2f9N3o0ZNxkrFBEZNXb8FNllyDsEitgEpHbZlzNE3kGrNto8UN5ZNuSdm/srs7wTmTVoyEgreefl7UOzkaPExMSRbdlrYJZ3y1eslUBD3p06fc7Ihsi7ouLiQYNHSoJkKe3rVXlHR6f//NbSatIoWMk7+eg2FJeU/P5nu6Bgba0MkXeATOQXS8Vvg5V3JCLao1oKD7NPl8S06n7tfkjRM0cnR0dncyDRHj169viJY3R0LOXm+PzFvfsOzi/cZHypIk3k5qxVKddS3pV+/GjeKGh9x0+GA5kr2Ay7qlKvvzjD1UJOoaITrAE1kHd46ZO1V3iYGPOrFQaVyDtf34A7d+4b75PKW+7IO7GqVlhe/te+7Ofi4ubs7EqIIe/u3ntw/vwlNnmmCsqVd/sPHJYMYOgRK7LAkMg7d/dXZBIFRgGWNcR+/gHTZ8578tSRWfiZsxd8ffVXYsvKu4TEJBknaiPvPD29jNtMyDtZFL2uqKW8w9EAFZ4x1DGu0w4PHjoqXit5Z3v8FH3NgHACRd7dvKV9OxspMH7iNBzALiXvBJF34ii36ilMWblTWLt+k3xVqHJ5d9jmmHRbwzgb8m7CpOkSAuZ1FQRDVzEPJBFZfcKQd/RWs7yTmmWyZCXvgkPC9u3XLyiSn+rKO5K1knekgJWQEJD3G4qKiugvc+YtvnnrTtkxCHln2ATKUD5iKSxfuVasmSHvmGSy+flpd4Erl3cLFy3XPdXBqtJrhiRSXYXH5Bk5Gx0do/stYG+ZiZmtcVBQ8MhRE47ZnhKvfAdf3FWn6Tx7Zy7l2sg7ijg+IZEeEh4eiYMNA5qcnKzv/mmkvH8/d/7il688dL8J+uphG1sqWPeXBx3m+fMXd+/ak2F6xfkLl40nruqE23fue3n56J4ageEYOXpC7V//NKiZvAO6ZW0UHkKHnoZWs7wb/k2tlpSW/rPvoLxFPmTYmFeW615M9JcsW4WYexMYjGKbMnUW9QJHj53ERL544Xro8NGy1haVwwhBNCwmCo8QpsWoKFJmnof9ZTt67ITUSGho+JSps2m3RGcuK2tGLli4LCo6hhQGDBqOEg0MCiZ7M2cvQDhu2rID0cMuyqSsvANGfZGDDOT+ls8AQnp6BiMEieDu1qPfy1evw8IiyDlZ4qzd3F/hXr9hKxN9QjC7+/YfJjI6csfOvZwFex8/eXbpynUru3/P3kHa6v4DNhMnzyj38l6Nqb28g4am8IyhLjY2rmWrjrQf8RryDlNw6PCxDp26xcXFU/LrN2x590572K5v/6EtWnWgQDAUeMuVdy6u7j16DoiJjaO30oa3bd9NfaWmprVt31XeeaSOmom8E3fZqqd4/+7e17KKpEbvvoMpZMKRa+fOX6JDPX7iKA8y0qoHDdYvpNEvpNva3bhF6ZEIneXqNTsc/OvI0RNEjoqKOX9B12oCe+3tH8p1Hbowx8XBoVetXk/VUEeYoPkLl/FLIJMrmUmSfw5B3pycXcaMm0yXJ+aAQSMYYjgKs0R5ZNOAmMg7hjmi5eTmLly0gt7KSa1dv5kUOGJWVtbJU+eIFhkZhbxD6hGIhbF/8JAEXV3dX770oKfLUxzQq8+glJT3Vj2dchg9Vp858PfuPfobZdi1Wx8OTeDDR08fP36GIzQsvF0H7TuWxEFlImcxhSSY/PYtbdVYAY4Z7I5d395irjpWlV4zjEQYX+g+RkurHGwLMvfy1wmAgLxjetBvwDDdjxpOS0e5Nm55R4lQLhRTjTdz9zMUXm3knWB0S4Fyr2jVscTqvAxfOfR8qln3mMBY7P3nIB1P95dHaWmpXN0Rps+YS7fRPZViXji6Ejxee7q4ajqjNtCsMRC6pwpgp6yq27xNnDTDKqTcbdr0OcwOrbof4TVWeJmZmdExsdgvjOmnT99MGOoK8cpRMPFEePsuhUAsFKPsu3cpIqSwYhmW7xpjtvBihRlKc8vYBYwXJphomVlZcjmQmDRvZt6XTYvPMVpQy2/fpsTExGKRyRhDfmqadgcE64yBIAWOyx/xcmhyha3HOEriclOvLKQpspVsf/yo6z/+heiUR63JXgpYbp3g5ugoAM4IVUdmKBPykJiYTGQ5KKdBb+IPlIncazYgMn8hM+SQ9D98+NbpqgsVaq53thmz5jMntgosu61avYFGYhVobh4NR+HdvHV31659bJQVJY+MI5BCRhxLOBuFTqCHhyduwo2VI6hNQq5dv4kWSEtP1+Pv3n/x4lVxM8AQLeBNoHjN16djY+Ml8NbXDNCw9X0NABlfrWqQbfLUWVSu4TaHi4O95nDZzFWPu2zV00Ewm3Q32i39TgJp8NJf4uIT6AiUD903Li6BHkdlmYfzjZu2U2UJCYlsYkCyc3K0zp6ZVbZQiSN1DXJtD7PDoINBEINGv6OngzFg0eWJQ3Jkkt73LkUzRARKhsmPRDMgnGgkToJ0c9wkRWY4yvv3qWK4sBVa/rKy2csmtZ+VrZkREYucEcZHS84y8Sivp38aO06/MQ00QsmkuQylTORkKTpCiENhspE3LCGTTI5O4lJuTk4vIiq96lERg4eMsqr0GmzmOSHZI0RqpHKYrnMWVhcdGffT0tL+81srilFC0NNLlq5q3PKOblm3dlMU3gsX97qVd9THgkXLzJdqDJgT665aU5G8e/XqNfIOxab7ywObPs0UgXKoorxjWqm7fj7VlXdy4U331DX0xhorvPoCi7ni660TuH/fAQure5o91KZx3aXOaWjX8BRmqPpaGvxKKFfhVQt//zerV39bUHr79j35BT/l+zENnMCgYLdqrl1fOYuXrDSkYbUwXlCrQ9B2P7yGF/AmyNvyYf0Ro8aHhHxb0kXG/VWr18vtdSdnF+TH0mVK3pWBNMdNmFZL1WIl79DUFC7zDyrm0aOnfv5vqCTmT0xuWrXpxMym0HInNDYu/saN27Tgs+cv4S0qKnrg8IjIrz29EYjEP3nyLBMRQl68cDPaJUmRDtuZsxfKlXeBgUGLl66UWRTTJsawiAjtq6+yV2BmM3jo6JGjJ5w4edbd/ZWxF8vC4Xx9/WlYeGNi427evOPm/vLM2YvEiY6O2bP3APm36iRMnsghf/Tx9Xv67DkhZ89dfObohCMnJ4dZ/r17Dwi/e9eeDAeHhJ09e5FzJ/7jJ44yaWMW+/KlByFPnjgyUdMS/SrvmJARyK6Hj56YC7ksP1XeQWNUeNoKCA8eUaFsNZu5NlV+qrwDBMTP0xBVQa4QNJCtKhcqfhk/Vd5B7RXevfsPrly1w+KxJdbpA6aNixOWdZd0Ty1gRGZsrXFJ/gx5Bz9UeMdsT3769Jlt85adZoki4z5Dp3wncJVlMqDkXTlQuDNmzd+2Y7furxEVybvrdrfCLF/azc7O6dxFW42iteWFJgFVNGbcZKLJLa0jR08g6XC8ffduyLAxOObMXUQcHLdu3X34SFvaJyYm9uo1fc2zYSPGlivvzp2/1LFzdwQibjRZRd+EYVfb9l3Y/mjRHqmH188/QBLEbXv8FNng6KPGTKR5yS0/OHjoiDjMcL7yxC5/lAhMvDZu1lbiuH37njy8b3fjdm5ubn5+AZEHDh5JTAK3bN350rIaanBwyKAho3AsWrzCWMJD5N3Dh/p3adzcX2VkZMqucvnZ8q7eB2xFHfKz5V1YeITVI3qKBsLPlndUPYpW9yhqAQpG5v+1pLi4pDbLZP48eUc7oTXq/u95/z6VIe/Dhzy2q1ftOv3VI9OyLBEY4/7kKTNT09LCLU8MK3lnjZTv8RNnatnby8q7fv2HIkomT50VHh6JpsnKypb3p6zknbE2GPQfOCwpKZnINMR797UVNIzhAXnnYHm05bnTCySjBJYr70JDqdkotOByy1JA5ES+AGsFhzZeO0cy/tGiHXLw5KmzzyzX3uD4idNe3j5EW7dBX4dMKFfeoVCZQKzfsHXb9t3ysV2aJuIMBycizZfJhzw0Rk9bsXIdDkDeub/U3w7h6OfOX+b0ZQU1EHkXH584fuI0Eke2GiqzXH6qvKOF1PnUQlGP/FR5J9quknm5oh75qfKOqv/hTTdF4+In3ZxFezCg6/4yvHkTJLfgICcnt0WrjoGB+mceGfflseNjtid79h4gja3RyztOY+iwMRRKjTezQTfKt/a93UreHbax3bBp27///rt46cp3lkfmQZ5RFXkXFR1DzSGeRAMJkybPyMnRX1yQCisr71zd3CuXd/MXLAsLiyDlSVNmHjl6ws+yYnhJSSkhEkEoLf34JlC7/SrMnrsQtUchmOWdv3+AVQ5B5N2zZ9qNVwPOJTQsnBN88sRRPqhqyDuOsu/A4YuXrly7fkOu2JUr73x8/Q7baM99A/JOnuIXeRcXF0/6ZG+T5fNKEqdcKK5yn5iWbfbcRVYh5W7zFixZuWo9Dj1RC0rbNT0YhidNnmmu+gULl61YudYcUu62es0GGolVoHleoQb4Bk7Zqpdt7rzFVK5VoHljL3GsAlXV1xV1OzlPr9HnZcvl7+59rSq9Blu52kP3lwElt3DR8uuW5d8ZhR2fO3fr0Z8QLy/vEssqMytXr2c8zc7OmTlbW+XqgcOjGbPmr1qz4eUrj/fv3zN8E//AARt5DKyK1LO8qz30PXGYy7du5Z2nlzcmQN4ZRLvIQlNFRUWyNPnvf7aNj0+IiIj8aFlwziye7B88tD2ufR0598MHcZSVd1lZWessF/xycnP79B1cVt4dOnTUy1tb/yIzM6tFqw7ymlJy8lurai4uLpabocLAQSNoNO/evWM2ICFTps6miMrKu3Hjte8Ayp1iA1IzltvdtUv7cvk3efcm6Nr1m2zaV24s35gvV96dPHX28ZNnOHbv2Y+8c3j0GPfU6XORd88cnWTxjhcubmnpNb/AblR95cglQPOll8q1XUZmpnEdtM5JSn5rTODqHRpDYFBw2ccfCwoKvLx8aAO6v9FCLTM8656KoWEYHVPAjBi2Ww3wjRdGhMorjr3E0T0WzNeAq1X1zOTDwr97FjYuPoH+xZhNR/Pw8Cqrdd6lvCec2bvu/wqBbPJyfSXk5OQa1xqqSMCbQEmcTQ+qmIruFFURT09vefvVID4hkVPWPRbCIyLJidViziUlpZLDj58+iYONmKWlpVeufls3oDZYVXrNMBKRVlSJtqsvmoi8syrf2sg7VFp0dIyrqztCBAebLK5oQI+KjIyOio6R8a+oqBiVExIaRjd+r60u7RwdHSufWIHkt2+JTCIkm52dzV60YPLbd/7+bwICAiWFjIxMIsTExLq5vyKCVa/GNJABIshjpPQEpCRiS/aa+ZCXFxsXT0xyjvqUQDJA/JCQMFkAhT9qObSsoyEROC4RjIcABM6FQC3nMbGiI4n22tM7IyPr6dPnefn5bIiACROnce4pKe9JMyUlNSYmjjhYNM6L8+UESSEvLx99jGkjkGgBb4KwSnR7diUlJddG69RM3tEwrMZyM5Tb5i07blRh9TUyP33GXOPqbFWg0pcsXUU16f76ZviIcQ4PH8u3Zc2sXbfZ7sYtD8vH0OoFesSQoaOrO3SVpfby7hdoOx8fv5mzF9g/eISDjWl6ub27IeDt7SvfN6sioaHhzBIrv0L/U2FQqLzuZODQPRYMeUfVGwupVgUbG9sdO63XYFu+cq0sFUma1KwEmjl77qLVdxdsj5/q0LGb7qkYTPQ/+w7JAnvVYsbMeVOmzdY9lWJ347ax5Hh1yc/P33/gsO75ytZtuzg73fOV4yfOMI7oHgtbtu1safq8SrsOXfcf0JdGDw+PePBA/7RubbCq9JohiUgTaoDaDpqCvCtbvrW/eqcoF7QsXVT3MBuLTzDWZ//11EDeVa7thMdPnlVF3tWMI0ePNxx517Fzd931Pcazko2dWsq7X6DthN9+b2N85yo1NY1pA4O3eOsd8xIezOjWrqveuEiDf/u23tZAZlyovPpk7NA9FkTeVVfbMXc9Znty4KARzMb1IAucvsg7Ly+fcuXda0+vHr0GiKAE6j0sLKJvvyHirZyXLz1+qrzDbFrd6qk6mzZvtyqKwsIiiqjL372t2nZUVAwCl5mAeC1XOj07ddY/0Qtmeffx48fdu2uoOM1YVXrNkNZlpT0aFI1e3skjF1blq+TdzyMqOmbWnAXX7W7KpofWB9WVdwz2lWu7vLy8nJzchw+fGPKuqKgoOCQ0OjpGTBLukJCwmJg4vKWlH3NztZvdWBzC09LSsMvyOGZpaSkhbMYl0pKS0uzsnMLCwrLyjuET08bfzfFxh0dEykHJAwctKCiUNDmuxOG4r197vUtJefMm6NFj7d562eOaSUxMYpdxK5b8dOjUjV/SkRDgdAhZsXIdv2/fpsTFxaM2+JfkhDEMNyOfcc313bsU8pCUlJxrWXM7Ni6Oo7A3Mio6IjLKeEUuMSmZP8pzBUSmPD98yCPEfH2OxClAKVtSyMnJMY6Snp5O5BTLbR1KGHdMDIey/jpnWWoj7wj5NdoOzPKOExwwaLgUeEMAW6q7tGZcUt3xvjHKO6q+WtoOnjk60Zhnz1lofLxLqIq88/bxPfr1WZHIyGjauVne0RFo8Fgw3W8hPz+fcBdXd0PepaVpfcS47oubjem3lcYCQ95Rm8SRDsVmNhpIMSxAXFyCUd0oe+LExmqdDpMlf8EdFRWNw+pJIQ46Z+4i3fOVW7fvUUSDh47mjPUgC8g7drlavqMNTs4unGxF8g727T+MEdI9NcWq0msG9qEhazto9PKuTbsuZctXybvmQHXlndUQbgZtYXPEFmmCndq9Z7/IOyyyfLUpIEB7YOXmrTvYOzTH3HmL+Q0MCu7cpWdiknbHfMXKtfsPHMYCenr58Dtn3mItUctzh4zT4RER2C+8zErXb9xS9urdqDETT5w8g+7Bsn/8+GnwkFGoJXSVYdQmTZlpfKziwEH9TefJU2eR+MFDR0+ePHvr9l2OS8ZklxxX3IC1PWxjK/PjS5evZWbqN8t+ePUOG81ZM4R8+vSZwpk6bQ6BGRkZlyydy8XFXRYpnTd/qd0N7fUg2+OnqBSs/6HDx7r16Jfy/j0FtWu3/jmEiZOmo54p6tFjJyFeCdm776CINjKPFCsqKnZ1e4m5Z1j67/9av3kTyFmcv3BZnh/dvXsfRbR9517GHiqUUYrAyqmxvMOAVHeArw1meZeenjF2/BROnLHw3n0H2t4Dh8dMD9iVkJDo7OxK3kQo42AqwohLoLePnwzklNWhw0f5LyV/2OaYlW1kCL9udwt57ebu8fiJI8eiWhlQHR21CgJ0g6+vv5agty8peLz2unHzDiO9jP0i72Jj4wjx/frISnZOjo+PH38xBDeZj42LJ05iYjLdqnHJOxrM0GFjqlX1lLzcMXz85NmgwSONB2OgKvKO3/4D9a9RycuVhryLjokV0+Hj68+IxoGYMlHFzJ0oZHoZfRPH2XMXT546SzRtoYP3qU+fOsrKGmvXb6a7aQmZMF+9o4MfPaY/nL1y1Xp+Se3Onfvy/PTVazfkS7UuLm7yVWu6fKDlsWnmk3917YWxovlZPbwEnEKf7y9A0oDFRNy4eXvmrPnmKSX9nX7ds/dA8YaHR9JyKpF3Tk4vSET31JQ6kXfde/ZvyNoOGr28Kxcl75oD1ZJ3uqcC/AMCT379tP/iJStF3mFWmNoiJoKCQtBSe/YcwJDhZaYuQ+nefw6IvMNyLV+pf0wCg2534zbR2Hr1GYSlmzZ9jnxNCPoPHF5W3pGOrJsDpDxw8AiO++FD3sjREyUQEy/f7IcVX99iWbVaM8eAEOGX416/fst8XNkLGZmZ23fsETfjx4IFy8T9Q3lnvokzZ94iyoHEU96ntu+ofQ5y+oy5xrVAyoffnNxco1JE3iW/fft//20pudqwcauXtw8ZW7pstcRB3r2wfFH30pVr8lIR2oLc4kChIu+Ki0t+/7Od/J2RjDFs+MhxiBIGuap8nbbG8u4Xg7zr0WvAwEEj2NB2ssIlLFq84vVrr8zMLAotMDBYFhhHPaPpcyA3Fx1QYLlwsmTpqqBgbdxlXO/eoz8O/oLox2EF4mDjpm2fPn+meBmhSZ+qXLN245Mnzxh0J0yazlhLtIOHjsiFKPOnyZB3RODAuMkwMxCmQwg4aQaubu7y3UJ3dw/UIY6CgsJ+A4Y1LnlXAygWkThoKYpU1joVqijvdu/ZJ+pq1RptPVtD3jGhMkwHu+gazHxkKQPYs/cAnSIvL9/oI8z0bt66e/rMefQZXk9Pb7OQEszyDolv1O+06dpHj+hxy1doK3CBt4/v+o1bcbRp95ekT2uZPWchIZzyyFETRo6ewHlZ4n4HbY9poe6xwJww0vK9osKiogGDhps/w4i8w+iNGz+VKS5eBu4fyjuxeLWhTuRdbaCIPn76RNfDLb9AiBFoQOFIOG5xiLsqKHmnaKzUobxDdsj0Fwx5h6VmfMKQAVaSafGWrbuwwrt2/yM9ENttyDvjleddu/ZhW+VfQMiwEeMql3ekY6gQbVS+dfeffYf4r2gmuH3nvnzQHQx5t2DRMob4vf8clLksx5Xv+gsSR0hLT9+ydYe4kXcTJ80Q9w/lHRkTBzAqxMcn6qmXlFACQ4aOLvj6YFZF8i4hMXHW7AX6fyz/4gQN22rIOwrf6p1xDo3+KLJc6tP/XFLCfx2fO2/avH3u/MXPHPUVfyqhEck7Gdus2LptJ8O5uKlfDw9Pcf+z72BYWDg6D/0nIcg7z68DbddufV69ep2ZleViKVsrrl674Wj5/AzF+5/fWknglq07L1y8woyCYTU/v4Dt8ZNnMnJbyTujSfzdve+7dynvUlKMh3H9A97I1e416zbJnTvo8nfvJi/vQkPDQ0LCMizfkh42YuyNm9+UDcVVFXn33OkFloEyl4tthrwzmw6qmI4cHRO7b7/+yoLIO2ptzLjJ0kEANZCekUGlrF6zgQkVXUYiG1jJu6dP9X4k8o4eN3+h3oUNeYeU11O3WELZiz1E4ZX76UsOav40NmB2oqNjpIgGDBzu8vVWLCDv+D1/4TLNhl9mbj+Ud9Nn/PjKfeXUo7zLzc09deqcq9vLuLiEZ45OISGhxvPrdGHMu7Ozq3gFmyPHR4+dLDeO+B04eKTHa6+yyx2Ui5J3P53snG9PEdULNIUm+WXSOpR3cXHx8vgL+mPh4uUi7zZs2iaXRgh88cINISIX7bZt3y0Viu0uK++w1O4vX4nb29uX/2Kb4uK1+1ZYxh/KOyy4rHvEHydPmRkbF4+9K1fenTt/CVvPJl6O6+auf8/Rx0c7rriBBrDv61tsAQFvzpzTn9eplrw7d+FyvGWkIeGnTx35Xbx0Ze7Xz5OLvCssLERu4qCgevYeiLzLzs5Gikm5iUSm6MrKu7PnLjk81JbOATd37cYQh0Z/8Jcly1bLN2D4Y3BwiL39Q0usLzNnaaVUOY1d3lEIGV/X+qq6vIuNjWPgP2Z7sqionKVtKpd3EydPlxADkXfyF0Z3o0lUIu8QFs1K3hkzQ2DAXrJsle758uXQ4aMi76i7cuWdfOMHWrTqeOSIXraGvDt0+Jih7+lZdOS0tDRD3tFokXcoQo4oV3ToI6glWb4eL3uRa5a436hc3tHjdu7W+4LlYRJN3q3bsEWeoyXZFy808cFBFyxafvfeA5oufVyLbYJsDB85TvdYkE8ACKgZ4ztGEPZ1XZg/W7Z3cNCMwA/l3UXLF5VqQ33JO4qX3uHr56/7v3y5eevumrUbxY16e+3pbfX+9aIlKygB3fPly//7TwuEsu75EUre/VyePnuOyZNlR1atXu/urg/81YV+ZR6wK0fEh8HKVesNi1AR6zdskQekGhF1KO8oMQyQg8Oje/cfMP8eNnxsenoG8+ZTp8/dv+9w6/a9iIhIbCUOvFeu2lEdKSnvhw4fc/78JQZRxA2iMDhYG9IYAolANDYGRWrtfWrqffuH9g8e2t243apNZ6uvhmjT2UHDGbDl0SUEJalpf7d/OH7CNIaEtLT0EaPGU4NoI8LHjJuCyuEoBw8dHTV6IhvTO7JqPi5jjLm14LZ/8OjChcvs4izy8wsIwd2iVQerdaQw4swshwwdffeePUchY8gpmRtkZWUjKPnX3bv2ryxj0pvAIDe3l6mpqbdu3ZVRn2I8deoscVxd3bFZsy2f4KNdaask3negS3Lo155e3Xv29/TyJpOzZi/YvWc/gWT++Ikzd+7cf/T4aVJSMgMGU3zSJJ9eXj78kb9fvWpHUSxYuIxB5YHDI0PnVUJTknehoWH37j0QNxZf+x5gBfIOGMLLvTMLlcg7xp4Jk6bLY1uhYeGyTPqmzdv5ZYLBb1l5V1hYtGKVPt94/vyFNAxyIjOcpOS3TVve0T7PX7g8ZNhoEWGyzigD8ImTZzELKK2/uvZyfuFKtOUr1gweMorOLn8UEhIThw4fe/uO9nQd0ehNRKBtd+7SE3NheQA39fTp8zR+uq2zs5YOHQpBwFzI3t4BlTZo8Ej6C5pA+ggWgLHm9u17xMHLtFOmRgZU/aAho+iAN2/dod8xNm3YuJVCwNoMHDxSHgAlGw8fPSEDu3bvIyYiDz0hltDy/HE8mcT2isho3/Hv1Ws2Rlk+U2lAPhl0xE2jOn3m/MTJMxCmeLFU23fsoYjOnr1YWvrR19ePEnjhoklG7TNIubkxsXHXr99s3aYzGaAM+UX4jhk3WR4HBCafEr821Je8e6mtm/vdZ6IYaI5+XZtWLtGZ30qhi9Exlbz7jgZ19Q6DKPKOXmTczKouRUVFdGzd8yOsjhIYGGxzxFb3VAATd3q77mkk1KG8M4NtEnR/xRhxqhgfyo1pFWi4jTBzBHEw0vv66lPA/Px8mXlXjukIOkaaZiyH+u5YlZCXl5ec/JZoxqhvpqK/m9P/YZza0CjknZe3DwMb+htTrgdZQOOePHX2xQvXTMuLFBAYGMRA7uTsIl/bDAkJJUJcfIK3ty8zB0enF3KJBR59v0q5AWMthvHmzTslJaWkw3FDQyMSk5JQGKgKudHm5xfAUfwD9KvFCAItD5YHIpnk2NndIkJgUDATG8YqAhnsnzu94C9Gg6TuwiMiOC8//wByaL649Yv52fIOaKdWbVVr06YWbjgMt5kqBlYUTdD9VcP8Fyu3OMyYI5gxAiuKYF4wyyqOeI2Aby4T5kBLZN3LBKPscno1oL7k3YYNW+/c1RfBMaAji0PkHYL7zl1tNgVoO6S2knffUV15h1lMSXlv3GkC7OD796nmyRYCiDiGciosLLT8RXvKGDATWEBGO3kTEKgnIst8V+QdXrG/xMnJySVxHPJOnEA05rv8Un9msUVSt27dxfqTgiHysrKz+XvZF6OIzIQyL4+5vR6ZeT95IKb573SY1NQ0Swa057QsWdVyYjmRLPOJmKFYCGfCwd9Lcb9PlaJj+/RJT1niMALpcVLeS2pyIDKAIMjK0jKv3WDKz2c+h5vIlAnunCovFPyT5F2jICQkjPOiuNiCQ0IZp/Ud9UFCYiImKb0WHyD5STSWq3d1CFKMrmQMD82cXyDvFOXC+Ghc1q1DQsPCnS13h2tJPco7T8/v1oUxI/KO8bFHrwEMlK9evcaoKnlnTbXkHbJjwaLliJLz5y8bzwEcPmKLnvN47fXe8n0qx+cvHj1+Spz9B21Evly5dgOHk5PLvfvaHRMqY/GSFcxT9/5zkBpCqZw4qc1ZHz1+xmwDwYTNbdWms7xnfvrM+RUr16HPSKFj5+6ie06fvYBUj42NW7d+MxPfRMv7awavPDyzLZ+dEGyPn/Lx9efvJGUoNoPHT56ZA5F3Bw7aeHp5kwfj+Y/tO/ZwaqQwddocIrdu+5eMcJYTWXnq9Nk9ew9wIhJZiI2N79qtD3/x8fFLTEwmteEjx12+cp0QOHvuInEomZ69B+INCQ0LDQ0nzpBhYzgoMs7HR5vcb9q8PTIy6uq1GzY2tpcuX0V0rlmzUSqLZOfNXyqPu1WF5izvgDkfrUU2q5r6xTA3oIVLM25QNDd5R3dbv2ELU8HKNU3zQcm7eiQ8PKLslaraEBwc8vKlB4OL7q8F9VXp5y9cNh7D+Pz53/T0DMa+c+cuyTfZRd7Ri3/7vbW3j++7dylFRUVK3llTLXlHyQ4eMqqgoDAwMPivrr0JQbE9faY9Kf/CxS0iIvLT58+rLe+rw6o1GyhxvMadlA0btzK84ZhqWkUC8WTUQe8+gxj8cHh4eMpX9kNCwmZ/vb/+vz/aivjoN0Bf+mj6jHJutJnlHQ60l7gD3gRu3qI9HGOmrLxbumwVjQamTdeWLkMaTp4yk/EYJk2ekZae7ucfYDwTYF4Ow0z/gcPlTdLg4NDJU7VX0+1u3JEVEGDVaq2Ips+Yl5aWThxOf8DA4YSg7YyllcB4bZMJijgSEhKXLF1FhmURtarTzOWd4oc0w6t3CjNK3inKpb4q3c8vwFidVEAMGEtcibzD4evnf/bcRTFKSt5ZUy15x2zg+XPnmbMXFBYWibw7eeqsyDuBEp9qUUUGo8dOMuQdIuPDB+1dSPPjRzNnLTBeaTTkXYjlJqnmCAlbt36LZec3ebd5y4779g5ZWdnz5i+RXWYMeYeU1Na52LZTwpF3EyZZv+xmyLtCywKbcnMWhyHvHJ87HzhoI0sZsREeFRVt5N98Ima6dO2NFJa/yEd1b92+J89fg8i7YSPGZWRkShxZ3hN5Z14bac68RUS4cfO2+f0gOltRUXF17ygpeaeoHCXvmjlK3inKpR4r/cpVO+O9E2AUM4ZC5N3jJ5rwyMnJ6d13cKHlRpaSd9ZUS955+/hu274biYMbeefj42eWd+gns7yjxBFY23bsMeTdth27pRrMqujAwSPVlXfIIDf3l2xl3zMHQ955enrn5ubKTV4ge8aaFwaGvAsMCua3rLwLCgoxXrZAmxYXF1dF3i1astK4CSjLeJaVd5u37jSeGvT21lSdlbwjcTlN3W8BXWtre0qWgKo6TVjemS++1glVTLDOj1u/KHnXzKlbeYf9pIPISGFQNsRMJbsMKkqh3HBCqthJq3JoK0hZ0P2VUoP0rYj7uuYLyHlZpSmZKXsgCTccQnBwiPHBxh9Sv5o+Kys7NDQclRYTE5uXlyevNEFAQKDjc2cEAKcsK2Elv3375k2Qk7MLkT9+/MgvETy9vM1PalWCknfaIwKr12ygNO/cte/Wvd/pM+cRcHv/Oejr6488evMmkDjXr9984eJGnOCQUPZGRETeu++A19XV/cbNO1RGcvLbWXMWECLiBu/Tp45+/gGPHz/r0KmbPJd25ep1WhWVd+78pSlTZxUUFFKdv/3e5ukzbZ0CV9eXf7RoxzZ3/uLnZdYoiYmJGz9xGuGyxvcx25MvX3pwOMub8N+9bA/Ix7Hjpjx99jwhUVuojDgcNzIyCk02ctSEuLj4kpKSXbv30UpI4dKlq0VFxZzF6jUbmTToJxLwJr/MS74oy6PHTmpn7fbS1c1de3h25bpr12/m5OQSOGPWfP4bEhJ65OgJvC9feaAyiTNv/pKTp869//oxRFKQ00RSGxKZvS1adeS8xFtFmqq8ow8/eepoPAZae9DuNIYfJsiEpG6PW+80N3mXYVk21mojnP5uDqG/S3wrYmLjdFdToa7kHSOru/ure/ce+Pm/8fTywdDZP3hIIJb/mO0pY8lGK65es5NlICsH67dx83bjrWcBY8gQs3HTdqlBgXpcv2HLnr0HIiK/W4ikLK9evT5kc8yYjVcFjsjZTZk2++q1G0mmD2+US0Rk1OatOxAfur+aUG4nTpwx3gv08fWzu3Gbo/v4+MXExj548Eg+u0fJHDp8bNf3fZPxi0wSubi4hN+Vq9bTeYNDtK/Q3r/vEBGhfRvjh9SvvPtlKHmnQbehr9K+xUEIHQOHuIHmSPfDa3QY2Usgu/Dy36LiYkLEC19T+ES4rC1JZEmfX1o2MXHj4Dc+PsG8zPfBg0dQkLrnK8QkBdFARn5A9lpB5JKvkeW4kp/irycoRwc5IyNvZU/EjFWxEJM/SlJli05S0HJSoufEw8PTUK5MWZavWItWJn5qapq8oVItmqq8szliy1zC+MZ27SFBlyokqB33RV0et95phlfvOnTqbkgK+uDceYulG6JORo+dLOFbtu60s7tt1cETEhJ79BpQuRhqdNSJvKOg7t6ztz1+Svdjx157Mk+meHG/9vSqaOWXhMSkKVPLf47ZiiNHj8sYYeblq9fHjp9680a/PQLU0cxZ81d/XQW3ckjTatG7qvD//tMCkaR7KoX0rZa7qzoBAUGyRCK8fu21eYv+TR14ExhEOzQuB9y3f/hny/YMLuKF2XMXde+pfXZPGDZirPGRQ7D5ujR05Sh514iprryrd5KSkidMnC4CDk6fOV/jntOQ8fX1X7Nmo5xjbu6HjZu2oeoeODxycnap6IpCJTRVeTdqzES53V8zUNJv3gS+8vA0Ehk9dpI8B2mGQYuR45XHayPaqNG1Om4DpJnLO2DiJI9PmOUd9Ok3xKpTODx8jBb09tHflGoa1Im8o6D69h9q1TX2Hzj8C+RdUvJb41v7QE6wmY1d3qHVevTS9Vl2Ts5//9faWF9M+GffIUPeBbwJPHj4qHF9tLCwKCwsohJ5d8z2mwqvBCXvGjGNTt5BYFAwQsfL24ct5utze02PFy6uL196yGl++JCHBfTw8DQ+vFMtmp68S0/POGZ7csDA4WfOnteDvkcWCNQ9ZcD2xcXFn79wiXQkJC09nQQHDhphThCpR4GfO080faW6jIxMovWv+LiNlOYs7xB2EhJpGYOt5F3f/t/Ju4KCwtu37125anfp0rfVaJsAdSLvHJ87mz8kICQkJFrJO+arAQGBLi7u4eERTF+1OIlJk6fMoktmZ+e8evVaFr9A3Hi89iIaMkXunEBF8o56QcpEfr0Ve+LkGRqqIe/QQK89vd3cX3EgCcECREXFcLjo6FibI7Yi79LS0t1fenBEP7+AT58/u7q9xO3v/6bc+5iGvHN1ffngwSNSI7Knl/6VM2BsIv2kpGTSF3nHiXt5+RCNzGCdomNir9vdsn/wiF3PnjmRjpSGQVpaGjNYcYeFhfcfqC8ZYZCZmWk8gI68CwoKNl4+8PN/Q6FVIu/I83W7bx9Aqwgl7xoxjVHeKapLk7x6d+/eg0uXr2El8/K0V15QwPn5+fJqNoHyEC6/JSWlBBrDAzx3eoGtlKctzdy5a2++70+0UKJ9/U65wb37Dy5dukqCclDzcUEOJ+MZDmy95KSB0zzl3a7d+9AQVssMIe8YUwkH9J+fn795noCMQG0wKxg8ZJQe1CSoK3l37foN3VMGkXd0HLSXXIWiJKXHobomTZmJ5MKNxhIViP6QI3bu0pNKwQEVyTt+3dxe7d6zH8fJ0+fogIa8Y+aGwKK30hnp4/KMOHM2+TB0fkHBhEnTkXeZmVkjRo2nfgmcMnVWeHjk8+cvcKPJ7pV3lc6Qd2S4W49+clNlxqz5oaHhOERN4uC4Xbv1Qd6RpZGjJ8rFtn/2HaKscOzbf3jIsNGouqvXyik3bJR8aRowWUwsxV0u8mqjcRH0uGXd1krknbePr9lbEUreNWKUvGsONEl5N2XabCdnF7QdM+AHDo9c3dyfPHF89Fj7xpSnp/fVq3YMz0ixK1ft/PwDmNcSjqXGnlb01tikyTOdX7jgINr+A0T7biZtMHX6HCdnTfkdP3kmMip6/YYtjFK7LB8Xx5TbHNFuxDg5vXifmhoaGoZxZECSPzZkmqe8Qw0w6BpLYwoWeTeJ8Jyc3KHDx1ot/X/y1LnU1DREABLw7r26XIq2fqkTeefl5TNi5HjdY1kn1dPL58bNO/L2osg7xJlhjt6+e4f0wYG8M9+cXbhoOVKMbnjl6nX/gDcIr8dPnsmuSuQdf2nTrgsOUU6GvPPy9jEWyH3y1HHXrn04+g0YalySHzl6AvKOaOMmTM3PL2BD/CUlJS1YuOzSpWs3b901boCaMeQdZzd3nr5EF/LO17JUwuYtO4yFTtu0+wubwGTvry69JP3Hj5/5+GhfrYV585f8/mdbedHQCvMUNC4uvmXrjrpHm7gWBAQEMr9N/vpuB/KO39VrNsrnlR0dtUKoXN4NHTZG91SMkneNGCXvmgNNT94hpEaPmaQ5SkrOnL2YmJgUGRkdERH53EmbcEdERDH1Z4S+c+c+k2ksozx2zQAQGBjs4uIWGWl9t4XJ98hRE8RtRCt7U0Y77ljtuIzxdje0h+7lbsg0ywrbzLD9/d9kZmUxzmVkZKAPrl2/WfZhvgZIc745a3z7SAZm883Z1NT04SPGydKV8Pnzv0wPaGlsk6fOGj9hqvnCXqOmTuRdfHxCj54DjNf8heEjx9ETcYi8owctWPhN3nXr3heHlbwbMXLcp0+f7t174OSkTbf4C/JOiroSeQeHDh1Nz8iwf/AQtyHvkFnfybvd2hW+fv2t5R3RNm/99uJCSWkpJ5KRkcl/7e21BK2ogbyjzUgIGC2HQhswaLisCGYF8s54r4Kc/Nmyffj37xF269HP/Owdv5gdCnzSZH1J/Mrl3eIlK3VPxSh514hR8q450PTkXWpa2oqV+ocamWTzu2rNBmSZzZHjHq89o6NjTp0+iwFdu34zuzj9Bw7aAy4Gt+/ce/nKwzwpT01NNRI0uHvX3v3ld9HS0tJXfv1AJNKtoKCQoSUmJvbc+cvs8rKsXzhl2uwnTxwTEhLv2z/Mycmpky8//mzUqxUgd9Osnr0bOmwMAj0tLY3WhWqXG4gQHBLao9cAwsXb2KkTeQfPn7/YtPmbSDp95vxCy3cscdMdkHcUo/F8GKLnylU7HMi7v7r0kkD6761bd3FMnzFXHjXesm0X8k5efEGxlSvv5Jboy5ceaCC5wWrIO6Zbi5ZoF7Tg1OnzYZb1jM5fuJxk+ZplRGSUyDtUGnUtFxrps2TJWD3+piU/VlQu75gfPnz0BMe7lBSRdxQCM0PRZ0FBITLh9PHxo1GR7bHjp+C14u3bFOPZO2DquHqN/jQhkD2KyLBOjx5pNy5g9pyFrq768waVyDvOUT17Z6DknaKx0vTkXURE1MZN23Awy39h0U9ypwzbmpmZyRQ2KCiYgURm/7du301M0h+pNmA2/Py5szE5Dg+P2LTZ+pt1QGqMWNhf8WKm5biER0XHlJSUIOPy8/OdnF2Ki0tycnN9ff3Z4uLiU96nMj/29PJuFO/YNjd59+5dSrsOf2/cvB0H29t37y5cvIK8Ky4ufvbMadbshYh1iTlg0AiGak9PbxrMtWs33qdqoylNiwgM8Ewbqv7p54ZMXck7yM7OfuboxFTH29s3+e07ykeuVN23dzhz9gLp5OZ+ePrU8bWnt6ub/jIBIcTx8fXjXwEB2lUoSE/PcHHV3kJAXts/eIRpQkwjEAMCtGVszdjduE3VEA0Jdf78JUIC3gTeuWuPLAsK1p7KePfuvYuL2ysPT7q55R/aextUqxyONEXAEciAyBFDQkKpYvKPGwUmr3qY4Y/8y87uNjNJ2gZuxCLxr1276eT0QuRsZFQU0bBI7GWTQsBLmnKOiYlJ5y9e0ZL78uXsuYtPnjpaPadLa+zVZ5DusfAhL+/Rk2eWsvWJT0ggw1KAUjIcDDc5wQBmZWUT7fz5y/ySLL+XLl9jrDeWj7E9fsp4sawSlLxrxCh51xxoSvIOPefnH4CxlpsRtUcSRCPKl0sqQqK5urobw09TornJOwY8xjaECA4hIzOTwfLjp08FBYVshigvKCjIy8tDOjDWounxEshQTYS09HQGURnLGzt1KO8UdQj6r9z7trVHrXtnRsk7RWOlKck7RtOTp87cv+9QV8OqnqD9DxKUaPfuP6ir4zYomuHNWYUZJe8aJiWlpQcOHtE9dYfdjVtyk/qHKHnXiFHyrjnQ9G7OKuoWJe+aOUreNVhSUt4bjwrUCW/fvo2Orup6sUreNWKUvGsOKHmnqJyqy7uhw8YQWW1NbKNafyjvVNU3w62KY0djp8nKO+rPqkbV1sS2SZNn6vVdKdIY1NY8t8pHdwMmAAhBtTW9Ta/gilFV3wy3ZjLhb5ryDqg/tTXtrYojt0KhUCgUzY0mK+8UCoVCoVAomidK3ikUCoVCoVA0KZS8UygUCoVCoWhSKHmnUCgUCoVC0aRQ8k6hUCgUCoWiSaHknUKhUCgUCkWTQsk7hUKhUCgUiiaFkncKhUKhUCgUTQol7xQKhUKhUCiaEF++/P8cFIo7t6zDDwAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small-Signal Stability Index Python Script\n",
    "\n",
    "### <br>NYPA Machine Learning Project\n",
    "### ALSET Lab Research Group\n",
    "### 12/12/2021\n",
    "\n",
    "<p style='text-align: justify;'> <i><br>Electrical Power Systems are complex dynamical systems that are in constant tug of war between generation and load. In a simplistic way, both ends are constantly being tracked in order to balance the generation side to the load side and reach a stable operating margin. Since load is a stochastic parameter, there is a certain randomness that turns power system stability into a challenging and important task in power system analysis.\n",
    "There are numerous stability indices that have been proposed throughout the years utilizing a steady state formulation. These formulations require knowledge of the system’s admittance matrix (YM), generator’s active power injection/consumption and the load’s active and reactive power. Alongside these requirements, there are also modeling limitations and simplifications that do not account for the dynamics of the system. \n",
    "<br>For real-time power system stability assessment, these static methods are not helpful, with the exception in planning phases. That being said, the proposition in [1] is a step forward into developing an online stability index generator that can be used to help system operator’s in decision making. The development of the small signal stability assessment tool comes from an estimation of the eigenvalues of the system using time-series data from dynamic simulations. \n",
    "<br>The stability index is calculated based on three layers of information:\n",
    "   \n",
    "* 1st layer - A scalar (positive or negative) that grades the overall damping ratio of all the modes.\n",
    "* 2nd layer - A vector that points which pre-defined damping ratios were violated.\n",
    "* 3rd layer - A matrix that has more information about each mode’s damping condition.\n",
    "\n",
    "\n",
    "#### <br>Algorithm Steps\n",
    "Figure shown below are the steps that the signal undergoes to finally generate the small-signal stability indices.<br></p>\n",
    "    ![Steps.png](attachment:Steps.png)\n",
    " \n",
    "### <br>References:\n",
    "\n",
    "[1] Sevilla, F. R. S., & Vanfretti, L. (2014, July). A small-signal stability index for power system dynamic impact assessment using time-domain simulations. In 2014 IEEE PES General Meeting| Conference & Exposition (pp. 1-5). IEEE.\n",
    "    \n",
    "\n"
   ]
  },
  {
   "attachments": {
    "table1.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn4AAACICAIAAAD25nEaAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAABw0SURBVHhe7Z3RlSO7rUUdiQNwDg7DX87EPzcNx+AEXiLOx2+rDvsYZrGqSypJzS6d/aEFgiAIkoDQPWvm3j/9J4QQQghvJK03hBBCeCtpvSGEEMJbSesNIYQQ3kpabwghhPBW0npDCCGEt5LWG0IIIbyVtN4QQgjhraT1hhBCCG9l3Hr/FEIIIYRn0DprYbP1/l8IYW5SpyHMT1pvCJcidRrC/KT1hnApUqchzE9abwiXInUawvyk9YZwKVKnIcxPWm8IlyJ1GsL8pPWGcClSpyHMT1pvCJcidRrC/KT1hnApUqchzE9abwiXInUawvyk9V4K3ijP9OEkAeZEtSn+/Oc/N234VEiD1lkLab2/kn/+85+vrurkwPzkjWbj73//O4+ypk0f4C7jLZ7iJDwLnqN11kJa76/kH//4B2/017/+tY2fykudhyeSOp0KFQ7QgJtqeSM4Uk1PqbsU74TwIq2zFtJ6fyWUFm9EhfvXX5c9NKOvsvdULUhpJNffofnUFOTPyiaHN2pSmABVTe27oOICDau8X3eSU7wXgOdonbWQ1vsrUY1Rlq7MiorfNV/RlFa5ODVUbS9WjfzsPDm8UZPCT9PVVGUppttLHa+7FO+V4DlaZy2k9f5KluL6E/VJNUqmAtFrWEsUMPOUCrKaAUqG8gC3NUmA30CeaR66mjJuosjH6y7FeyV4kdZZC2m9v49azCq8YR3WcgUVs35Y7soVJUMVeTULk8NLNSn8NF0lmlqGx+suxXsleJTWWQtpvb8PFZjKshZeN+xKtBYzAkgPddjVfJgZv1r4cbYK51ZaX2UoWXqow255ivdK8CitsxbSen8fKjA+kRH8WHTcOpTsOtSQYu7MVOT+Sbmr+TAzfsTw41Ayt6IqteNCU3HdVXeLYYr3IvAorbMW0np/Hy4w1aELD42HLtHKcAolny5yDaWXJkwLz9SkMAEqzA7X0fG6S/FeDF6kddZCWu/vQwVGEarXuvDqb8OekgC1Gv0dgbFWYaYpDTUlTZgWnqlJYQ5cbsJlJQ7WnZykeC8DL9I6ayGt95qoCFOB1yZ1eklSvBcjrfeD0G/G3Q/d4WKkTi9JivdipPV+ELxgHvHy5IkvSYr3YvCarbMW0npD+K2kTkOYn7TeEC5F6jSE+UnrDeFSpE5DmJ+03hAuReo0hPlJ6w3hUqROQ5iftN4QLkXqNIT5SesN4VKkTkOYn7TeEC5F6jSE+bmv9YYQQgjhPK2zFjZbb5NCCLOSOg1hftJ6Q7gUqdMQ5ietN4RLkToNYX7SekO4FKnTEOYnrTeES5E6DWF+0npDuBSp0xDmJ603hEuROg1hftJ6Q7gUqdMQ5ietN4RLkToNYX7SekO4FKnTEOZn3tZ7ZjvWDvnjjz/+9a9/ITS7n0ZRAYH95S9/QWMB/va3v0l4Lrg97rnG8xSe7jB0kE5Negvnt/v3v/+9FEGDDGkT33GmQLbWtiBSkuHFkGNNKoxraWj6ClSKT8kGErrzc1eWv5ThDwEO73X1gFuct8Hb8QHDi/hddUoy4IRaaOMl/iM+zxTI1tqUZHgbwzodl+7Q9OmQBGwET8kGcrrz87NZXrnl++qMDo+pV8Spr8v6Tfdm5rn/q8L7NumVkJ9LmZ6qUzIBD+RkG3+B8tskOVMgW2vRQxt8kZIMr4BHb1Lhx1ovSaDUf1YurmtYntkIAVz2KgCxlYiakg0emnbReznCQVdGQy/k0+HpNqQX1SFDT90Wb5hV/fDnepDesER6y3ZCbJKtlJm4rVmoDv11Jns5DC+CG27Sy3hWnW4tJ3tVX13CWN8VCAJT0sgADq41UgoNWWsn3SrpRa01hp66Ld4wq/qU5GfCtTepMFDB0PQVPCsb1n6ksZKUrYXkb4EqG9moegFZtWSf0h9xJZhSGCqJKsiJZrccSg/Sb5nJZ5V9hAp6bQfciWxkb0EG1YlkQAaUupNqD8hVL2V4EW+7YaWcX/ledpKBbOkyUJCZtWC1tWRnNfIw2YZrO6z3Wgt1lWSmkMGy9CD9lpl8VtnBV9BrOyDy7kIkyKA6kQzIgDIlOSfDax+/xNte6FnZQG51fjoNealqdIKKpUjbT4Wmi8o2XU4fcQXVG/YqGwt1dsvhwX07PUPpKzjBVfUmHM+WE5SskhI6M2O9HYbXUV/kpXQVcS87yeCpzgaZTRHq1l0YS+K35Px2baXq1wHUWTToJYN3lE1KMhyhPpMZqGBo+gqelQ3kYufnlvslyzkR1eIsr1Qz0a0lSA2r/qArqGe0h7Ww49A2O2bS8ykzcMl1VA+21xY7TmQgJWDGfUrGYHHWGNqHV8BtN+nF1Bx+gJ3lzhMLUPOw6qsMuNWw6rfWVmo8tlkLctXR2eyY1UhESvIz4QmaVPjh1ruVi/eCny6rqmenb/15docuKjuv+i1X0gvs0SwJ32KzB2ZVJ9bsxHZ83zZY8BZbYACqam2x46TGUL8OEJiSHmxvIbwOLrlJL6a+/gOwtiaJUSKtE6zmYdV3YTDcL8zO3qQkwzvh5ptUGCfW0PQVPCsb1n6qBkEnkqAE3QGbLqflau3zW1dAGWiVy8MCSvvccWgb2DKTvg02zDobkGfHM3SCYINOf9A+vIh6+S+FjXjcNrgf5cPaA6XhFlgNULp5VD3yt4W5tbaSkgzvpL6IGahgaPp0npUNSq/qp/N8+xG3VLhlynh9Uq119frbQXopxbeuoEayLo86C8hrh8f3RdAXk5ZYb6SXDfjeHA/YQMa6h2oA7K4AZMOnZZl19uFFvOeS6ys/jBKVxGjjJXhXWd1CllsJtlWY+2srVe9EtdCtQpY3kHME2UgphmaAwBBBS6w30ssGFAbKWkE2kHFK8tcxvPnxY7znkdbZwBBlGxym5rpA4xIFZCc3YCxcLRV502e1IbDqU8gGhq6gntFR1fC03KfWEB7YV3uBzrteBa5GoRKtxjZAg3P577zJv2QE2YMPW+2llxyei2771fhZzWNv2uVe54FMk17J4wSTEmNlmvNNmSm+XauhqMeRfRWgW6UheEembGyaUUoyrNATdIxLd2j6ITi5Q63nMCEfVacpTEhJ/kbSeg/xycnNV5t/NPbPy2FaPuqBPrMwU5IXIK33EJz9k/80huMb/dlXmBbeqEkfAIf9zMJcarGRkvyN8HBNKqT1hvBbSZ2GMD9pvSFcitRpCPOT1hvCpUidhjA/ab0hXIrUaQjzM13r/cvX/5zLf4XvDH8s/5pt/58feBbjp2z6Ik7+IwrdKrTxT5C/jXmch5/7PTf8rDrFw11/T4prOVkIz0WXAG38hYOcM+fP3OE7vyfvivPpgb30pMOsGCfKGxKIc/qoCCdrTP/qvPpcM3m7NSfjnKT+edCTb/ohnHnuX1SnqtC7/oIu23E5bfDTbJVVfb4Jc/7kl8nbTnQyzvO89KTDzBmX7qtLuqvD892CW+Pl9v1gM08l73AyzlsGTVD/U31vzsyZ5/5FdcoZ71qrrdmxjX+arbKqzzdhztfwHuBtJzoZ53leetJh5o+L4eECe4yTrddfEDt+uFmmQPWDgL0XLjM3mJITYMmy9Ea1Gb4QbuXZsFyanbWOCjrN/lqGnkJo2kVvGN67HLaW+FpgfdI2US4fmSWe6nYBbVS3Q3nX5W+FVPWgrb2dbPDjLaRfbG+CZfB2UtYplMejlRJBepDe99Nd6UHs5z3ovG1wJxxQt7R1FVCvru41vNL6goB/ZrFEHtrDMDE6hmvbeEEa0T0fAltYaf9b8RgZVDOUx7MLhkfrwqvIRrMIlsFupeTTsrABbjvP7ChNjcc2coIgPUi/E+fWk1l2YDix/62NYCcwOdza8Qz4aVJhoIKh6evwgz0Ga7W83mDH+nItSAZlwNqVZIwXq9vlWDakI8vboGyxsxa5JrGW133lZG0vG5C+A/3DyzVVl0iWK7mV3pFXuTsFaMnwfeUTNLz38vdDst4hyV5KQF/38tT+diD9Y9GC9Cyvb6G1D2CH72H4jgfxkbeuQvp6XUwhb10pYCAZD+j3n0Cyr5qDyH8Fy85+mFRGe9V9PcR5TY9hPEb7goasRV5n18HwfLS6tiI91OE6Qm2HsGXA7vUObS/B+yLXJwbpWSu9/Nu+UvU+13AjydUAkMEbVXtArnor1zuexM4rAxUMTV8EZzuz3fo20Uiu1Mt1xiDUtSiV7oDSU345gY3NTPUPDmNrLZ+KQXi76qdb67Bl47ArJ5d3Sxxth822TsGswwBs1n5k4zCqjf1If1dIWls9C2ZrqMh6oO42trY7Ga2WSwk2q+/1AGfW3gsxnwxVF751Fd3VMdTVdXpfKeCQWd2hnMOWvZ4J/IhrurXIaBC2nqnqMa42LNS+W/FUtPbe7HJ4W0fbD7tuN4zQ/rcMOv/IaNqg4OV8Vnvrt+LcOtdWYAwdeXXYmRnr7XBrx5PUYMxABUPTV8CBT+7l6xZbt7xkSzOzbEEQiVPH2eDHqNRVQmZ6sPVDVrQWYZimS0Q3g+pQ2K1t1nRb37V8J1pgeVMtMESJsH8KMTTbsTl4+cOQoK7yDdTtZKCptb5Ds9UM0N8V7dbyTn8v+GnSiyFzzuylK1pfOOgq6osIduRBd65UKLDh0xvb19m6l9Bs1eNWJWOhYzlKc15lwJXP1VHNxHCtZF0dwk54UDeyTefWVP1OhDLbMdCUtuuuiKEsBUOUcigDQK9jdvpK3V0bgezr7oIA9jeCbwMb7ngSXDWpMC6noenTYZdhQh9HedlRL934Vapcld1DoldsTv1vwQzj6mdrbbdXxSGt1zo/atgdtnlg+c5J0es2BMPupB11iy2zHRv02u7ekNrgC5Qgt3W76nZL33Ey2q3lVf8AW9s9F3bRAR+GMzrU4VWsr44hyp0rFcxWg2/thVb5BcV6rUvGQkc9y8FzDRmu1RC9Ln8nvAo2Xl7dVqp+J0L0zO4YgGy6mJEVs5ANwtYxt+KsYNzZrwNjuL8RwjAwCxWUoIUnwU+TCpuX3qSXwRbDbL6LtRNuvF6uwVKX272EbxyBoWRYsvrmWfojD6D39kLYWltj6HBIXTxaoikLa4jh4eVaso6qc+XhkVNAt9zs2PgOpT8eUqcHhiirADh3klS9lt97om+j7S6qLq9uH6CG8SLYQqc7Ax504VtX0V2pb9JCm/hfmJJnf89u2Xf+gSHKNljobBSqbFxWHdVJle1Kwlb8ZrhW3HKrZJeU4PA6PdibhY6q1/J1hPLP55aB0OU7SJC9ZPDQDjs9IDCUbDo/ILMusDZ3YKNv7TsDYIiyDU7QuRUDFQxNnwjv5O++h+HJdWVtvLC+PqiPsb7xxaTv2bWeMXNiaVPJHTpU9Qlba7HskpXPLiRkx4A9IMhGyo6Ty4GpdbTVrWTpAZ/fngID25jOhr0UnkB25Jh5+bchSfZahYRSej5Ryol8Sr/YNjwFw+2kvytahSElYCCbzu0DVLevgDjrSR/GV7R1FYBeNroWmyHYxlcK+FRsnc+hvXzKPzgxNDTVxv63nqnqO5vlWP891zqeSrcWG+0rkB0SZlvhWT/M+Yr0bbDAcB1hvdWhgUBPGGi8Ud1XMiBvPX21r0i/PlcXmAxkrAs5spFkmdleyvWOGp5B/jsGKhiaPgsfu6IpBA4seZ/umoz06/u67bE4Z8n+CwkZt8HXcqg2HXjGYB2SFkK3tmkXmqrEWYfgtUwp/jXdiUBr4chy0Rb8b7Q6muh2kVJI4xsWyOs76ZywV92OKQza4M6QuuxyJrBWGoWnkIa3ITPwdiej3bkQWWq5heOwpEkvoLtJoSmE46HKj+x3rgKDZYfbF2hnJj34SjWUDBh3LyKqcisxOtr09utXZInBzrlAZlDjMZ1/bKoZU7o6sbi5UW22jqZhXQ4Ma6hClmC33z6BwAxlPSxIKXy6zmG9IltqaIbnqn5sgOZ2a0tsOxshyB6Gge3c5Dq847C8SYXNlGpSCGFWUqchiNpiZyOtN4RLkToNHwu/4/q3Vf8KOydpvSFcitRp+GTIf+M/HJ4QwmtSIa03hN9K6jSE+UnrDeFSpE5DmJ+03hAuReo0hPlJ6w3hUqROQ5if6Vqv/5VV/WdY9yIPa3D+g3/trf5bsXeif9kGLzr7vW637mH5N3j/8w8EwwO8J72fVaf4aYMvlE4H/4LMT9XUFsuV3GjjL55Sg6847F0V9/QAZnu+dzJMg3FmPJwxx+ElvAtPcv6LGA/d06I57/YXUZP7RWe/1+2WPXESbRuER/lFdTpc+3vTYKutvqEGH6MG9iNMdRtvZpgq49J9Q0mzBekr+cyPh2Zd27+3sB+D4/u8Lzr7vW6H9vpPxvj1w8P8ojolOUmGNlh4StX/FLc2MmokKF9dg49RA/sRprqNNzPM83Hqv7kknlKEeOieFg2eeXIE8J9r6atfDLNhKau+i0ujUEW1YWi3ErTd0F4GCNKD9MIBQ1MdiNmrtAvC8OySrbQM1W0NGzCTHvn4lUopm2og57IJZ3jzNZ55OLKiW8tQqVKTbaqagqFNGy9IIw7W4JF9mcJMlhoKzQKyrhTYQkrZ86khBprqAqvUywSvtSyfgBPJVspM3NYs7LyOHH4gnL1JhYEKhqavoz7SY6yfVhoryRttIT35sVjdtrZsnLJCyWTBuyCrcryX9N/aSw/Ss5f0YBsgYIVxJGbZaC/Ho2F3drgt2HWLrLXQheEpu9Vx6tp6TAsykCyH4SS63rfBdnrxB1DyOKnIEAWvfLAeWcnjZJO+y6W1vfQgPQnmOrIN3FVT8lltqk/HYORTeskeulikP7KvBcngQ8mJiwi5XoKUwI51U0XSUfUYy2cXgAwkVwNABgdW7QF5GNinMTz7+Dredk1sBDxMGz8KD9zF3GlIDmWhs0Qsydl/m3SJshXhVsKhVIJ22H4dm/REUhfiUG5tIL6Nee1f9l2cW261b7UUW247P8hoqtAZeGE4SX2Ol7Kkw9k6rR62vDlVulx1LnXYnk/spQTrybS6ULnNpw3ELfVXOdnZIMuVYpOyUvXreOT/yL7eaO1Ea7sA7MQLBbIueSvgnUqXH+8oGHqj6rAzM9bb4WdS78oMVDA0fRFKC6XIw5AQ3dMu2fjfnNYWTrVKNROuT+QuaZRzRonV7VWHR+zR6/gWKgdjrj6rDNjLbdXvu62zugeoy4EpX6ltwDcm+7WBCzKchItt0uvhrfXibXw/fndniGCIZyMbJY8MoA6P2KNXqBYqysmOuhx2EttCR42hyoArF0tHNRNea0FgrLN0euLRsOpr/FXfUUPyYWVfPQi/oAykBMx8yRgszhpD+0+De2hS4edbL/hFHwYP3dNWn84h8uPg0TBzqTj5kGvJyQahi9/DI/bewoL05mDM1efQf6c/6Bab4XK7XftxmaFndm0gfRuEE3QX+2pqAjxAlxhSIuNWMniq28vDI/ZOTgvSm3VOrlnbOH4LHTWGYTxH9gWvHTrp9MBQ8VR93auzH4Jx538dLUOUCNVhDQyBKenB9hY+E47fpMI4D4amz6I+lTiSGfusn7ZqEHQiCXXrLRRSrTE7EXVY9/Lptuy742/pjQzW+g5sHEOV7R/W+q3t2mDBq9bLqyB0CpQ+Tmeg4bfHCUeoF/t0/IJtfLpOWSsP/nYe5oZkBIaS93MJoQt1S29ksNZXur3kCiUy8UvosEEn25WE/X0dswXp7QQQ/BAy014WgC8u33PVm+pQyMz7dgYebgV20P4DqddiBioYmj4RcsItDeHkdnrj+rTdY9/659d26C2TvltbY0OQ1Un1KRksLyY3hglX7bsEXUJr8dRrkRmfyAjWD2Ne7yUZFve3tdJLKYZuZebC9t1Kz6f0dgvVniMAQj2mDeTE+nCSV99kTUiEk9spJUBZDTWpam5IXkxuOJe27GuywZKbLex6BJlpdwTrhzUFKNeJXWOorGOzzRLOf4tlf1+fZetQcq5gAKHqtak8r/UV6X1AbYey7msDGWvT/cC0kWSZdfYfyPD44xt5w03xireX2fgDirtYZzCa6hbZGQbLtjeUMUPks64CKYXzCaHuVbce2lcD2IoNmmqhqbZj1iy7bPnv4hRaBdWty0aolg6GbT/VXmcHrbIeDVOSwwNwgU16GTzW8nRPqFPlVfUD5MPi/obyBGWXqzVnhvbVAJDRtEFJTmiqhaba/R5oFsXGm66RJQYH4xnu67VbTvhkrT47J8hSaq3sQUoC09B8W+k2QINz7VUNoG6EIHsYvo70kj8H3UPHZg41KYQwK6nTD8Qt8M3UFhvuIq03hEuROv1A3tYCafD+bXXnd/3wLWm9IVyK1OkHwqO/7c9s2cvoj6PDA3B7TSqk9YbwW0mdhjA/ab0hXIrUaQjzk9YbwqVInYYwP2m9IVyK1GkI8zNv6z3z1+WH/1hQoPRB6r8tG/KtwY/w8LUcYc4jh+Ok9YYwP5O2XhoA2z381+X/9fV/xGzjL/S34X91a0lrDPu8s05DCI8xrNMfbr1qkPDwX5fXb8zrgGla8NLfGl8NwT/8E0n4BNZpH0KYjWGd/mTr1Z8Vq/s+/I/G6K/67bB6QEPf0pQ02gJ0ND6FZgHZBvoUTPnnAxzK2GYa+tfTg8s7bCDkFmMN9dND9Vn7sYbS3+UfJEtZkcHWjmEeeJomhRBmZVin49J9T0mzC10BtnrGEeSE/sRnU30p9clQsgXJwL5qKjsGyLX5qS3ZXmAgm4PLO6oeY91GtZeM88Wknc563x7ysEfaD9h/dwSQtxrtescwFd0LhhAmZFin49J9Q0nTANQn+NTX/QOoQ/CJB3cdOaytBY36DUI9Wo2hGrhRoXRs3gvZ9gJZneng8or066kaP34Up8Cn3FYbsL6y5b87AjC0BmG4Y5iK+vohhDkZ1um4dF9d0vWrvPuWvwv3DwQ7VPC1tXg7CwJLtUzrLQgbQO1z1Uy9TY2t6mFreYe7I7hB2lWdNZqyjaiXUNn3L5CZlbyzY5gK3qVJIYRZGdbpuHRfWtJP/Ga/dY9lIY3NPRgQGNqnu7sFUBhqRd8aAHpmJVez2lMPLt8Cey+xq52eXbcDhj7ykKF/QLAednYMU5FnCmF+hnU6Lt23lXTXn+7F/QMPyAiOHMFT2sKCDNxgqh4B/TLfd6Cly7fGVs1Qams4uNx0NmAPFmTjsCvo3T51Ci0xW/7rkWVTF0oz3DFMRfe4IYQJGdbpuHTfVtLr3nAXrHXPQKYPabhuLVUQ7oXW11WAN/dUqG3eZmiQ5efg8oqWWK9IUErPp/TeArQjgmy8BYJtjGzW/rsjrwNDud4xzEbeJYT5GdbpuHTfVtJ8v9eGwb50hTb4DrUNPjWk97gPubUAnUN6CwJZLcf6ugp2YkPPELRWfo4vr+gUxsfR0Es0BPtkU4b6rPqOoX8fmVVt4osuADhyivAj8BxNCiHMyrBOx6Wbkp4fOmJtiuEDSZ2GMD9pvZfCv22HjyV1GsL8pPVeCt4of/b74aROQ5iftN4QLkXqNIT5SesN4VKkTkOYn7TeEC5F6jSE+UnrDeFSpE5DmJ+03hAuReo0hPlJ6w3hUqROQ5iftN4QLkXqNIT5SesN4VKkTkOYn/tabwghhBDO0zprIT81hxBCCG8lrTeEEEJ4K2m9IYQQwltJ6w0hhBDeSlpvCCGE8FbSekMIIYQ38p///D+/iqx5ZqKeEAAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 01 - Preprocessing the Signal\n",
    "\n",
    "<p style='text-align: justify;'><i><br>Time-series data coming from dynamic simulations provide useful information that can be used to determine some aspects of the system. Before estimating the modes, we must first capture the signals that have the highest oscillation content, which is done by calculating the energy of the signals. We have the active power transfer as data, therefore we can determine the energy by simply taking the integral of the active power transfer per line.\n",
    "The equation is given by:<br>\n",
    "\n",
    "\\begin{equation}\n",
    "\\label{Preprocessing}\n",
    "    \\hat{y}(t)=y(t)-y_{\\operatorname{mean}}(t), \\quad E(t)=\\int_{-\\infty}^{\\infty}|\\hat{y}(t)|^{2} d t\n",
    "\\end{equation}\n",
    "    \n",
    "This step is done on the “test_filter” function below using the variance formulation. Variance measures how far a set of values are far from their average and is expressed as:<br>\n",
    "    \n",
    "\\begin{equation}    \n",
    "\\operatorname{Var}(x)=\\frac{1}{N} \\sum_{i=1}^{N}\\left(x_{i}-\\mu\\right)^{2}\n",
    "\\end{equation}\n",
    "    \n",
    "<p style='text-align: justify;'>where $N$ is the total number of samples, $x_i$ is the i-th element and $\\mu$  is the average of the total number of samples. Comparing equations above, you notice that both express the same concept with the exception of $(1/N)$ term present in the variance formulation, which in our case does not impact the outcome of the simulation.\n",
    "“test_filter” fucntion , which can be found below, sorts the input signal according to the variance. The signals with a threshold lower than a specific value are removed. \n",
    "<br><br>The script has the following configuration:</p>\n",
    "    \n",
    "![table1.png](attachment:table1.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Empirical Mode Decomposition\n",
    "\n",
    "<p style='text-align: justify;'> <i><br> The empirical Mode Decomposition is a numerical method used in signal processing to determine intrinsic modes that combined can recreate the original signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_filter(sig,t,var_thr,step_min):\n",
    "    \n",
    "    '''\n",
    "    DATA_FILTER\n",
    "    DESCRIPTION:\n",
    "    Filters the input signal based on a pre-determined variance. The data points with a threshold lower than var_thr are removed. \n",
    "    The time vector T is required, which has variable step size. If no oscillation is detected, the function gives empty matrices.\n",
    "    INPUTS:\n",
    "    -      `sig`: The signal matrix based on the number of buses in the system in study.\n",
    "    -        `t`: The time array.\n",
    "    -  `var_thr`: Variance threshold.\n",
    "    - `step_min`: Minimum varable step for signal analysis.\n",
    "    \n",
    "    OUTPUTS:\n",
    "    -    `yh`: Output signal after filter.\n",
    "    -   `ydh`: Output signal concatinated with the time array.\n",
    "    - `var_s`: Variance of the signal in ascending order.\n",
    "    \n",
    "    LAST MODIFICATION DATE:\n",
    "    12/13/2021 BY FERNANDO FACHINI\n",
    "   \n",
    "    '''\n",
    "    \n",
    "    from PyEMD import EMD\n",
    "    import numpy as np\n",
    "    import matplotlib.pyplot as plt\n",
    "    import pandas as pd\n",
    "    \n",
    "    \n",
    "    nL = len(sig[0])\n",
    "    ts = t\n",
    "    trshold = var_thr\n",
    "    tsmin = step_min\n",
    "    \n",
    "    double_dif = np.diff(np.diff(sig[:,0]),axis = 0) # The double_dif array is the double sequenced difference of the elements\n",
    "    std_dd = np.std(double_dif)\n",
    "    \n",
    "    if double_dif[0] > std_dd:\n",
    "        double_dif = np.delete(double_dif, [0,1])\n",
    "    elif double_dif[-1] > std_dd:\n",
    "        double_dif = np.delete(double_dif, [-1])\n",
    "    \n",
    "    loc_peak = np.where(abs(double_dif) > std_dd) # Position of the peaks in the sequenced difference array\n",
    "    n_pk = len(loc_peak[0])\n",
    "    print(\"The location of the peaks are {} and the length of the array is {}.\".format(loc_peak[0], len(loc_peak[0])))\n",
    "    print(\"The last element of the loc_peak array is: \", loc_peak[0][-1])\n",
    "    \n",
    "    if (len(loc_peak[0]) % 2) == 0:\n",
    "        new_time = ts[loc_peak[0][-1]+2::]\n",
    "        new_signal = sig[loc_peak[0][-1]+2::]    \n",
    "    elif n_pk == 0:\n",
    "        new_time = ts\n",
    "        new_signal = sig\n",
    "        \n",
    "    \n",
    "    dt0 = np.diff(new_time, axis = 0)\n",
    "    gg=dt0.size\n",
    "    timex = pd.DataFrame(); signalx = pd.DataFrame(); dt0x = pd.DataFrame()\n",
    "    \n",
    "    flag_EDM = input('Would you like to run the filtering process with EDM? [Y/N]:')\n",
    "    \n",
    "    j0 = 0\n",
    "    for k in range(gg):\n",
    "        if dt0[k] >= tsmin:\n",
    "            dt0x.loc[j0,0] = dt0[k]\n",
    "            timex.loc[j0,0] = new_time[k]\n",
    "            signalx[k] = (new_signal[k,:])\n",
    "            j0 = j0 + 1\n",
    "            \n",
    "    dt0x = dt0x.to_numpy(); timex = timex.to_numpy(); signalx = np.transpose(signalx.to_numpy())  # Converting to numpy array\n",
    "           \n",
    "    signalx0 = pd.DataFrame(); osc_signal=pd.DataFrame(); var_s = pd.DataFrame()\n",
    "    \n",
    "    for k in range(nL):\n",
    "            if flag_EDM == 'N' or flag_EDM == 'n':\n",
    "                signalx0[k] = (signalx[:,k] - np.mean(signalx[:,k])).tolist()\n",
    "            elif flag_EDM == 'Y':\n",
    "                emd = EMD()\n",
    "                emd.emd(signalx[:,k])\n",
    "                imfs, res = emd.get_imfs_and_residue()\n",
    "                signalx0[k] = np.transpose(imfs[0,:])          # Detrended signal\n",
    "  \n",
    "     \n",
    "    signalx0 = signalx0.to_numpy()\n",
    "    ydh = np.concatenate((timex, signalx0), axis=1)            # Array containing the time component and the signal\n",
    "    sig_var = np.transpose(np.var(signalx0,0))\n",
    "    sig_std = np.transpose(np.std(signalx0,0))\n",
    "    \n",
    "    var_sort = np.sort(sig_var)\n",
    "    var_sort_desc = var_sort[::-1]\n",
    "    s_idx = np.argsort(sig_var)\n",
    "    s_idx_desc = s_idx[::-1] \n",
    "    \n",
    "    std_sort = np.sort(sig_std)\n",
    "    std_sort_desc = std_sort[::-1]\n",
    "    std_idx = np.argsort(sig_std)\n",
    "    std_idx_desc = std_idx[::-1] \n",
    "    \n",
    "    nL0 = len(var_sort_desc)\n",
    "    vA = np.max(var_sort_desc)\n",
    "    var_sort0 = var_sort_desc/vA\n",
    "    \n",
    "    for i in range(nL):\n",
    "        if var_sort0[i] >= trshold:\n",
    "            osc_signal.loc[i,0] = s_idx_desc[i]\n",
    "            var_s.loc[i,0] = var_sort_desc[i]\n",
    "   \n",
    "    osc_signal = osc_signal.to_numpy()\n",
    "    var_s = var_s.to_numpy()\n",
    "    \n",
    "    yh = sig[:,np.transpose(osc_signal.astype(int))].reshape(len(sig),len(np.transpose(osc_signal.astype(int))[0]))\n",
    "    \n",
    "    return yh, osc_signal, var_s, ydh"
   ]
  },
  {
   "attachments": {
    "table2.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoYAAAEJCAIAAACG5nG+AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAD1tSURBVHhe7Z3NmTS5rayPJTJgfNB6LJjV8UQbuSEb5MB1RP6cG53BDqHBn2JmZVazKuNd5AOAIMgkCaLqU9fof/7PGGOMMQvgkmyMMcYsgUuyMcYYswQuycYYY8wSuCQbY4wxS+CSbIwxxiyBS7IxxhizBC7JxhhjzBK4JBtjjDFL4JJsjDHGLMFUSf4fY4wxxpxBqawtZkvy/zPGrI3z1Jj1cUk25hY4T41ZH5dkY26B89SY9XFJNuYWOE+NWR+XZGNugfPUmPVxSTbmFjhPjVkfl2RjboHz1Jj1cUk25hY4T41ZH5dkY26B89SY9XFJNuYWOE+NWR+X5FuAPfI23RwfgDVhbpK//e1vxWruCo5BqawtXJI/gX/9619XZ7vPwPp4j1bjf//3f7EpNaV5gl3OPU4JYs4C21EqawuX5E/gH//4B/bo73//e9FP5dLg5kScp0vBxAEozMW07RGYyaZT8s7JuyDYkVJZW7gkfwJIOewRMl9fl3UdgOL0fR2oKSYqLZTjd2482QT8b26Lgz0qklkAZk2sx4DJBahGeZx3lJ28HwC2o1TWFi7JnwBzD+mqjI3wUtBdEGETeylpqTLnN6+CP2svDvaoSOa3STkV2ZLpa6fm887J+0lgO0plbeGS/AlsSfc/yFtkKWVkJuxUY+oCuKmJiRrdAIxQGQF89fEBeAe8TeuQckqouEKezzsn7yeBHSmVtYVL8tsTk5wJ2czPmMaASc4P1ymNYYTK5I9uZnGwU0Uyv03KRBHTcD7vnLyfBDalVNYWLslvDxOP6RoTMqkpdWOSQwC0g6imu8CsjHbN/Dq9xPlKre80pEw7iGrq7uT9JLAppbK2cEl+e5h4eEKGoM1CJY4qZeUnVSR5cmPy65N1ugvMymgTza+DlPlKqpA7SjQm16682xydvB8CNqVU1hYuyW+PEo/5qYSERapSN9JsghFPJT9V2mkxy4JtKpJZACZmQnk0n3dO3g8DO1IqawuX5LeHiYfkZA1WQsZvz2qiAGKW6u6AM3vBjU1U2USLWRZsU5HMGijdiNKKTOYdgzh5PwbsSKmsLVySbwGT05n52ThPPxIn74fhkmzKN+n0Id18GM7Tj8TJ+2G4JJuvQ+BN/Hi8xR+Jk/fDwG6WytrCJdmYD8F5asz6uCQbcwucp8asj0uyMbfAeWrM+rgkG3MLnKfGrI9LsjG3wHlqzPq4JBtzC5ynxqyPS7Ixt8B5asz6uCQbcwucp8aszzkl2RhjjDHPUypri9mSXCRjzKo4T41ZH5dkY26B89SY9XFJNuYWOE+NWR+XZGNugfPUmPVxSTbmFjhPjVkfl2RjboHz1Jj1cUk25hY4T41ZH5dkY26B89SY9XFJNuYWOE+NWR+XZGNugfPUmPV5p5KMgf7zn/8UZSd//PEHutf8+eefeB4OeymY1fMTw4v/+9//Lso3sOCti7IAp7zpJL133zUHeC57bHq8bMefXJlmqv7zn/9k61kr/7IUaL7OeP6nvCODgKJ/AwumVJQO6NscnTGbTeYs6i2LTB3ZcYiz+Ouvv045De9yqk6ZJy4yrFtRAkuV5BfvSPPd987hxXM+hdfs+PN5ioKhAky42ti4op/BK0tyeh2o46FPOV0MkuLQOC7Jp4xuDoPFL1KLqSM7DnEKGII8f1De5cCdMs9eBJfkonyzdw4vnvMpvGDHMQR5ZmXqGgZQ6ZsfLg/ziyUZYOjBJ4xTTheDYNHi6JBhcUlemfGxnDqy4xDPg7OLY3TWQUlxpFJgohIcXxpJTCGc6WJtvXtKQshy63WUnfdOHFehkg+gG+MDGgXs6FKUDbqBOCVQRybN2ca5gfpmqVeyuXSDNyWTcUAxbfTsxVS9O2jOobcmhF0YCsC5NGz05vC7XD0Z7A7WiiuDZ7HuB6uXTgJAZG6E4lPoHY9i+nYo1j0pEKfBIIrPyYwnIJqvE53roRkZT6pyALQA+oh6aE0P3YvpezKypCBoihapEGihiictEBgHMiyUzZNgMYvUYiqHxyHOQqeh6EdJcaRS0GHlmUsqZRiVPLDLRyQjZJ7XXsdoxyjMLk2M9ugjmT6yJ+DGcQk8pUIGlJuRk6zZpllRTdcBjfQH6DseK8UUk3E0N9CzoyNUyYoj0hxinCiLNLeef5zDr1O/9RWklTwAVkybRXDAEDPlBQUtL7roBSE0t7tnb25f3DtY4qwgYzKDCURiR4Jo46EZGc9oBJpSdJDaTEMY2YUWyNESe2n+MThlTSA2wUi7BHMK3IUeozYxDnEW8TQ8Q4ojlYIOaFPFExYINJLoRuQc5V7H6ByJ9tRXKn3S6CI2pQhKv17kZAdQYeSIxdQhzUov0htLDrSLyTiAbpP25tUZ55D8k0r2zq0ov0r9FlcQV/IYKBiIkKiXmsLDLegd9YcpwIA0YkrwZyWTvTcBqqL5OmwaD41ncgBQYaRDMXWgD5wxczxhkawX2RwLGlSCjOwefWITnmw1pzBez6m1fs2WxNPwDCmOVAnJnlRmcgJGugnmMAQ8mQC9jr0zHSegLCJqklAafoIm5ZKmQdgRQi9yb7bwwSdiqr2PxmlWUve+RbJL7c1tYG++eySO1Ztn0TfqIFDRsTeH4vSrYCZFupLmcu0Ci68VY7S4gIovIdl7292zD7YbdrRCpkOKEz1rVcA5zj/yMCPQEUKC0SbTEEMA+uBJNQ4KlXHINvh/XyTKtcrpIQJVcwpY0iK1mMrhcYizSKfhMCmOVAnJnlQl5BideyVkryNToiiBOIGURWqSUBp+giZlSxqdHSH0Is+8JjxJ0b9Js5K69y2SXWpvbpN2xinKN3Gs3jyLvlEHgYqOM+v2W9RvfQXN5doFFhDLWJTvBMGTquJLSPbedvfsg+1GDWMZwxMq3OjM6UXPWhXyr+kNLSHNuQk8SdG/YRCuG4PQJw669Stz1qASonFzzyqWBSoXx5wFlrRILUZtYhziLNJpOEyKI1VCsicVBzraB8CNzlR7HRW56N9Ee+qrsL2+BE3wpBxnAqT2Iif7ALjBuSgbaVZSe2Mlf5HsUntzG9uLUqlEwSGnOA/9o9qbwwrUb3EFaWUOUNcwXv2UFV9Csve2e2xXnOgGGZNhVYaKWXEmdNaIm2/3xevXEb2hFSo5DIAbJykYhEZMm5OHDAtLsgSisdgRTxijnFT6Q8CTo5hT4Kr2GLWJcYizSIfjMCmOVAnJXqs4xzrKOsdUI8xepgHpdYRRblGGg8560yfOqgY+8S6Ik4EMKPdGhwwoa7YUaATNCSRjVFP8OB+9qRjHAbSnlWzaITTfPQKj5oAgzXkKTkZjJX/Z4xx+neZbnw5X5plXxuqlGsaYXGHFlxB9qEJQBMhAcvMYYMTedtONYbmb2tzBBCL160SaQ8dQMIKt/Uga8khjApA5DVgYUNHkLBUCO9JOY1LlwzhbuzmB8WJOLfRr9iMdDoCDNTjrPVIcqRKSvaliaKhExgRPKk+t6HWUPV4HkKOl9kmzSmBodCnKBrsDpmixdkYHzdnyvUQ9eppVUmfelEzGATKCnr2YqncXaQ69NSGcDEPVPr05/C6YTJGuhCuTdmRXnjb9udR4Kr4EOiQVMmFHGgGNINl72w0VTUXZumtu4wmIh69fD51CyQHIOJmGvH9iwHgtYMSt9xfRn3Y8Y0cgNS0LZE4ewvhlzUOwwkVqMZXD4xDmF4npZG7OPfM0ViBj1scl+ZPBJ1Z+ejXmJnnKr4/6JKovcMa8BS7JHw6uJP5jlLk598lT/qM0cT027wUObZFauCQb8yE4T41ZH5dkY26B89SY9XFJNuYWOE+NWR+XZGNugfPUmPVxST4H/WJPQmlYHv6FKuX3mrnZxQvyFEOIYtrPH+FnuKL3a+DDxGN/KePX2Uvv9U95HQav48Dy8Idk6FvPCjBms8k0qdc/MrXH4xA3QSdPQmlYnpjMkB/mnnlTrs5TxNefN0M4fJDQMVUs5hQOZ9HP4JUl+V1ehxMD8fqicbyb9Im9zGGwkkVqMbXH4xA3QYdSQmlYnpTMyL1zLwuzCJfmKU+Rjv0zWVDXMIAar3p/CqfUsBne6HW4a5hYnDBkWFySX8Z4H6f2+DUn+xRwsDBbUkwbxbR9dMUTZysdsqhSFjy+cpCwJd1/s44nuyjf0FmwFtIIfxpjMiT/mDl6tTTKw1fmQMW65fY4/cybEnf5apQFRd8Djl882GRLpq+DrcgUmLCE6UOK6duhWPvHvk6fOA0GUXxOZjwBMX4dotFBMQ0vBzzpwyag10kOSW0OJOiMgeBWTN/zl4U+Ak3RIpXzAVTxpAUC40CGhbKJYGWK1GLUJsYh1gGnSmmA06BDhvnLDhng3AAKtEcVgnKP5wyCHCQwV+kGMFzKWHlGFT4UAO1x2nSgrKFBz2fwysoHyIAySLMyH0Pc5avBqdNh2ws66nASphJPtc4nBY0S0wFC83j37DFNJMNTwfk66gtZeSqfOIFI7Eji6wA41EnK4HjSTlWD0g5BkSEDCM2OVJsDReiMUdQEC+RooQNlvXIchbIGik0w0i7B1HBJe4zaxDjEIuAYpXlChTHZqeIAxZMEpMpCor3uCAEB5bb1+C9NI6CdHaXiCWgh0U6hNHyTXg0wbLIrrwTdimI+hbTL18ETVR/ISXD7o3tCB1KnnUJt7x3vSbtUBqQRU4I/y5LsvQlQFePXSaMDtsbRIxolddTrpGn0/AFUTYPQGUZEY5NkvfvmWFBwCTIqcrMJT7aamvHiTC3cW6wvj2wCRsDTRnSAJCQ7VZ4q8eXd6ojg/DAoIQHjFuDHvzbHCASqjng9dO+IY1D6RL5euPXKRdlAK3yKYj6FtMsXgZODgXRcDxCPHw9nPI3KDgnJ3jvePTumWtvxhAw7WiHTIcWJnrUq4Dx4HciwJOgwvhw0DUK7BDyjnf4QEnEmgM54ZcBB8aQax4LK7gS92BFPtEa5VjkNRKBqarA+RWoxahPjEIuQTrBonmw8JSQ7ZAiSZa8FtULAEONTCDcCOUYgUNl9c8lDM0M2xx/seuWibKAVPkUxn0LzkJwLq0g8ugdIx4/HW+mjYy8h2XvHu2dH5NqOJ2S8DmsSnlDhRmdOL3rWqpA/Sa+TZtUE/gSyRum9TppGz78JnTk3OjMm3/rLo3MFSYjGzT2rWEmoXE/TBOtTpBajNjEOsQjMBJ0MQXtRghuI/lLj6QS1vwQ6wBnJAAvVMXBDQEaAQCNVCA+HLg3fyKHo39BelEoFUGEsivkU0i6fzvw5H8OUKcoG73HKOu0Skr13vMd2xYlukDEZVmWomBVnQmeNuPlmVYxfJ40+AG5w1ihxnkBqmkbyl70JneEJGZPk+0KGhTePBKKY7IgnjFFOKv0h4MlRTA2XqMeoTYxDrAMOk86TDhPtPHkARtkhKJfgQ3vsiOeX96ZSjgI7IgJUxY/ogBJ1pKAmTa83dPSpZUA5doeQXpkyYGS6mU8i7vLpnHhscGJTDWNwnlgNJCH6UIWgCJCB5Oaxx4jN9AF0Y1hmkBJqMIHI+HUAHBRTSTq+HPCEMQaBDCTXFxfleiCqhMHRBJkXF+PAwo6xF52lQmBH2mlMqnwYZ2s3mfHKTK3aGy0uDhZmS3RogOw8iGziYSLxLCINaAT0QSsFPCVsgX+c8hqGFezFLpwJiBdEc2g26RWiP+i9cjF9D1SsP/8lDQJz0nwAcZdPR8c1wsO59xQ1/RkfTx57PCXQIamQCTvSCGgEya40SekDVekA4KC5jScgxq9DVaMDRRhcDvJhE4ivQx/CIPJvDiTYkbtGmT6waBF6VxDteMaOQGpaSciwUKjX585guYrUYiqHxyHei3SenocBizLH6XPYBTKEOWY+jE/K03liOTFmfcZ56pL8LPgAuPcz4C+WZN9fH8xNSjLOcEwffRsz5i0Y56lL8nEY6kCFO3EOe/mtcc0L+KQ8HYMPwXhZ4nps3gsc2iK1uF1JNuZTcZ4asz4uycbcAuepMevjkmzMLXCeGrM+LsnG3ALnqTHr8x4lGUOIYtrPH+E3eUJ/Dn3WH1XxDz6LsjB40+bLwnjKOpjVeMGxVIrt/ZVB5CZ5qr9Bw7tEuTRPgy4Hev0ibzfhF4NjUKQWU0d2HOJ5EF9/Npl+b74LdEw3BU4GgiM5i34Gb1GS+eJOjFtx9bFEbipPMdbhqnyTPI1vdPjt3i6RffM8ZHwsp47sOMSTMHO0hc/saJ3qIN4jp/DrqT6DE+OGXHos04lCop340Rl8WJ7G5UpLt4tn+v4Kbzfh1zM+llNHdhziXJ7Z0XGqKzIFZiyJH2CL6duhWIMdQ+BZrNugtOtCidNgEMXnZMYTILUPYtJI1CUaQXKTymkDqnjSAoFxIMNC2bwj2MEiXQ+OykUlWeeTQi9NiunboVhfnqckBWdH8ueffxap1QpiTMUBUKNnWjE28R0BjQfekTKJzoDBKaeJiTRic8L1rJqTvw/jV55ajleuGrYNW1iUnaQjBeIp5DnAk4JGiScPgiJAfmhHEJ0zyfBUcL6O+kLGZAYTEE2fZhcISjPZ2R1PySkf2AQj7RLM+8Ktfw0YS6duLzEjyPvmKVDAKLMvngNZKl88xtHckr+gXf6gOY2H70iZRGcgZwj1xEC06y3ShKOP5HrytwLvXqQWozYxDnEi2O+4nXvBlqN7QsdOZ4VCbYcFAo2Ak4EwaZfKgDRiSjrEsvcmQJU0feougBYS7YpJuTmcmvBkq3lfXrOJzDIe6WN8Up6Og9O5lr9cf5LiAKgwxr4R2jW38TRobL5jJI4l+djEaB/PCirtdyOuSc2oTYxDnAVT6Jl94pmjzF2XCnRWJCS7ziuhHULPjqnWdjwhw45WyHRIcaJnrZKxT1J51gXs0SHKtYpZQb1tenwS2MciXQ+PXFF2gixQYvI0SgU6nxKSvZePPfuledoLLiEa6YNvh1BB/JqIQWmMwJj6imQ//I4JGNEEQQ69iSEmhK3TD+LQM4tzQ5rrJkZtYhziFHhMn9wknSfCQ4MnVZ0DCcmuU0hoh9Cz9w4cZLwOgAOeUOFGZ04vetYqGftEFYJk2SUk51rFDKFynuatwT4W6SVgOCXXLpQIBEFiKJ1PCcm+VJ72gkuIRvoIGAnkNHnR65vsh98xoThy6E0Mnpx5Ig7dm5WE0nAzmusmRm1iHOJ5sOunDFGfM5xCRdY5kJDs6ZBJHdsVJ7pBxmSYCVAxK86Ezhpx880qGftI5UA0Ak1JDjBGOamaM56cqnlfuJUXoaMiDp+ZT8rTXvDo3Owo0IReKY7o9U321J2q5ME71qApdW86pwmIaE99FbbX9yZwEXqM2sQ4xJOcuD11qjM4jqBkPCVEH6oQFAEykMwglGXnWa9lQDeG5VmEA5sGExBjH6nx0NMoFQJao/NXz5+qfJQt5n25egdxgHlaAASd572g48fkKWgGj85RTokWm9AX0E432iFo2UXsSDR0kgE85ZzesQYd4RC7f02rNTEYmyPCQRNu+tSTvxV49yK1mMrhcYgnQWohfoI7ii1MqTum6c/4eOocSKBDUiETdqQR0AiSHYPSHg8xgIqmomzdNbfxBMjYJ6oYCDKhnatHO569vmmGkGGhoKmaNwLbWqTL0GlPJ2fXgWn6M63w1PmUQIekQibsSCOgESS7Zs5DLlIWwEFzG08gUgePzqkjq5qIARUHyI6YUNO0U0xy4B1rOD3VVNKcGGiOmCZc+6TJw2Ewn88D716kFlM5PA7xqeBQxnNszOI4T41ZH5fkWfjxMH52ix/9jFkc56kx6+OSvAP+Yxdxnpv3Aoe2SJ+O89S8Lzi0RWrhkmzMh+A8NWZ9XJKNuQXOU2PWxyXZmFvgPDVmfVySjbkFzlNj1uc9SvIf379de/4HavyDTAG1NDyCHYtyHvpTlPh7PgCVFjzr1uvYhj153C3kjjgXLXXk4RCY8CuX/QVcvaTimaVj33qqsOz6LdMz+xv/OgzIZ9DlGM2ANIKifwPLwxVA3xSNMGazyaxGvfWRqRweh3ievzYoY6xnqjLipNlCnUz1i+oEYiJyUb6JKfTKdLpirAMxXZKv4OolFfwAfWzpuOypO42vKcmsx0X5jlNn6Ck050BjstM4XoFmNPN2YBOL1GIqh8chniSdMyTMrsyM9PrCOFPmr6gTvSyK9p7PFVwx1oGYVyx14vCV/b5cvaSEWXZ46bjs+OgcUxIyLLsS//D+YpRUgPcOPU9zDodXoPdG5r0Yn9upHB6HOBcmfFF2UicbUfbyTGMIPGkBlAHtNAJEox3Q0uweURdkF1T6E1pItCNaCgviW9TTqGn6xFEAYjbHxZMCV4molWg+0QiSG1S69eZcTN8vW6zfzE8DpKUWNII0RD0lRsaT6gegV7sOLdrhpWNfbCV2pJi23cF+JQvcwOn7y7GK8pPUZQv2Bc+kHOIR1Zlkk+AQKSBRkPS+cQXqaNEiFQItVPGkBQLjQIaFslkEbEqRWkzl8DjEuWAsnfJd6FAWPaAmCjHJoerIQgaUkRtygwNTpe4eiV0kswuetItop6xsZFJRbk4j0fRJ41LFwka7ZAppAkmlDEG7I7vi0N6cD4APVMmAspifRhwiynBoDhF9NKU07Q9A73sdWLp0ivbCvrEgwcKwskC4bn8xEPsqlIhdosPm/jUWHTRPONAOIDxMDULjeAUeRqOst45NMNIuwSwFd7PHqE2MQ5wFjiMG0qHcSzyUNWjCEaePzjqE+Go6+skOmt0jqYvU3qyiPYVVU4oJohvp+TBIMQUUPMoU0gSSunn9eItopwBjbz7JHu8ywTjwHKh4plBSe0MkO4AKowIW6/uTXvN0sKS84p9ZOvbF+iManrBIZu5DiC8iNdmf3F92J6pb6tIbS5OnXf6AFhLtFErDBo0I0luBpj+eEmSsZxKb8GSrWYrxvkzt2Su39vBJioeyhk3JB2kQPwGwFUJMVwFj6h5ROhF5SigN30R78pHamwbdyMAHFw3V+Ek5jiVZQu1Tq9wgAXt06M0HNJc6EuMM1N5S94aAHUICRnWk/weA1ynSBWg9JR9bOvbFJgKeTDypcvskEI113f7CjZEfjpViJhUzhypgTw6Exu2N2ytAoG5hCluw/0aLcq1yTRCBqlkKbE2RWkzl8DjE6Rw+TDjQzY483BDqg9tL8mgXqXsEQ9Sh8JRQGr6J9uQjtTeNyIwPohHIcSzJEjb37nwgQ5AsuwQYe/PpLXUkxhmovaXuDdGbkjoW/f3h+14ElhHxEwfylMvOjtwXqHhqW1+/v/KU0BtLDtFOFYJk2SVs7gUaBytAC2BHBZEQjZt7VlHjobLem9XA1hSpxahNjEM8CQ5iig+V53UvvDWKEsBBRxOEdHDT0FIpyE2k7pHURaF6XaI9+UjtTSMy40PgBmcFh0WyBHr2VHSP95SGlkM00kfQXpRKJTHOQE1DKFRviOQvFLDo7098/Ut5ZunYF5sCGTUDicnKAQtPV9qvtI80gmSvJ9Oc5NgooTeWHGiP/g9TQ9AIH8i9FXgYLcpJpT8EPDmKWQruTo9RmxiHeB6cPx2ddBz3gsOdZgtVAdM5BlCZD5QBZXRRr2ZK1MBfoSQPusDOt04+Ud1mkadBVTR9KNAIYkwIaVwJm293PnECNEqFwJigOR/IEJpLLRiTzmMV8eulBnBoDrHNaN9uviN636t5ZunYl6eFn6H5cRkW7RGE6/YXEZJdQ8QuEDgxygBCiik1nnMaqVLGk/6ExsEK9KJBhcCOtNOYVPkwztZuFmK8KVMb9oJ9xVnEKECpBSArK+bhQRQ8nSSdY0I3wPQo1jAlwC7N7hF10cUx6MKrAc/kk9R6GjVNn7QOstfjSqDDQGVfQjuXVzG3Ht05F1O11CQO9FCtl5rQCHbtJloPnLTVwOsU6WLqvZhfPfblsYlxYEGczeWLc/c3wb5Ck09d2Ao0VnKIajM1kj9Rq2Q6xBVoRpMdz9gRSEVTWkZYKMzvkbkU7FSRWkzl8DiEMWYFnKcXkT4uGPMMLsnG3ALn6VmgBmMx9R1U3zWNeR6XZGNugfP0ROI/brsemxPBiSpSC5dkYz4E56kx6+OSbMwtcJ4asz4uycbcAuepMevjkmzMLXCeGrM+71SS04/tDpB+cRhDQd4bWV3wTNGOscU7M2BCr38s8tfkrpzeLo7NYW+vJ1dsKfAWRboGrq049rsg9Cr9f/Lnn3/iueYuYFbPTOyKM4ZQiwckCNiMCeMVw70FePEitVirJDNdD+8Tj35Rvn/JwJ/YHzgBpx+a0wPW6H0PEKf3gqk+5DVzeGbFVgPvUqRrwEKd+PPcFc7YDM/P8y3O2BXbcUXMDwBrUqQWC5VkFNQnSzK6p6Ov/5bNgcNx+nk6PWDiyfix+5OhTuEFc1jhNU8E71Kka0CGnvhzoHdZ/CfneZPXbPIu7/5isCZFarFKSdbmPbOFqL64NYoSYFhCh2gB0QiZRgGLJkYBhX9r+SJ+CCimb4di3WBHMg6oVhLj8yMLKaZvYhfdm/KXhW56RxpB7K4JyA1MToPEaABxZIxvqpgDfzwxec0foJVq6sVo6tVzELFV8edXbEGunh7WhFtzCnGnokohnhMMSiOJ+zg+imiNE+YmUu51TAcgjqtQ84ckdj98xmIQUJ9zwCbAdZNDXEatWwrI90oBSR0EzrH7IGa0SNU7UsWTFgiMAxkWyp8K3rFILUZtYhziFHBMsbvap2LdiY5OvakpMmQdJp4JCPRJecIukilgtvRRXwBB40KWXShIlCmkgM34MGpusMtHKCbV6C+ZPrJHYnfKx6YB4Nxb3mbMgT+e3Fa2AkSARa00Uo12CcmBKkk+8dUk06e5YquBeRbpGhAfy4InKdajpMWXSiGdk6RSjvsFu3xEMkKGhUKzY7RjFJ4WTYz26COZPrJHzu1ONZ1/CHwvygACHTAE7XCgHUBQIsgeA4pmkGMxKesdYxOMtEv4bLg4PUZtYhziebCF3Im4T4fhgSDa4Bg5xVcTBR0s2Wd8YIFAI4iHVcg5yhRSwIfxQXQjcoac/KWm+JHYPbmpKYUFdTRGEOpLoY4JaCHRTgFGdaQxCgn16jlE5Aw5vZpU+qR3XJOH7/sMXAcc7KhSPkZc/KhSSOckqXim/QLRjcg5yr2O0TkS7amvVPqk0ckp3YsSUNgUM5VDxZQ/oIVEO4XSsEFjCnIsJuW6b2zCk62fzfg1p5bg0pXifkSZ+3QKiMbPdHVkngCBpuQTVckSkg8yQR8eAe1F+Sb2lSyh9olq/JwhdD+S2BdvV88HTwmlIRCbkpvUmWmQXcsLxv4Ygp+uJAAIdJYFxF5Nh0h0PrBiq4F5FuklYDgsWlH2kxZWqoRkT+rkUcSe0ogn97fXkSdw6/SDOIEDhyQ2HTtj43Ou9yK0S1DMpPJlBezJgSTjWB3HjHKtclMQgepngzctUoupHB6HeJJmhpy1N9p4CbRvg+SzIoE+UZUsIfn0EiMS+0qWUPtENcVvEvseSP7YlNykzkwDwFnd1VcCfaL65T30pwwBo9dn46vzBuTYS7AVFP2b6HxgxVajfsFLwXDP5GlaWKkSkj2pk0dR24onukDodYQnIhclECdw4JDEpifPGHwIZHVJr0O7BMWMKgTJskvY3AvJOFAhSJZdQnKuVX7s6H10/jDwpkVqMZXD4xAnkvZpF82+MkqAMSUGUzH5gKhKlpB8Uj4nlcS+kiXUPlHVJGlvEvsmf80n+iRiU3KTOjMN+Oxa3kl/+OD2gYVqDZrQN/WK0KEoG9FZ47KJKoRBwNXghC9CCyKeXJa0sFIlJHtS034NgFucfK+jIhf9m2hPfRW21xc82b0GnuioLgpCejGj/0y6kWTsqTMxo5xUzRlPyFv7J8OX7TFqE+MQJ5K2bS/8qBW746Dok5f2WycGMkekShnPzf0LqOyiJgl0iCoEfgynDChHYJwPmFS8i859fAWR+sZ3l5x8EmhK06M9qtssRtOIRnakShlPukmd9Gc91huxF2UgZwk9h6JvJKNWKcrNjmsS3/cKsCZYVcoQdAyOkRZWqoRkr1VMQHOIpyiBfUSTdhb0OsJYHwAAB7140yfOKpGa9nbn9IoSPCXACKG+eaJDVOP70kiVMp70J8nYU3sxoULg0sn5q+dPVT7pZT+V8TtOvf/LliltG46sjtokvLhF7I7TDwtzgDLhoDgNaXSgLmqSQIekQiacBo2RvQGTigWBSmQUyRnIHyPSUvtE6unRntTxNADjEPatlzeq8/5MXcIEFvSMvZoOkehM9q7YUmCeRbqMbW2+wEIV06E8BWlhpUpI9qb68CgCHoN4ckCvY30AAM+nLLVPmlWkbtrVHTw854BNQDdPcogqX4fQjiGSP0nGgdqMKTuevb5owoLQCCDDQuHAoXoL8OJFajGVw+MQpgbHMZ4z8zzM4aKYFl4f45tnfVySXwE/yeozoD7rmbPAR+ZP/dR8Fs7TG+Kb5+1wSX4RKBhYKOKsOBF+P/Zn/4dglYpk7oRvnvcC21SkFi7JxnwIzlNj1scl2Zhb4Dw1Zn1cko25Bc5TY9bHJdmYW+A8NWZ9blqS+TdB+kPEw/APGotizMK84KD+9f3z02f+2g59GSQyTtVT0vmsO6HH1XfF1fM3r2F8SKYO0KXnbHFcks27cPVBRT3W79AgH67K6Jh+z8a/GS5Ki7eoRr4rzAzjQzJ1gO58zpxm5l249KCmovhMjaxLMkA05FpRKp4Z7mX4rjAzjA/J1AF6i3PGpNVP9JTDFJgtJGZ+MX07FGuwp8/v+jc3/f4v3i8MovjwAeMJGHMWOFpFugAc2sNfixMPS3KdZUwiPKnKAdAC6CPqLFMQCg9Tkm6i6QNK88RdUY+LLnEUDZGG5nLRiCeFh/M3a4L9KlKLqRweh1gEHtM6hynoNolpA0FXA+SHdgRRfMnwVHBYIKsvZKTKYALGnMil5wrnFsc7loHSsJ+YIwSRFRCtdZYxifCMRqDsiw5SU6GSDwV2BM2UlHNU68oHo94FMqA8eIs0bnMaEDSW7JpSMw5lsz7jzZrayLfYbx5TnWOq2+lt22GJ76VjPWmXyoA0Ikngz1SRvTcBqsachc7hFaCoIL7KDM6zSsJe0BGhEmwaZxmeyQFAhZEOxdRBQSigV7JTJTQWpUOazORdkcatpwFoIdEehbojVbM42KwitRi1iXGIRUjnUqqEZFftJLRD6NmRALUdT8iwoxUyHVKc6FmrxpwFD95FoBjH8w8wnKrCLhAH2VGUn/SyTALLXoLR+KEB6HNDQkEkJDtV8TDg3rtCQrI3VQSBKmCXgwR6JtUsDjarSC2mcngcYhF6x1RCsvfSaW+aQUbGAjgwdeFGZ94U0bNWjTkLHtSLwGHm8RYYDue8KHtQatT0skxCSs8m8CRF/0ZBJCQ71ZqvWBtF/2bvXSEh2WsVgmTZa2Hrl+OYxcFmFanFVA6PQyxC75hKSHakTXwvqWO74kQ3yMhAVmWovLzkrBE3X+ePuQodyCuIB54cPsZIlkFJjmE1qLImOQyAG/NRKIiEZKfaow6o6RGpaZKyj8eVCv9Y0RVNDhLokFSzONisIrUYtYlxiEXoHVMJyQ4Zgq4GyECyvhBEO+tuLQO6MSxTSEk1mIAxJ4JzVaRriKU0fUHcRYxTEzNLcswaGMHW/t9ypbJHmlkmo4Rkp0pmAgIY5++KFKSn6qVkpEo5CqkjVbM42KwitZjK4XGIRegdUwnJThUywR2BJ42ARpDsyC7alW8Eqq4JAAddOoMJoMvgbjJmFzhXRboMnf942vce44f+GkVZlpJIDkBGVjIhu1AQCclOVTwMSErzxF0xHjeq6PLVc4N2TEYOEuqOe/fCvB5sVpFaTOXwOMRngBOP01wUY96QO+SpMe+OS3Kb+K9DAPU4ffE15r1wSTZmfVySu/AfmojrsXl3cIyLZIxZlXGe3rokG/NJOE+NWR+XZGNugfPUmPVxSTbmFjhPjVkfl2RjboHz1Jj1eYOSzN/VicM/RtIPAROleQL+GTaE+FO/Z0AEBsHzlIAJ/YXaschfk7tyers4NocVZr4IWIciGWNWZZynUzl8daqf9ZtgBHnyZ/IqyafwgmqB+JhzUXYSp/eCqT5khTm8NVi9IhljVmWcp0uUZNTRU36DdLeS/GT82P3JUKewwhzeGqxekYwxqzLO0yVKMurxk6WUDEoyr3tWXBK/XBZT+E/iqTxMdqxrOTsShKW6xcsB1Upi/PhP8cX0TeyiDzTyl4VufC9AI4jdNQG5gclpkBgNII6M8U0Vc+CPJyav+QO0Uk29GE29eg73Aa9cJGPMqozzdCqHr051xH946c/wsCTDgSprD2UI6vU1/GbXRb+3Y0RBokwhBWzGh1GVCXb5CMWkGv0l0ydWOBG7Uz42DQBn1T/1HcQc+OOJJnkCRIBFrTRSjXYJyYHqHcD7FskYsyrjPJ3K4UtTnVcn7uWoUt4L7m70TTByuqCp4plu/7o8UJjvGJFzlCmkgA/jg+hG5Aw5+UtN8SOxe3JTUwoL6miMINSXQh0T0EKinQKM6khjFBLq1XO4Dzd/fWPegnGeTuXwi1Mdw6VLfxKUZJX2hC7upMI/fu3TtS4HCdEBz17HSOwrWULtE1XW+ER6u9gXK1bPB08JpSEQm5Kb1JlpEEygNG+gr4LQIaljfwzBr+YSAAQ6xy/9sVfT4T7gxYtkjFmVcZ5O5fCLUx3DuSSn+E1i3+tK8sNpADiru/pKoE9Uv7yH/pQhYPT6MHx13oAcewm2gqLfg7u9rzHvyDhPp3L40lTnt6WibEBNN+wkuL73luQ0ulQ5SKCD1F7HSOwrWULtE1UGlL1J7Jv86xehPRKbkpvUmWnAJ5ZtdQG9mDP+3E1YqNagCX1TrwgdinIDBmtljFmEcZ5O5fDVqR6/CaX7ehe8xIvyk3RxRxWC/p0TMoAgBwl0iCoEDQcZUI7AyFdTRwl0GKh4HS0FgkQ3kvrCWS8iOfkk0JSmR3tUt1mMphGN7EiVMp50kzrpz3qsN2IvykDOEnoORb8B8fWNMWsyztOpHH5BqmMIotsfQO6V2CbwL1F+8nVt/7ygk7p5faGvZXKQQM+kQibqmEBFgR1PdZRAh7Ea30hGkZyB/FXJap9IPT3akzqeBmAcwr4okClIVOf9YacKWHQFPWOvpgPYe5DeFLxykYwxqzLO06kcdqo/BMUA935RzBmw1hbFTODlMmZ9XJIvgV/I4vcwfTE1p4DvtXf4ansizlNj1scl+SpQMLAyxPX4RPj92P/qsBcsWpGMMasyzlOXZGM+BOepMevjkmzMLXCeGrM+LsnG3ALnqTHr45JszC1wnhqzPu9RkvVD1Wf+qAd9GSTy7+F/4Ol0ZsZC69ghwb/uLkpgb5xnwEDN9+rNzbweb4Qx6zPO06kcvjrVUY/1cxfIh6syOjZ/NtMrJ1fwcKwDk2mWvQNxrsAleR28EcaszzhPp3L40lRPpeWZSuOS/HpcktfBG2HM+ozzdCqHL0113OnP/GN1ZKYkUxbRv5i2f+vGs1i/YUf4by6lFSNSBbSMx4oWDV0HIcXU+q91juM0fyTNLg/nD9K0sRoy4kkfNgHNLTkkdX4gcxisYZGMMasyztOpHL401XGno4SwCpLSsB/c+6pPkVgeIOjqVzkBENQXsuyCQWK1w3BS0RcqhIdjRQfQDALgMzOfZpwoC/rvmr9UvEW0Q6jn1uxIdddAVM0BuBHGmJUZ5+lUDl+a6risEV9XNi5lXtkHQEeEirBy6PYH9CSyY1AIxRrKZ4TOqhmpC2CrYoLSsBHtFGDsBUn2wXyacZJK6I8mqrUPW+lWTAENlzpqbnKgvecPoA4GMofxehqzPuM8ncrhS1MdxTjVYF7ZRdkD4rAGJ1K1YJEQsKNXnEOzWqQgLEUJGJNbPVZ06AUBu+aDUWp/Nolk7A2NJn5IAvqcBNS9N7cUP/pDSAwGMofBShbJGLMq4zydyuFLUx23c7qRMdx1JRmC5Fg2DpTk2EVENwiSZZcAYy/I3vkcK8nNoSPwJ5DVvTe3FL/n3wSepOjmEF5AY9ZnnKdTOXxpqvN7ZFE2oOpm3wWufhSAogRUHlLp4tC0xznUUwIKQlV9qQq59caSQzTSR6QJJJUM4jz0B72ha+AGZ3VPwaWm+Mlf9gFwg3NRzH64EcaYlRnn6VQOX53qsZTOfKnqMVOSKcgoFYL60k5ZKEjRt+E0VUWW23gs1Z5mEMgQ9I8HXz1bWwBjjCP/KAvOgcHJNnIemgKNQL0kwAihOTcIWkNEhkr/XQMV3ewnrqcxZk3GeTqVwy9Idd7gQHc3gNwssT16/vG6Ry3hQIB2VTUaAYLgSaNo1gxNG7ApuvXGol1VrQ5CiqkzH9CLI0skTkw0h2axFLSn7mwCcW70ISq9bJofCMB5174bgjUskjFmVcZ5OpXDd0t1lApUhaIY8ya4JBuzPi7Jj6m/0jW/aBqzMi7JxqyPS/IU/AdY4nps3hEc3SIZY1ZlnKcuycZ8CM5TY9bHJdmYW+A8NWZ9XJKNuQXOU2PWxyXZmFvgPDVmfVYvyenHqQTG0ryT+FdaQH9EfRGI/4JRyN6xuLBFCWiJXjPtJsfW7Viv+4DFKZIxZlXGefr7JTnxV/X/QjFP/M9WANakw9V9hpWLRK8kX70mM7i4XkFzu40xSzHO06kcflmqs4ocvqlRy1OxeabAz/B2JXmRCa+8bu9Lvd3GmNUY5+laJfnJ/0YHug/+Q4wsA4KeNLJ6EdijJ2s8LfonX5V52vGkCjsdAC2RXhBS903+vbHSitEI2LFYNxiBqFcdJ41Lo4hBANxkjMuoz0YDfzwxqMYFaKWaesVdwLPncHOwDkUyxqzKOE+ncvg1qc4LnRfuMVQSeO8nYNfFzXoDgTe7qiPtSYWgAkA7HGLl4JxlBOioIKIXJMnqS3/Ze2NFGQ7oLhlQFjEIaMZJ4ybQNLmMlAf+eHLL2AoQARa10kg12iUkB6q3BYtQJGPMqozzdCqHX5PqqAG9MrAL3vtEAXV9E93pFHSbN9XNa2THE00Q2EqiP6FzHaTXd9JfarLHuigUBHIvTho3wo5C0SioS7TTQqKdAozqSGMUEurVc7g5XhNj1mecp1M5/JpUxyjNMvAMiKmvboBVR+hyx5MOPVUC7QAqoskePwcIGIv3Ri9Ir2/yl4ou8aVkR5faXpRv5Ay5F0dCaajYu4xjf0ybn5wkAAh0jp/SYq+mw83BahTJGLMq4zxdpSTz1i7KecRLHIJk2SVs7rmWSJVAO4CKOcueymGTXpBe3+QvFV1+sSSjSa3ylECfqH55D/0pQ8BkMCUIka/OG5BjL8FWUPQb40UwZn3GeTqVwy9I9ZmSNqZ5X8uYyg9U2uVAe0+loIJBNTooIB2a9IL0+tJBdqnJn2oUSFJJjNmLE30S8Nm1jJP+8MEBgIVqDZrQN/WK0KEod2WwgMaYRRjn6VQOvyDV/9ooylEQIV3ZuOsZVsUAMm92qpTx3NxzNZJKAdCusHKgEWztP4YTvSCU6770x5P2qKa+kuEQZUBZpJjNOMknsncZJ/1ZjzUT9qIM5Cyh51D0uxLXxBizJuM8ncrhF6Q66gHu5aJ80zSO4eUuYnfc+MW6Xd944mZPt3lPpaDgKh5yoIoJ0wHIKHpBSN03BU+q/FMcGgEHKtZvUhBQx6l9InuXcd4fdqqARVfQM/ZqOoADZ+ZjwDoUyRizKuM8XaUkL04sBoc5JchHwpUpijmK19CY9RnnqUvyFC7Jl4Lvtbf9ansizlNj1scl+QRcki+Ca/JH+BMwcxjnqTHrM85Tl2RjPgTnqTHr45JszC1wnhqzPi7JxtwC56kx6+OSbMwtcJ4asz7vUZL169gn/9Lnnz9/lxz/lgpyVGdQFzxTtGNs8c4MmNDrnx55L/zdcFFa7FqBs5ZrJg5a6YDnCis5z3jBjTErMM7TJUoyyrAq8V9//XW4KrMgFeW7KvA/QHHgej39Rj49YI3e99c5tySfxcNBf2VWZzFecGPMCozzdCqHL031dAk+cyeilqeCpAJ/IOwzM2lyesDE1fF34ZL8esYLboxZgXGeflpJbv4XJxiT0CFaQDTq334FLJoVBdYbEj8EFNO3Q7FusCMZB1QrifHxgsVa7Ujsov80pvxloZvekUZybCaKTGgEHKJYWzNnZDypijgi4KBypvBwkrTguYXM3fGUUaQ4UukMeosZJxOX6PVgAkUyxqzKOE+ncvjqVMdlByjre+0BdDniPi2mb9L1Clm3p4oHfdKFyy6SKWiG6gsgaFzIsgsFiTKFFLAZH0bNDfZ6lRSTavSXTB/ZI8dmEmU4wE1y018zT7MlyUgVOyU7hTSrpEKQP+1SJcAIoXcG6BDl+AqS6VCP/lv87ujGmBnGeTqVwy9IddxxGAXogjsMb0aiazRdtTQSNVHQNS37jA8sEGgEzdtZzlGmkAI+jA+iG5Ez5OQvNcWPpKamimcvcrJrBZIdQIVRAYt1g8aiBORMAd2jPambVxF6dsAmEu3qKDm9glQ61KNTfT1xksaYNRnn6VQOX53qiB+/gpw4HEKxxtfXJS9WgabkE1XJEpIPKlD8MEF7Ub6JfSVLqH2iGj9nCBjpRmJfvF09HzwllIZAauqpvci9FejNXB3pL/ThTEcCyFlCsidVwtg+PgOSJxczqa8HoxfJGLMq4zydyuFLU53XYlHOvtcUTQLtkKWqSQJ9oipZQvLpFaRI7CtZQu0T1RS/Sey7VEluzlwdi16BVgJZzhLo01MlDOwQJMsuIRpdko0xZzHO06kcvjTVzyrJzY4ySoAx3bCcQPIBUZUsIfmkt0gqiX0lS6h9oqpJ0t4k9k3+mk/0SaSmntqLLIEku/yFAha9A3wQQc4S2NpTJfTsiPnwDEhOr0A1OtCe1NfDWRljVmacp1M5fHWqI77+lRJC80vVDPw3z3gnIpQiowmXKYR4w/IapUoZz839C6jsoiYJdIgqBHwjpB0yoByBcT5gUvEuWpn4CiL1je8uOflEUtNAbUYGcIgyoAwfQFkzjwEFW4sSBq2F5FCrELQdGD11jwtIo1QIaY8Y4eFiJvX1YPQiGWNWZZynUzn8glTHEER3N4CsW3US+JdAG7E7rlFYeJlSJrxJcQvXV6q6qEkCHZIKmXAaNEb2BkwqSwuRUSRnIH+MSEvtI1LTWK0jExpBWoF65imgYLEUyVkCnQcqZaICTCOecOBeENpZies9+go9sZhRhfPeo/s8GL1IxphVGefpVA471feCyx2XclGMeQnOU2PWxyX5Feh7GFXU4/T10ZircZ4asz4uyS8i/pu567F5PTh4RTLGrMo4T12SjfkQnKfGrI9LsjG3wHlqzPq4JBtzC5ynxqyPS7Ixt8B5asz6vEdJxhDkmV9zxh+/RkrzBPzDaQjpJ6eHQQQGwfOUgAn9TdnpkfeipeuxawXOWq6ZOGilA54rrORhxutvjFmBcZ5O5fDVqY5SCiQfrsrP9CUP68ouXnDFIz7mXJRf5dySfBYPB/2VWV3EeP2NMSswztOpHL401XmV6058pijerSQvVU5ckn+d8fobY1ZgnKdTOXxpqqer/JkrclCSGZZjEcilbXtBwn8HhkXTmOyY3gKwI0FYqlu8HFCtJMbHGxVrtQWxi34GLf/0H33kewEaybGZpJ9c0wi0dKSeOSPjSVXEEQEHlTOFh5OkBc8tZO6Op4wixZFKZ9BbzDiZuES/DuZTJGPMqozzdCqHL011XnC6B6keu+lwh+JiLcpPeJnCgWqsHxDUCzLtup33dowoSJQppIDN+DCqGMAuH6GYVKO/ZPqkOkqOzSTKcICb5Ka/Zp5mS5KRKnZfdgppVkmFIH/apUqAEYLOVbNjlOMrSKZDPfoiLDUZY0yTcZ5O5fDVqY748fqDergko28C9yaaeJkqrC5fWCDQCOqbmsJ8x4ico0whBXwYH0Q3ImfIyV9qih9JTU0Vz17kZNcKJDuACqMCFusGjUUJyJkCukd7UjevIvTsgE0k2tVRcnoFqXSoR6f668Q5G2PWZJynUzn8glTHEITX37FrDiWZBbgm3Z5S4a8vPYB2CXhKiA549jpGYl/JEmqfqLLCJdLbxb5Yt3o+eEooDYHU1FN7kXsr0Ju5OtJf4NMYffSxDMhZQrInVcLYztMlkoPkycVM6q+DyRTJGLMq4zydyuFXpvoz1xyu0VS0RAortVdX5CAhOuDZ6xiJfSVLqH2imuI3iX2XKsnNmatj0SvQSiDLWQJ9eqqEgR2CZNklRKNLsjHmIsZ5OpXDl6Z6utdmSlEPdET3ovwkjSKVX5toBFLlIIEOUnsdI7GvZAm1T1QZUPYmsW/yr1+E9khq6qm9yBJIsstfKGDRO8AHEeQsga09VULPjpjxXGmScoBRcnoFqtGB9qT+OpykMWZlxnk6lcNXpzruyvjnM7gBad/LgZIMGYL+yRQygCAHCXSIKgQNBxlQjsDI11FHCXQYqHgdQHsqEiT1hXP8n+TjkqaOJDUN1GZkAId66QB8AGXNPAYUqnZEPrWQHGoVgrYDo6fucQFplAoh7REjPFzMpP46mEyRjDGrMs7TqRx+QapjCMLLkeAq7JXYJryIa3BpptszqZvXFxgOz+gggZ5JhUzUMYGrHHY81VECHcZqfCMZRXIG8lc5qX1EahqrdWRCI0grUM88BRQsliI5S6DzQKVMVIBpxBMO3AtCOw9bvUdfoScWM6pw3nVWrwCTKZIxZlXGeTqVw071h+Byx6VcFGN+A+epMevjknwJ+h5GFfU4fX005sU4T41ZH5fkq+A/1RLXY/Pr4BwWyRizKuM8dUk25kNwnhqzPi7JxtwC56kx6+OSbMwtcJ4asz4uycbcAuepMevjkmzMLXCeGrM+LsnG3ALnqTHr45JszC1wnhqzPi7JxtwC56kx63NOSTbGGGPM85TK2sIfq40xxpglcEk2xhhjlsAl2RhjjFkCl2RjjDFmCVySjTHGmCVwSTbGGGOWwCXZGGOMWQKXZGOMMWYB/u///j+4YUL9wzteEAAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 02 - Mode Estimation Method\n",
    "\n",
    "<br><i>Prony analysis is an extended version of the Fourier analysis where you can estimate the frequency, damping component, and the relative phase of modal components present in a signal. Prony method estimates a discrete linear prediction model that fits the records [2].\n",
    "Suppose that you have a system that is affected by a stimulus and once the disturbance is removed, the output response decays in the form:<br>\n",
    "    \n",
    "\\begin{equation}\n",
    "\\dot{x} = Ax\n",
    "\\end{equation}\n",
    "\n",
    "where x are the states of the system and A is the dynamics of the differential equation.\n",
    "The solution of (3) can be determined by the eigenvalues and right and left eigenvectors as follows:<br>\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned} x(t) &=\\sum_{i=1}^{n}\\left(q_{i}^{T} x_{0}\\right) p_{i} \\exp \\left(\\lambda_{i} t\\right) \\\\ &=\\sum_{i=1}^{n} R_{i} x_{0} \\exp \\left(\\lambda_{i} t\\right) \\end{aligned}\n",
    "\\end{equation}\n",
    "    \n",
    "where $\\lambda_i$, $p_i$, $q_i$ are respectively the eigenvalues,and the right and left eigenvectors from matrix A. The idea then is to estimate the exponential terms in (4) so that we can determine the best fitting function.<br>\n",
    "    \n",
    "#### <br>Prony Python Script\n",
    "    \n",
    "The original script was written by D. Trudnowski in 1999 and was adapted to fulfill the requirements of this project. “pronyitesla” script has the following configuration:<br>\n",
    "![table2.png](attachment:table2.png)\n",
    "\n",
    "#### <br>Prony Python Script Variables\n",
    "* acoef: is the A matrix.\n",
    "* zPoles: are the poles in z domain.\n",
    "* sPoles: are the poles in s domain.\n",
    "* Res: residue matrix that consists of the product of the right and left eigenvectors. The residue matrix is called as such because it is linked to the exponential part of (4).\n",
    "* E: is the energy described in equation (1)\n",
    "* x: sorted energy array from the filtered signal.\n",
    "* ii: is the index from the sorted energy array.\n",
    "* R: rearranged Res matrix following the ii indexing.\n",
    "* P: rearranged sPoles matrix following the ii indexing.\n",
    "* that: is the time vector used to plot the Prony output.\n",
    "* yhat: is the Prony output regarding the pre-defined and filtered signal.\n",
    "\n",
    " \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pronyitesla(t,y,XX,tstart,tend,T,shift,tstplot,tedplot,plotFlag):\n",
    "    \n",
    "    '''\n",
    "    PRONYITESLA\n",
    "    DESCRIPTION:\n",
    "    Prony analysis program for fitting to a ringdown.\n",
    "    INPUTS:\n",
    "    -        `y`: The ringdown matrix of order N by NChan corresponding to t, each column is a different signal.\n",
    "    -        `t`: The time array.\n",
    "    -        `n`: The order of the estimated model.\n",
    "    -   `tstart`: The tstart(k) is the starting time for analysis for y(:,k).\n",
    "    -     `tend`: The ending times for analysis (same dimension as tstart).\n",
    "    -    `shift`: shift = flag; if shift = 1, residues are shifted to t = 0. If the data is noisy, reccomend shift = 0.\n",
    "    -        `T`: Sample period for analysis.\n",
    "    - `tsttplot`: The starting time for model simulation.\n",
    "    -  `tedplot`: The ending time for model simulation.\n",
    "    - `plotFlag`: if = 1, plot results.\n",
    "    OUTPUTS:\n",
    "    - `model_Poles`: ringdown pole matrix -- column k is for column k of y.\n",
    "    -   `model_Res`: ringdown residue matrix.\n",
    "    -     `model_K`: ringdown offset row vector.\n",
    "    -  `model_yhat`: model signal matrix.\n",
    "    -  `model_that`: time vector for yhat (starts at tstart).\n",
    "    -      `lambda`: Poles in the S domain.\n",
    "    \n",
    "    LAST MODIFICATION DATE:\n",
    "    WRITTEN BY D. TRUDNOWSKY IN 1999 AND LAST UPDATED BY FERNANDO FACHINI IN 12/13/2021\n",
    "    \n",
    "    '''\n",
    " \n",
    "    \n",
    "    import sys\n",
    "    import numpy as np\n",
    "    import math\n",
    "    import pandas as pd\n",
    "    import numpy.linalg as lin\n",
    "    import matplotlib.pyplot as plt\n",
    "    \n",
    "    #from scipy.interpolate import CubicSpline\n",
    "    from scipy import interpolate\n",
    "    # from scipy.interpolate import UnivariateSpline\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ------------------------------------ BASIC ERROR CHECKS ---------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    \n",
    "    if ((len(y[0])) != (len(tstart[0]))) or ((len(y[0])) != (len(tend[0]))):\n",
    "        sys.exit('Dimension error in y, tstart, tend.')\n",
    "    if ((len(tstart)) != 1) or ((len(tend)) != 1):\n",
    "        sys.exit('Dimension error in tstart or tend.')\n",
    "    if (T <= 0) or (tstplot >= tedplot) or (np.max(t) < np.max(np.transpose(tend))) or (np.min(t) > np.min(np.transpose(tstart))):\n",
    "        sys.exit('Data Error.')\n",
    "    if (len(t) != len(y)) or (len(t[0]) != 1):\n",
    "        sys.exit('Dimension error in y or t')\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ----------------------------------- DATA PARAMETERS -------------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "        \n",
    "    NChannels = len(y[0])\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ----------------------------------- SHIFT TIME PARAMETERS -------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    \n",
    "    tstart = tstart-t[0]                                  \n",
    "    tend = tend-t[0]                                     \n",
    "    tstplot = tstplot-t[0]                                  \n",
    "    tedplot = tedplot-t[0]                                  \n",
    "    t = t-t[0]                                              \n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ---------------------------- Set up analysis data and calculate offset ------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "        \n",
    "    tanalysis = T*np.transpose(np.array(range(math.ceil(np.max(np.transpose(tend))/T)+2)))     \n",
    "    Nstart = np.zeros([1,NChannels])                                                         \n",
    "    Nend = np.zeros([1,NChannels])                                                         \n",
    "    yanal = np.zeros([len(tanalysis),NChannels])                                         \n",
    "\n",
    "    K = np.zeros([1,NChannels])                                                              \n",
    "    for i in range(NChannels):\n",
    "        Nstart[0,i] = math.floor(tstart[:,i]/T)\n",
    "        Nend[0,i] = math.ceil(tend[:,i]/T)+1\n",
    "        spl = interpolate.splrep(t.flatten(), y[:,i].flatten())\n",
    "        yanal[:,i] = interpolate.splev(tanalysis, spl)\n",
    "        K[0,i] = np.mean(yanal[int(Nstart[0,i]):int((Nend[0,i])),i])\n",
    "        yanal[:,i] = yanal[:,i] - K[0,i]\n",
    "          \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # --------------------------------------- FIND MODEL ORDER --------------------------------------------\n",
    "    # ----------------------------------------------------------------------------------------------------- \n",
    "    \n",
    "    NdataPoints = Nend - Nstart\n",
    "    Ntotal = np.sum(np.transpose(NdataPoints))\n",
    "    nOrder = XX\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ------------------------------------- BUILD MATRIX AND VECTOR ---------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "     \n",
    "    for k in range(NChannels):\n",
    "        Ym = np.zeros([int(NdataPoints[0,k])-nOrder,nOrder])\n",
    "        for kk in range(nOrder):\n",
    "            Ym[:,kk] = yanal[(int(Nstart[0,k]) + nOrder-kk-1):(int(Nstart[0,k])-kk+int(NdataPoints[0,k]))-1,k]\n",
    "        yv = yanal[int(Nstart[0,k])+nOrder:int(Nstart[0,k])+int(NdataPoints[0,k]),k] # + 1\n",
    "        if k == 0:\n",
    "            Ymatrix = Ym\n",
    "            yvector = yv\n",
    "        else:\n",
    "            Ymatrix = np.r_[Ymatrix,Ym]\n",
    "            yvector = np.r_[yvector,yv]\n",
    "            \n",
    "    \n",
    "    yvector = yvector.reshape(-1,1)\n",
    "    acoef = lin.pinv(Ymatrix).dot(yvector)\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ----------------------------------------- FIND POLES ------------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "             \n",
    "    zPoles = np.roots(np.insert(-acoef,0,1)).reshape(-1,1)\n",
    "    sPoles = np.log(zPoles).reshape(-1,1)/T\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ------------------------------------ SOLVE FOR RESIDUALS --------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    \n",
    "    Res = np.zeros([nOrder,NChannels], dtype=complex)\n",
    "    for k in range(NChannels):\n",
    "        ZMatrix = np.zeros([int(NdataPoints[0,k]), nOrder], dtype=complex)\n",
    "        for kk in range(int(NdataPoints[0,k])):\n",
    "            ZMatrix[kk,:] = np.transpose(zPoles)**kk\n",
    "        Res[:,k] = np.matmul(lin.pinv(ZMatrix),(yanal[int(Nstart[0,k]):int(Nend[0,k]),k]))\n",
    "        if shift == 1:\n",
    "            Res[:,k] = (Res[:,k:k+1]*(zPoles**(-Nstart[0,k]))).flatten()\n",
    "            \n",
    "    del k, kk\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ------------------------------------ REORDER USING ENERGY -------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------       \n",
    "             \n",
    "    P = np.zeros([nOrder,NChannels], dtype = complex)\n",
    "    R = np.zeros([len(Res), len(Res[0])], dtype = complex)\n",
    "    for k in range(NChannels):\n",
    "        E = pd.DataFrame()\n",
    "        for kk in range(nOrder):\n",
    "            if np.abs(np.real(sPoles[kk])) < 1e-8:\n",
    "                E.loc[kk,0] = (Res[kk,k]**2)*(tend[0,k]-tstart[0,k])\n",
    "            else:\n",
    "                E.loc[kk,0] = complex((Res[kk,k]**2/(2*np.real(sPoles[kk,0])))*(np.exp(2*np.real(sPoles[kk])*(tend[0,k]-tstart[0,k]))-1))\n",
    "        \n",
    "        E = E.to_numpy()\n",
    "        x = sorted(E, key = abs)\n",
    "        ii = np.argsort(np.transpose(abs(E)))\n",
    "        \n",
    "        R[:,k] = Res[ii,k]\n",
    "        P[:,k]=sPoles[ii,0]\n",
    "        M = np.transpose(np.array(range(len(ii[0]),0,-1)))\n",
    "        R[:,k] = R[M-1,k]\n",
    "        P[:,k] = P[M-1,k]\n",
    "        \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ----------------------------------------- SIMULATE --------------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------      \n",
    "              \n",
    "    that = (np.transpose(np.array(np.arange(tstplot,tedplot,T))).reshape(-1,1)).reshape(-1,1)\n",
    "    yhat = np.zeros([len(that), NChannels], dtype = complex)\n",
    "    for k in range(NChannels):\n",
    "        yhat[:,k] = K[0,k]*np.ones([len(that)])\n",
    "        for kk in range(len(that)):\n",
    "            for n in range(nOrder):\n",
    "                yhat[kk,k] = yhat[kk,k] + R[n,k]*np.exp(P[n,k]*that[kk,0])\n",
    "    yhat = np.real(yhat)\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # --------------------------------- Place output in structured array ----------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    \n",
    "    model_Poles = P\n",
    "    model_Res = R\n",
    "    model_K = K\n",
    "    model_that = that\n",
    "    model_yhat = yhat\n",
    "    lamda=sPoles\n",
    "    \n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    # ---------------------------------------- PLOT RESULTS -----------------------------------------------\n",
    "    # -----------------------------------------------------------------------------------------------------\n",
    "    \n",
    "    if plotFlag == 1:\n",
    "        plt.figure()\n",
    "        plt.plot(t,y)\n",
    "        plt.plot(that, yhat,'--*')\n",
    "        plt.xlabel('Time (sec)')\n",
    "        plt.ylabel('Intensity')\n",
    "        plt.title(\"Prony Signal Estimation x Real Input Signal\")\n",
    "        plt.savefig(\"prony_output\")\n",
    "        plt.show()\n",
    "        \n",
    "\n",
    "    \n",
    "    return lamda, model_Poles, model_Res, model_K, model_that, model_yhat"
   ]
  },
  {
   "attachments": {
    "table3.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoUAAAEXCAIAAABUDalWAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAADttSURBVHhe7Z1NliM5kqT7JH2AvENvZtMnyFXdpDaxmjv0GWoVu7lIHqX3NUIKXFJd8UPQSDPCzOV79eypqgkUMBgVGowK9/yPfxtjjDHm07gfG2OMMZ/H/dgYY4z5PO7HxhhjzOdxPzbGGGM+j/uxMcYY83ncj40xxpjPM+rH/2GMMcaYd1A6a58H/fj/GWPWxnVqzPq4HxtzfVynxqyP+7Ex18d1asz6uB8bc31cp8asj/uxMdfHdWrM+rgfG3N9XKfGrI/7sTHXx3VqzPq4HxtzfVynxqyP+7Ex18d1asz6uB8bc31cp8asj/vx9cE78mv64fgDsCasTfKf//mfJWp+KvgYlM7ax/34xPzP//zP3qXuz8D6+B2txj/+8Q+8lJpye4KnxD3eksS8C7yO0ln7uB+fmH/+8594R//1X/9V/Leya3LzRlynS8HCAejKJXR/R2Cmmt5Sdy7eBcEbKZ21j/vxiUG94R2h7PVFWWcBKKKvs0C3YpUyQjt+28aVt4D/qm1x8I6KZRaAVRObMWBxAbrRHtcdbRfvBcDrKJ21j/vxiWHhoVZVrhGeCDoIIrzFUapYuiz4u6rgP2UvDt5RscynSTUVuRfT7U3N152L90rgdZTO2sf9+MTcK+4/ULQoUdooS8TpxroFkOkWqzTKAIJwmQHcxvgDcAb8mtYh1ZRQZ4U9X3cu3iuBN1I6ax/347MSK5zV2CzOWMOAFc4/VqcaRhAuKz/KzOLgTRXLfJpUiSKW4XzduXivBF5K6ax93I/PCquOtRqrMbmpbmOFwwCMg+img8CsjN6a+Ti9wrmV1lcZ0mYcRDcNd/FeCbyU0ln7uB+fFVYdrrBh6GWhDUeXtoqTLio8yVj5+jN1OgjMyuglmo+DkrkVVagdFRqL66m6uwtdvBcBL6V01j7ux2dFVcfiVDUiIld1G2neQhBXVT5dxhkxy4LXVCyzACzMhOpovu5cvBcDb6R01j7ux2eFVYfKZANWNcbvzbpFA8QS1cEBMUdBxlt0eYsRsyx4TcUya6ByIyorMll3TOLivQx4I6Wz9nE/vjKsTJfltXGdXhIX78VwP/7p8Dt0+uO5uRiu00vi4r0Y7sc/HbxBv8TL41d8SVy8FwNvs3TWPu7Hxpwb16kx6+N+bMz1cZ0asz7ux8ZcH9epMevjfmzM9XGdGrM+7sfGXB/XqTHr435szPVxnRqzPu7Hxlwf16kx6+N+bMz1cZ0asz5v6MfGGGOMeZ3SWfs86MfFMsasiuvUmPVxPzbm+rhOjVkf92Njro/r1Jj1cT825vq4To1ZH/djY66P69SY9XE/Nub6uE6NWR/3Y2Ouj+vUmPVxPzbm+rhOjVkf92Njro/r1Jj1cT825vq4To1Zn9P0Y0z0119/FWcrv379Qh5ciz8E03FGXF+f/Y8//kCSxIs5t8FNAL9//z5mDZpxw1xv2fwZMEU9CyKcXUa5cTaw+GLtzCu71KwR8Oeffx6w/2d/xVq/jHLjEc/qH7JhAWMx7jYTKl4bvEvjRGDxxerz+X6MgtQuvwIKnqmK3wdzvWVGganTnwPYpYpzIJj0X//6F4y3P2MPzbiBYxbZm+WY2Q8AT1GsPXlXndbbfsCLOGCKXdH6ZZQbh/PUAmbEDzW14OGQNcGai9VnpJgZ/yKYgry4uWgJTILrw/bw9tdZ92Mws5L3Ep/r7c/Y5MVZPrvIY2Y/ADxFsXYDU5DXt6ve9gNexAFT7IrWL6PcOJynFjAjfqipBQ+HrAnWXKw+n+zHaFfv+tsq5EFThMFvyQwSRJAfMM7pCJqoZk89NX7BVQbASCKNJRCrH6c1MKIh/MOExNBEGQcCRrhgLk9BEJ8Lw/VcvLt5AZGUJM1IDWkusvcs9SKBInFjmRBGnBpQwyCfhcCOypgKaPZkpAxF3Voe4BAR9YeBeYu1D3govF8+Ka4lupU6DyP6tABuY4pT3HwLoBcncVLagp+K2yc+fIYRlFtnbi5MpPzxWWAzCDA8KvXJaS6PQVxlUBzhLRHnlZ63AFciQVzY5EqoEUmMJDESx4oYhM0gdvue7++JksH/J45g1ODFrQYWXKw+I8XM+NfRdhd/E3w3MPheGQR4u3o9iNefUdkYqI8CgM2EMUPSCIkFhmgZMYPsmAqRmAE21xkHSs8FKx5pPhfsmEf2zAIizSRxlgjj0oM4PD1LvUgJ4iIBbC4So7RCRODCYDbpU5yzRBRPRp0BIDhev9x66/ZGi9yV9LCbqfMw0ntxDz9FoBcXcVIY9YcHERoEGahpZq4XJuJEcpGKhhbGeZNLm3raiiutDApEitPVvIzDQEIKYAMYFGxYCeMkBemm2QGDtFNCLQAGN1ZjBwZkSAj7PvQGhmuK1Yjr7DFSzIx/nbi524ivhNn4SnqZY1y2jBhMLxvArd83PgSIJ3grZZDLKRjEcHw6+YlUvDc1BfUagJYd7c0LiIyT4Mq4YFyLTMMB72r4QwEisumKGOeoZhxXxoXiyagzjJdXQp/jmDVoN4q/lToPI70Xp3jvLfTixbkTc5bQHcVhaxSDMHqZKUhTEI1NpCFNF1fACInxZBTFF4wXJyB9epbUDjespNy4w2BxAlHcHIIrjeYCHhocApvDGWdwQWbWNlIc82xpczeAP0/FP6uir9BNH0ERZ4w22xIMNSd+ahPURDSwBmtgKpKmw13YFCCO68OpY4ZEvCV78wIivSQyyo0vUvzhs/QEGIt5aeAa14AlFd0dJFE2CuTKYFwoXhtJMFgePml04yfwYDB7sfYkbc5m6jwpIlcG47230ItzFEmp6g8PghjClxiNogggmLIlmh+JNGTsDj7bMqiMjOfFsusSlqGEyX1YZWLmqUEvYbl9By5kjEvQNKivX9yaYM3F6jNSzIx/nbS5G8DwGsT57qmJxBmjDT0/sr0e0EP6GuUkcTp8dAAE/AxBRvF46pghEW/J3ryASC+JjHLjixR/+Cw9AdDUcWEYpfxKIoMauTIYF4rXRhIMliegJMU/kGMmTZuzmTpPisiVwXjvLcy8nZgKhuwYpw0D2fCpg9HLHEcNuE1zB3YaMnBvA74vCdfawNpuujuw72kKJfp93vQsjMvANcbpwpCtuIy7vMF90A3YScw4Xd2iQQGBi/1PgqZBPV0YenFrwkWOGSlmxr9O2txn4eeyOHeYML7RcuOLGE8aDsSVLu06QwKfg1QSImVIyTGQHREuMsCWuDc1Ir0lxVuyNy8g0kuiWRgXKZ6GC8l6AoJbmhHAjieLxoKYRK4MxoXitZEE4+VFIIO4OEeBSYu1J2lzNlPnSRG5MhjvvYWZt6NUEDc/PHRZyIjQ7WVWtuIPgRJ50pCe21ueBDKKog9kcV6mKve+MsNICaN+psp6QBNnR2QyodwZg0NAenFrMrO8keKYx6s39ynwJvj1LqJgvBttzIgPBIw0OxtSTIhRgLY+Q3QFBL1+DHC3uQaAbErI5JoLwJbb+/hG4q1oI8m2BUSaSeIskTqOIYB281maArqYC65mj3eZgS5tXCmLLgyMYlxIUBtJALu5PBoMgjT8MOIa9uNdT1fnSRG5MhgHzbcwiAulincZlAt4pscCaWbmQI2K9D4SMlI8ub3l0Y4GB4qZeWHomIINYERBdB+uhHrSmx02DNyFMU6oTYaRzpbawF0YTEvqF7cgWGGx+owUM+NfJ+4ywSsZtLcIx8YXQ/h6mBPZYIP4tmAzkmbnhyYlVAYQ1ykeLri5BgAXt4pz3/CUp546LTgSbyXZ5gVE6iRpFtGMP3yWWkDql4IFUAaYBHdTtuhSn55dgtpIArrN5XFtIoonP8Ovg3mLtSdpN8C2Z6zzpIhcGYyT3oekFycxVfPDE2XxkwbqzDFbTfMjkYYM3PFnWwYHRmbm5S2AF4drLYjuTJWJ5uwglt4gIdcDqASaqDaUKomRje6aYIXF6vP5fmyMeQXX6bvgsV6cS4PWhT9nFOf8nOLFuR8bc31cp+8C39JAca4Fv7/yyyVAM9aXywtwihfnfmzM9XGdvg6/YF3pK2MNOhaekVymGZ/oxWGdxerjfmzMuXGdGrM+7sfGXB/XqTHr435szPVxnRqzPu7Hxlwf16kx63OCfowpRAk9SfzpwMh///d/46p/T7gUWNV+C0NmPfh9nmLsuhsP8/cEDwdOggyvJ+mBzLvu3otgbcVamGf3kP8eGMZnNz/OjquMySXFf0IFNIQZAF2BCP9p0rNTxH9d/Mbk8+h9NUn7oB8U/shSP0X9mDUjxcz4V0B+/TM/GC/+G7mzvLnD1hknOmzSHr0FvGVhb0lyXvDsxboQ4/P9eOJnbPLzxiZUnK8nYitiBhCTMPhsH4IehyfExX9r8nkG7wuTxltU8g8QH1nqp8CDFKvPSDEzfjN8K9ro1/f9LG/usHXGiQ6btEdvAW9Z2FuSnBc8e7EuxOB8/wjxMzb5eUNHYfcV+tbBDHDj91rYSfBwCp2iuGqudyV/it77asaxmPfuwylo7k9ipJgZ/y5e3/eUQS4NfiwI3jeDRJ9jgA9BibaeHXfT50ay3kDF8QmDG+dVqqQBlDE/YJDENfChtH4MBxz7+/fv+9Ab0KeEID416T1dXDOgJiWki+ukHnMhqDgHAm0FYCSRkuMpYoRzgTpPbwGJOr+CuKa7RDv5cPE7cfB020h7iE27bdOd+FEsoa83hYgGUqBNVrGQ5uZzrIgTEYzSZ4ZLkubFagIxeYIZMAqaEvrSM0IBrrzVAyukHlfYDG5Ozrjg4pWNQRAftoTC+0pw3uJUbF7qGWnuT2KkmBn/LvTB2kx6c3JpKDk/N8mljaA+0/pARFIQNiI0mgNjHLPwc6yFMR41sqlRPBLzQwCbawCwMYXyy8At2hqIIXBpi5gZKDOUqkANZEKtMM410Cs/jPikHKggSIshUSw3PjLjzTzUAMajRqQ8dOv8IiZpTnoMeqiV0R7S0P7oEwJgwJXNuAbCjpuMoD5mzc2PA+VqCIlvChlgawGw46uXgVu0NTA+QgTDEQfKKZgBAiVBhDMyEqcboORxDZuTcxRtJaRYqeJEMPRosBWPIKicNZuXekaa+5MYKWbGvwW+4xc3Pb05uTT0mWi6uCICg0ESZUTiaPcGRnEkxtNYudSk2Qlv0cZHFlsXP7gycJWhoBLGW6LW40pXxHgzIWCExHhPT0PPLqKeUFycgJLA7uWhpl4AXcJgcQJNMY9v2jOL34/mmldDe0ijfhFpD1M/uA1rvQXQ23zqS6hD1Ly3mgQfhOgPDcqAu8wjOy6gl5NADA1tJYz2U8lrlxoazBPjiMBgEPAxi/OFxNEVcBnR8qCRPb8PZwEPUqw+I8XM+Nfhi+TLeIX05uTKSPHkchkJBCkTLFoYuPIT0xuIJ4JxH/SNuAB97IhuySg3voMhGIi7HMtZtB6NlRGDsGtXIEN6OsJnERiYMiS3p+ddApdPQUFvG4v6C5xovKWjDSgJ7F6eqCFwsYDifPEwP+EsxelPWm7vDOYq1sJoD2WkOLYrfuQYl4ErP1S8GxlsfvNtJjApP4ecHWJctRjNLiMGYdfuAMhiWswLuDZc6aZ5b8M6YAjHEgyk+0pyaHBLQJPEcnvvK4EgchbnCyWhcV/a9n04C839SYwUM+NfBLv/ru1Ob06ujBRPbvp49dAHBVfWfG8glMhcnEBcgLIR3ZJRbnwHm8ZPLa5wkYF5uB6NlRGDt/GVK+qnA1BKrIEyqInuTd3R37UFuJhOgsn9FxhFYCsJ7F6eqCFwsYDiVNxS34GdxmJUdMGzi38vXOTiaA9lpHjaQ8Zl4Mpt593IzOZjICl+YL9qqpGSBj9+XDxcXDl1VMLuAUEN4puT38Z/xaWRQY3c3vtKQMOdjCgJjRf34Sw09ycxUsyMfwW8pzdOkd6cXBkpnlx8AmJ8AGQU0+0NVObifxHjaazS9sYSyPBh5SECF9sIW3qNlRGDt/HD/IhrGYBz0Qa8hYEpg9xn9TKkpGASDMFAJUGkl4ca3I0u7QGT+XuTHsPMg3wc7aGMFOceMgjkSiCDAjG/+ZBBXJwvENmjmsZBGpwR03FS2FxMVN6GtagPT+Xcllx3CVxqklguBQyC5Ip6nUBJaGAsgtv24UQ09ycxUsyM38zbNzollCsjxWsXnwB9IvVxpBvBJwa3+LkhvYEIShZtCPgRBE1NXFUT3JWAM2oBcSwMTpQSJjeCBeCWlhQfh6Po0saVMrljvRbZfFIEJYh5BIPF+T4WBp8UNPNQDBiHQM8oevmToYkizUnp7k1c87Jw6+5bWIwUhw1DLwU2gBEF8a0lG9DW5vfeZvEDiOsWRylbHAWDrz6lSq7A8lIcablmDmE2dixcYSPCqXs5hVJFGNyWHHcVpIAubVwpiy4MreGm7nwOuQ/F+b7DzPbsUk9K3IQeI8XM+M1w9xN8MXgTfCtPkd6cXBkp3nQxNVyiYIKfJy5V9AYqHouHH1BFak1aVQ2UGFWc+5vSjsWxmiglTG6kfjomIRyIu4OEA71eOjT1QPBw/7k8IY2elG6dhxPVC0g083MsrvHRBDWguXgEN3yYnwLTFWthtIcyUpwubMI3VQu0yekNNje/92lJIBWGF+et1aTPG0lpsTzZHI4IVxKD9UcoDo9wOv501mTySPx4a4okTu5NeodTM1iDPJQRrVyzyGbm5lLrfTgXeJBi9RkpZsYbcwpiYV8M1+m1wYf21H3IEPdjYwrux+ak4Mui+/EFcD82puB+bE5K+pt5c1Lcj425Pq5TY9bH/diY6+M6NWZ93I+NuT6uU2PWx/3YmOvjOjVmfU7Qj//4+qnBzf+AkP9OBxT/C0TGP2wnTvovfTY/18OBkyDD60kAf3LxLUsSSLV4wjeChRXLGLMqM3X6yX785x3amGtbS+ZBCeJZySD78Y+l10J68ad4SxKCPPr9AMvyxud9O1hYsYwxqzJTpx/rx+mAQzPe1j6ZB309tnPYiLgfxx0WvfhTvCUJeFeevVl5nVhYsYwxqzJTpx/rx4kX+7F+vxqBrYQ6SaXElfBr2UCAJAySpL/N9PzwCNdZnK+/uYURRwFqGKQG0MV1Uq/9iQMB4hQARhIpOZ4iRjgXqPP0FiBiHv1NifIokvIwKOrlKYgrNbwF+HYkiC9LbycljJuphEuBhRXLGLMqM3U6UhxZ55hLB+JT6GCN/Qa2IjpJaUjG8/2hYKBnXC6Nh8MjCEoAYCMCA0rtRpo3dSlcYQ/0yg+DY+NABUFaDIliuZgrxZt5qAGMR40Y5JFNTT0WpOF00/JgcFdpAxgUcJ1AmwZg9DaTCVdDKzfGLMtMnY4Ux9Q5zkRMpJPxWXhQ4gDF0cljVDZz6iSV8j4ux2UkwUDfjD8cTpfEoGy6IsabCQEjJMZ7ehq4BYN3SdQTiosTUBLYvTzU1AugSwZ55KY8Ed4qTkBpU87UX+u1AUZIjNMoN1aCT2SMWZmZOh0pjqzzdG7Ow4MSwwG/QuFKt9mPcb2Pa8RlJEFyZUzGm67AItEkYOAa/1CC9UMvMDBlSG5Pz7sELmSM48rmlOBiIthP3orfUJUEdi9P1BC4WEBx7kSNXhnRLRnlxnfGy0u7yrgM5Uzuw81fCiysWMaYVZmp05Hi4DrHdOmwnoEHJQfy5OWydbjrJJWBYDMuIwmSK2My3nSF1omreiGUEmugDGqie1N39HdtAS6mkyD1qodgFIGtJLB7eaKGwE2vOGq0FUS3ZJQbHaAhsDVkQz+GIVtxGXf5WmBhxTLGrMpMnY4Uu9Y5Dt+UH246rGfgQcmB+IaE85ffk3S46ySVcRvWistIguTKmIw33QjicTe0csJbGJgyyH1WL0NKCibBEAxUEkR6eajB3ejSFoM8dGFEzUOgjMtTEtLLGfUzm7kUfCJjzMrM1OlIsXed4+DDeUc7nYPz8KBkHjRj2LjCVkKdpDJuw1pxGUlQuzA4C8Aszw5P8G9c9det6gGwOYoubVwpkzvWa1dhcArGcWVQgphHMFic72Nh4C7jzTwUA8Yh0DMKanClGzWykybSW54MBGHoZcEGMKIgunETGKRLG1fqlwILK5YxZlVm6nSkOKDOceZiFqDTHMDWAfoQHpRsDPHQRIQ5FZRxG9aKy0iC2qVNUvvBtdbXboQZ+AgEfQgRwoG4O0g40PPPKEB9Lg4EegVAwQiXJ6ThpEpb5+FE9QIi1MR5lae34ERzeWkIbwEuphZEl89FGK83fymwsGIZY1Zlpk4/3I/NhVmwh6Gz8k9pV8J1asz6uB+bT7JCP+a3Z60Bzbj5Nf3UuE6NWR/3Y/NJFvl+rL8wB9drxgDPVSxjzKrM1Kn7sTHnxnVqzPq4HxtzfVynxqyP+7Ex18d1asz6uB8bc31cp8asj/vxE2z7x0eL/JMl85NxPzZmfdyPn8Cd1ZwU92Nj1sf9+Ancj81JcT82Zn1+Sj9mKxX8XZsMxl+mCPuhHtc/71AJcJduGsVsGtUTGLM3+LAVyxizKjN1eoV+HJtf+gXF+uWIioOBHlfckhIgAyK6yyDdGJeRBHSN2Q980opljFmVmTodKU5R52qBRH2RhjpijDNCYpwGghrIYDQSGtUTGLM3/uAZsz4zdTpSnKjO0T6xWqEGiSsFyR3r9XfUMgAMiuPfZsdRTYExe4OPXLGMMasyU6dX6Me3HvjVFNUgZVAT3Zt6qKcNg39ZDSNyG3wHdhwleBcU35g98SfNmPWZqdOR4hR1jn4Z/wt6/OKLBpk6pdxJPTTx/3KuwS2MTaMiFBTHmN0YfEqNMYswU6dX6MfqiOyOdGnjSpncST2bsf7mmaNoA4ll9ATFN2Y34gfPGLMmM3U6UpylzvV/3AI2QnTH1BGjO69HnC5gxxVUxlFNAeBXbdrGvB182IpljFmVmTq9Qj/eAzba4hizMP6gGrM+7sfbwTdaf6k1p8D92Jj1cT/eAr8Zx3/zZczKuB8bsz7ux8ZcH9epMevjfmzM9XGdGrM+7sfGXB/XqTHr435szPVxnRqzPqfpx/GneDfwxx9/YHgi/kixDIjpcuC7QELlfwvvTRiz4SrjvWtO7J2/yUcmBZiRk+J6/AIwY7GMMasyU6dL9GM21M2n2MNfuKFTUka5sSo7rTOm3WmKz/KRh4qTfmQBmLFYxphVmanTz/djtFL348RO64xpd5ris3zkoeKkH1kAZiyWMWZVZur0w/1Y59crp1ivHystjd+/f+NKpOcfBQgjFEMQg4J3BX+hJoO4UsNbgL9BU4L4CzU5ULcEFxYTJoEGivj4nEKaP+8wQ3p8BvWYYJwZUAwjLYmalJAurmM9F0ziAvResP4SulO/r0icFNRJHu4V7XoWZtbTMQjio+FukgHlj0owswmTYFSxjDGrMlOnI8UBdY6DDwcQTyVcS/RJ4iEbUdraoAADdQQjA1wY1CgeScPppvXD0GJgAxgUMD/geU0bho5gxZVQBgV005GtlQMsO+4G7Lg8GbhFu7kkETMDZY5r0EAm1L7FuQb65gIQVJ44NsbT2kicNIplx1GIwEaELmxOFAdKz8yKR+KktDVFfCgYT23CPBuGGGMOZqZOR4q96xynDw+4eKJtAMcZhkeQGXGlrQ3cxeEI+56gABdBanR0RnirOAGlTTnTsauc0gNGSIxHo9zuEDXYCkwaW4iMW64vQ8F6SXRJrceVrojxZkLACInxnp4G4yLtLYjDicYmsVwKGGzuVW8WCtJ0RJPKlky3ACMkxpt6upNgSLGMMasyU6cjxa51zqMn2s8eQ4IHa3ECSlsbuMtmmUAwamrwBwgq41clDdH5ThiXoZzJZQ8QiEdBc8YEJmXP4OwQ45qazS1pSBvt2hXa2/Ro4zWD5M7rqWQwggXcRn6HaxMxSVyt4rARH+xVb5aYIRFvJVlyn920STCkWMaYVZmp05Fi1zpvHnw4sMrtZ8BJms5lotOtNnBXR3AiagZAQ2BrSMrJuAzljC4M2YrLuMsLd+GN4gfQqgF2jz0ba2BD4rYom4wYvI2vXME8MJQNQCmxBsqgJro39bSeTYvBSO99RWKSKFYc9niverPEDIl4K8miC0O24jLu8tEsAzCkWMaYVZmp05HisDrfdgwJHaYJpa0N3OW5TzsSNQ+BEnk0hDnLva8pYKScUR9Pfy1JgnIjgDhkxfmCedhm4GI3YCuDssmIwdv44VMjzoXR5Vy0AW9hYMog91m9DMaFBha/hcYmMV3Zg73qzaLMxQ/EW0kml5MyCDSLBIwndxIMKZYxZlVm6nSkOKzOtx1DAifdZD9GEAYPYoCBOiV7R2SEmuL082sxsAGMKIiuJlWQLm1cezMWP3Ab+XWLo/RocRQM3E3B2o2wXeFK9+GaKZP7rB42Vq7pkg1ox7RiJgmARrK0V+A2xzOfCoBbUMJIMrkbNmEeDCmWMWZVZup0pDisztMxhNOw2V979PRKKwNBnMuw40kNl1AQxTU8WEVzCG8BrArXWhBdrocwjimioDljDfKoiwAotScxmx4/BpMmwQXgWvyJNYPoPqsHei8Yywip31dkMgnc3l6RepaUOcEHxDXJoksNYfzhJkyCIcUyxqzKTJ0u0Y+vCg7ceO4bsweuU2PWx/34aNCAsWn6foNmnL6ZGfN2XKfGrI/78Qfg31ETN2NzAPikFcsYsyozdep+bMy5cZ0asz7ux8ZcH9epMevjfmzM9XGdGrM+7sfGXB/XqTHr4378Nw9/snPbj35+nNWWrX/OtuxOYmFcG667rpP/2L44e3LMLMaYV5ip05/Sj80x4DODPlSc9di7B0fcj40xYqZO3Y/N2ziy223D/dgY8xFm6vT0/ZgnLM8+8uvrv45H+HVNBzGNqH8oGCe8LeL54ZE/vv++T/6VL4w4ClDDIDWALq5jfVxPXIB+MWT6Uen4CyNL6Dv1XDGibHV+yrR+Buvp6mWnKZpPARiph8OOw5UNVw7pLTUl4S3Q26IS+nrGEt2TY2YxxrzCTJ2OFKeocx6aOBzp8hBMLgwdvk39Q8FAz7hcGg+HRxCUAMBGBAaUagBp3tQzcIU90DcXgKDyxLExntYmBnNxMSDmkU2N4qA5XXPZyaU9P1xxrjDaMYnsXhIQ9ZoUQABXtvS7cswsxphXmKnTkeIUdc5DU72h6eKajGcFA30z/nA4XRKDsumKGG8mBIyQGO/paTAuoES8OHficJJGKZUMBFMeudQoYW+6JGu6uD41HFcZMZiSyKWgTpL0gLIUj/17V46ZxRjzCjN1OlKcos51So7d2nhWkFwZk/GmK/AFi9+r4pctwCNeYGDKkNx5fWoegl0koe98kfFcuBsfRLdkMN6bLsl67rPDZcTg5FLl9iYFdZ7i7MkxsxhjXmGmTkeKU9R579BMbm08K0iujMl40xVqCbjiWGcQSok1UAY10b2pp/WDfhybSg+MZR7YyikDwfl+3JwuyXrus8NlxOCGftycNMWpL86eHDOLMeYVZup0pDhFnfcOzeTWxrOC5MqYjDfdCOKxR6YmwVsYmDLIfVYvg3GhgcVvMTNXykMXRtQoLlckWc99driMGJxcak8vNJAkdz+OmcUY8wozdTpSnKLOe4dmcmvjWUHtwsBXIsbRn54dnvjz679pTzce+hxFlzaulMl9Vg8ba67/ERNtQLvZeybnSjlpJw24z5anS7KB++xwGJClIDI8XGrSA8bjbsBQHtiA9q4cM4sx5hVm6tT9eFZQu7TJhkaSYAa2CoKTHRHCgbg7SPisHqCp3NTf/8EzUBxIHJmZC9T5aw2op0uysfvUcK4c16R5uNSeHigISuj+ZzVcGYRYf3R7O5rFGLMsM3V6+n5szA/HdWrM+rgfG3N9XKfGrI/7sTHXx3VqzPq4HxtzfVynxqyP+7Ex18d1asz6uB8bc31cp8asj/uxMdfHdWrM+qzej/kDneKP8Cuf5sEo/cCo+PXrVx18CNYTf5D0vfBhmf8+TzEUHHCX77WwMZh3ZoWb4Q/pbp7ixeHXAI9fLGPMqszU6Sf78b++/9rFbcTfuiCQNv5ijRn2bjwiTjQz6YxmP/aeHcmffVORF4dfA2xCsYwxqzJTp5/sx9u+xSbqhsFIcabZu/GIONHMpDOa/dh19heT77q2E4FNKJYxZlVm6vST/RjNGC25OC+Q8sQ2jy/KeArCCFGcSp7sRKmSBlAGAeMMEog1EN/YcFff2zAccOzv37/vQ29AnxKC9G0vLezhLLTrlUdiEsDZYcS5ADUM4irjPujvOF3NCBgB1IjBow1WThkXCRgEaXhT1lwVKKH7knDFWA7HlYLk1nko4HASn65+ChiyAdZJF0ru8ytgomIZY1Zlpk5Hir3rHPl7J+ZT4ChEnuLczzgejjB0COLUkybGMS/F9REsjWxqFI/E/BDA1jkLG1Movwzcoq2BGAKXtoj6h7PQqFceiUmAkmAWZgBaiWaXQUF04yxK3tQrP0mamEc2NfVTgDi8ljVXBaIMNsBYDseV8eg281CgnNouEPUI8pFxlQBAk7biFWJmY8yazNTpSLFrnfNEw0EWXdobwFiensqTTkAAF0EKKI7EeBorl5rmScpbtHHapoNbBq4yFFTCeEvUesabs/RWHqkT4kpXxHgykqCeAi6CFJRQh5izt3Jq4DIeicOTLGUDvJvidDGWw3FlXG7SA7gIUgCDQellMB6RnhoG38J7sxlj9mCmTkeKg+tcB9YG8I2ErR1XfjuBgYQJBOsTlsSTFBr2OaJbMsqN72AIBuIux3KW1DJvKUKSaNcuScHxLL2VF/8LaDAEhgYSbo7AQGWQQaVcZLhJv8PkeBF09X0xEXNu2PN4K8l6qwIzs8jt5enpuYEMJjBKn8zehmyjN6MxZh1m6nSkOLjOMR2Os+I8iU5zXJkknbyid2LGE1bZiG7JKDe+g0MWYCxPW64EV6wErsbKiMHb+MolKTiepbfy4n8hmQYCKCXWwNqgWG5vnyNQkuJ/EXP2Vi6j3AjEW0nWW1WKa5SMFO/l6enxFDAYTFADAwkhY/At9GY0xqzDTJ2OFLvWeX1ywdUBtwEMjzlp1wl1dBb/ixhPY5W2N5ZAhqOWzRIujnLY0musjBi8je/kT8HxLFxqFMOlnUA83mVa2kB5QDIokJtmHAAZxMW5E3OmPHRhRE0i3kqy3qqUlkgGol5uL09PL4PxBHYY7wuC4r+Jtyc0xrydmTodKfau8/hFAUbsBxtAW2KjKv49v3LGgzXKog2B1tPUjE9bgLsScEYtII6FwYlSwuQK6QlcydIsoLnyGsQxUHeZhzm5DLq0cUUcBtrJXX7LrDhswLjy0GAQxDwiBZGkXnlzIIm3atl9UQ/ePoKKw5h/OiABiC7E9VMQNuMYeQvIWSxjzKrM1Okn+zHAFERHHoCtk3EenpW4Fv8OUjE/0OkJFI+HI+wYqTXpFK6BMj4IxHqQOFYTpYTJFdLL7c1Cmk+XqLeLsxCuBHfjkmgTDtdSNSNQkBqhuGDCGH9qz+Otpqy5KqA4G2TMQB4+HcUSJLd+CkJZ3HMoN3zUE8hZLGPMqszU6Yf7sTEfJPXRveF0xXkfrlNj1sf92JgRB/djfBV+/dtwjevUmPVxPzZmxGH9mBP9Ef5fhjfiOjVmfdyPjbk+rlNj1sf92Jjr4zo1Zn3cj425Pq5TY9bH/diY6+M6NWZ9fko//vX1S4Zf/Ic5h/3rHsCfcC1O4F3PYn4OzQ+SMWYpZur0Cv0Y64y/Y+EU9PrxGZ/FfJaz1KkxP5mZOj19Pz7yS+0bafbjkz6L+SynqFNjfjgzdXrufswGRga/G5IytkDy6+s/mUf4lVTtkEbUN7+zUib4qx7GY0vo6zc1luidmO31Z5FScHmcF7eoUdycF7zEYhljVmWmTs/djwG7jhpME2r02xjYk5ILQ6maetoRBNX8UpLmWBhwZSsuOBbX4rdo5u9NVy8PoNmz38swp0Zv1hizLDN1OlKcos7ne1j64phcXJNRC+iS2h2PRQQGgyA2SCFx8Vs08zefhRGiuOy0HnNe/B6NWZ+ZOh0pTlHnsdP0SJqeWxtJQDfCriZuIztj0YD1LRYwXpwv0tgmSTN26+Uxzj8NqIubU4NXWSxjzKrM1Kn78d9ubSQBXYGg4tLIoEbu8f0Yhuwk+/P7f8DRnBq8ymIZY1Zlpk7dj/92ayMJ6BJ8uYz9ld9EbyM7YylgECSXpLFNkqbn9pYnGwausO/3zYnh2zTGrMxMnY4Up6jz1JCaJE3PrY0koEtih6OALm1cKYsuDH0lval37seKMygXBtuwGrM5NX6JxqzPTJ2OFKeoc3Wg4rdImp5bG0lAV/BvfQk16HDjsTfpnRf/PZc0A7e5PATj92bY/CMCDCyJQXMu8GaLZYxZlZk6PX0/NuaH4zo1Zn3cj425Pq5TY9bnZ/Vj/pVsE/117lm40rOYvcGnoljGmFWZqdPr9GNjfiauU2PWx/3YmOvjOjVmfdyPjbk+rlNj1sf92Jjr4zo1Zn1O0I/1M7Lxh2KfQj9BG/n161cdfMhfd4rzbpAZj8n893mKoeCAu/wJ/WHsvR7+lPZbptj7l5988NXs+lzGmLcwU6ef7MdomfodFLC3teTmL9ZAKv4KqnkOO0/jRDOTPqs/kr3Xg+TPvscee/fjD3LV5zLmSszU6Uixa52no3zzyV4PZKQ402xewLPEiWYmfVZ/JLuu573J3Y+NMR9kpk5Hil3rHOfj5r+jTsTv2SD+ZTWmwFMQRojiVPLoJ0qVNIAyfiMHDBKINZBHv77YYTjg2N+/f9+H3oA+JQQaRdLCevrmwh4+OyihQHwQwJww4koANQziKuM+6O843eaM1IjBg2v/lUcRytKDRxgHehCQpo7PwndHEI/K5lZjSfeUJY4rjZgnPloJfQlK9DXelccYsx8zdTpS7FrnONFwqsZjq9x4HiTRsQhg8wSEoYMb00kT45iXYp2njEeNbGoUj8T8EMBGhC5sTKH8MnCLtgbylKctZvSMx4VBI7f37DEuUhA2IjCQnxsF0ry4yqAgus0Zm3rlJ0kT88imRvEEbnHxtIHs3rNweYonFwZltEFaCa406oEARnM9L/KuPMaY/Zip05Fi1zrHKYb8OklxPuoI2wBS4RyEwdMQBhKm9cNFkAKKIzGexsqlBi7jEd6ijQfBscvHUVz5ZSiohPGWmNGnuBYseLcXL84d5Yw2XRHjyUiC3owUlFCHmDPlkUtNegSShsSGygjRLDSUquneVaO4jFrQW8/rvCuPMWY/Zup0pNi1ztGJUwPGdDrFngXZ+M0DV/Z4nncJBNOxKHRuwoYmrk23ZJQb38EQDMRdjuUsmDE25luKkCTatUtiMAnkymC89+y9OEcJLJhBXOM+cOtEnFcGlXIHM/IPZEB/JkvEnBveSFo8lcUZPgsFPVcG4wAu33tTIHe8nld4Vx5jzH7M1OlIsWud43hKBzGmw7lWnCfReY0rk6TjT/AgLk4gHqPKRnRLRrnxHTwOwFg+F1eCK9uPxsqIwdv4yiUxmARyZTDee/ZePKHH1+IB8msKTVcbFMudmRFKUvwvYs4NbyRNTSVtGBqlDDKo6bkyGAdwsTzFZfCu3MF6XuRdeYwx+zFTpyPFrnWOIyzlh6tTbAMYHnPSrhPqfCz+FzGexiptbyyBDAcuWzJcnL+wpddYGTF4G9/JH4NJIFcG471n78VrINNTAz4abaA8IBkUyH12xuLciTlTHrowoiYhDZELY/wsjPdcGhgV47UA1yjAtbee13lXHmPMfszU6Uixd53jWNTXr5kvUmPQ/NgOi3/Pr5w6dhmXLNoQ6JxtatJRW4O7EnBGLSCOhcGJUsLkiod6GYwDzKup07M34wk8L25pB6KSc9GljSviMPQ2MYXizRlpMAhiHpGCSPLsG8EtDYENYMw8i27VLg3AeL0SGRREF4a2CDag/SLvymOM2Y+ZOh0pDqhzHGe3Yyn0LQBbx9Y8PGdxLf4d5Qc6IoHiOq8B7BipNemorYEyPgjEepA4VhOlhMkVD/UyGCcPnx2kIaLeTK6BcDrcjfPSJhyu5M0ZqRGKCyaM8Q1vhHqAF4Ergw+fBfRcGswG6pXISHG6tzF34npe5F15jDH7MVOnI4Xr3JhE6q+bwZ8A4h/dXsF1asz6uB8b82Y29+P67wz03fpFXKfGrI/7sTFv5pXvx/pbbvCuZgyQrVjGmFWZqVP3Y2POjevUmPVxPzbm+rhOjVkf92Njro/r1Jj1cT825vq4To1ZH/fjb/x1hwYejfYB6Kdd48+36Kdv9QPKAMHJf+bDRxDK/GLaV+DUxWmBu2nP99iZHwj2qljGmFWZqdOf0o8P7sECTUXNJv6ICxaDHpNWNb/hGBt7mHgx7Suwjxangq037v9OO/MD8eYYsz4zdTpSXKnO0/l+GJgU3YW2OhYXwyBaEQXoT4DBh0BZf198Pe0r6OlqECdx/+HusTM/EO2YMWZZZup0pDhFnfMEx3mNKxfMiMCtGJGLKzPg3OctdTgK2CSIOkdKrjiS1D0jTQQ4hHFG1HUUmQFL7U1He1va9HScgkFkYxDoqUEJfb2CEg1AjNUyCa4MJhfAhZJxRrY9wg/E+2PM+szU6UhxijrnCR6/LMJVw1CTiA0g2jj0NVY2BXAZbyaRG5tTArc4UEhPQwkxBaBmBozC8nAlJfqOtHocDGRmZqt3A8BQftiK12hJdPfbmR9I2kljzILM1OlIcYo658GtLqITn+hYlxGDqSvIpSDmpJ4GgzOk/AAuuwtvyYWBK/40wCAmuqk7cBnqUnFVL6Yt1h2mxZVGvRvp6TBjdBMaRXennfmZcIuMMSszU6cjxSnqPB30RGc6wd0okw2ZvvYBxWWkOGw1Bn6THtPsOggW5ws0HsBZ4KZVzfDGtFyzwHBmwJUCuUhe715xKlISzkKbwH3XI/w0BttujFmEmTodKU5R5+mgB3AV0V0ZMZjOd8VlpDhdgggpfot6INy663ANam/N6ca8Ky0E0kgvgxq5ykwYL05FLwldAHePnfkJDLbdGLMIM3U6UpyiztMZnVosXN6NMtm6G8VRwHhyI4jXXSQSBcofwfdsCua7Tp2n1m9L+3DrgNy0jHpVkZQEwOUKQXPshkf4mQy23RizCDN1OlKcos7TGc2TnS5vyYXB8z0OwUGvv3mWHQXRTW0jyZogoZpcnItgrO4yGwzMomAPCPgsoNZvS8un4+NwFF3auFIWXRh6opu6/4FJScBOO/MDGWy7MWYRZup0pDhFnTcPekQI77J1MY5rGoIj/iatft5JguiyaYmYBF/jaCe0HuUXGMW1EWSgciYtlaBuUZvTaqmAT408NDQ2uTfpHc7CYE0aRXbamZ8GtqVYxphVmanT0/fja4P2s0fX2Smt+QiuU2PWx/349OA74h6Nc6e05iO4To1Zn5/Sj//P//1f/O+z9k7Uf5H7FnZKaz6C+7Ex6/Mj+vHDTnmATcOYj+B+bMz6/Ih+bMwPx3VqzPq4HxtzfVynxqyP+7Ex18d1asz6LN2P//X9B3lJ/JHTeTAw/QMlJi/OCzR/avbt/Kp+uPbtvGtD9qC5tmd3/o/wX7sSe784gYkOm6sGT1osY8yqzNTpSHFknaOh4kgtzpPcjt7vvfxc/Tgtfg9O14+fpf71IPxTTnH25JgPyYBl36wxRszU6UhxWJ3zRN58omEsDt/Yzk/Uj485zX9gPwZIi+TF2Y1j3uCAZd+sMUbM1OlIcVid4zB95Sdieewig05kHfHprJRLgzKCsQwSnuOM8JsWSN/g4TIOGEl6BiMawueNM6YdmFwhSWkFg4BLKtHWykFMDmJ+QgGutTsY25wLlFC1NqLkNOI+1AsDmAV5ivMFxcyQXkq9Y72JIIi7ijzR5SiiBfRex05gomIZY1Zlpk5HimPqnCcgzrXiPw+G69hlhDlhMKjkcmng3GSch3VyYSQZDJ2w0YaeGup7p3BzOIfgyriYXCFopgUQQCa7qU8r1zLosiGJpgZXGTHOsc25AATNtQnlpKGB8cEjECghwbxUMoOWAeKqZPcm0meJQDDelmbyXYnLM8asyUydjhTH1DkOrBfPLKyTRyTOUB7KOkPTWSmXhg7WpntXtePKLyhL+kgaIlc5GRcpVdPFtZc2xXvdBcBFkAlLqIMmTS4NBiO9uVJca4uk5BiS4nQFOh/iCd5KGdLscgcTKc7g/f7fRGUv+a4cMIUx5kVm6nSkOKbOMYsOwW3EDLBxMuocjGdldGWkeHJpMEjgIjlbSAJBDSzqAEbpuxeQUka58UWK99xeWiymjsPorRy38Kciur0/Hil57TbH9uYCzbVFlFxGitMVSIi0xflOGjL5IqKLzHwuGZGo7CUv/j5gimIZY1Zlpk5HigPqHOfX67MgA/LQ5omptOlAlCsjxZNLg0ECF8lTRxEaWPxA76SWUW58keI9t5c2rZBxGL2VR6Akxf9CyZsuuY+7Abs3V29tESWXkeJ0BRIibXG+k4ZMvojo0oaBgRh+v/83UdlLXvx94NqMMSszU6cjxQF1PtMeHoJ1xlOSRzMXnw5EuTJSPLkyUhzTxbhI+kgaQhdGb0iK99xeWhkkxaUfABnExbmjSZtuhGN7czFenMolSi4jxekKvvTifCcNSavS7EmW3PihSkRlL/muHDCFMeZFZup0pDigzvFdtv4LwGfBOnHwFefrENTiYeikxqkKF8cloMF4z6WhPzHA0GphK65TmHpcGU+k4bR7Q1J84DbTAgiiDWhDA2hr5alzpOkEgvVmDsY254INo7k2oSQyUpyuwCyT/RhAXO9YkiWXzVijEriFp6PdTL4rmL1YxphVmanTkeKAOm8eo83ggHgaEhyCWjwPVqKWwCCuUVO7NHgWg3S2Yp2MAygRoZ52Ew1Rqt6QFB+7dVrCIOAjlGhr5YCbIxSPcF4S++tgbHMuUELV2ggnwlVGitMVmKX3mWkOefgimm76mAl+3pSq9zp2AhMVyxizKjN1OlK4zo0h7MfFWQzXqTHr435szHvAl+/e9++P4zo1Zn3cj415FX4z/uPlf3W4H65TY9bH/diY6+M6NWZ93I+NuT6uU2PWx/3YmOvjOjVmfdyPjbk+rlNj1sf9+FD+ulOc12imQgRvpI6/wn2e2YRQvncBb0/4M3GdGrM+M3Xqfvwe3thaDutSh03U4+MLuAauU2PWZ6ZO3Y/fwxtby2Fd6rCJenx8AdfAdWrM+szU6en7Mc/0X1+/0pLB+DsaGSG9X2TY++2J8ddA6ncl8pZAPEawkrQkurjG4XLT1LxL0i+g0EAa9doY17z6kVnGcY3u79+/76obaSJAjYj5lYe3AFciQb0w3RKcMSU028AeFssYsyozdTpSnKLOeabH/oo+JBfnvtpSjGOIWkWMy2ZajWWTU1wthC5SxTht5Yy3kttcUtILxWn01kYbpGfBlXG5MhgXKU43PSMMdXHYAAYF9cIADG14bzPNNrTJxphlmanTkeIUdc4zXWc9jLRs3u0d/UkvN6XVcBoMRiSQXY9NrgzGxcM4jd7aenFcx3HBeHEC0muXSOqvzQUwQmKcRrlhNsHNN8aszEydjhSnqPN0prM3JBBMLUQgru9zQNlkpDhsfOmEDfTVFkRBtAdub0lJLxSXMY4DuJglxeXKYDwyfkbsZ71pMpQwuXxegXgSmG1gD4tljFmVmTodKU5R5+lMT61CvLEfC0QI7ChI4p67cj8WEBDY0m/oxzBkKy7jLjcbwR4WyxizKjN1OlKcos7Tmc4mVx/xvaM/6dUjk743HCCeel4S91wZjIuHcRnNOBYT41EQ47jKYHwA00qvXSIPNw2C2L+pR1yCcsNsgptvjFmZmTodKU5R5/WZjqNfp7+Ofsb1t68P7ZRWbmpFUQYDd1OQwMV3StqYQnd7y1CqiNLKaMYB4ylhcwEw6ol6zygDQRhKCBvAiILoMiHjDNKljSv1ZhvYw2IZY1Zlpk5HilPUefNMZ8sh8ZbialSkjqe00WV3EdJgLFxc01jACInNCTSXpFTFv6O0MppxtElcQRzOWyQuoDkRaD4jk9AGvAU4Yy2ILicijGOKKMA+qMGbp8AeFssYsyozdXr6fmxEbG9Hgs6KbloccziuU2PWx/34Z3FYP66/4tffsM1huE6NWR/345/Fkd+P9bfiwM34s+AVFMsYsyozdep+bMy5cZ0asz7ux8ZcH9epMevjfmzM9XGdGrM+7sfGXB/XqTHr435szPVxnRqzPu7Hxlwf16kx6+N+bMz1cZ0asz7ux8ZcH9epMevzhn5sjDHGmNcpnbWP/2RtjDHGfB73Y2OMMebzuB8bY4wxn8f92BhjjPk87sfGGGPM53E/NsYYYz6P+7Exxhjzaf797/8Pt+VnNwi7nFcAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 03 - Index Calculation\n",
    "\n",
    "<i><br>“After estimating the modes from time-series, it is possible to calculate the distance from each mode to a predefined damping ratio to evaluate the stability condition of the system. Typically, a damping ratio below 5% is unacceptable in power systems because following a disturbance, oscillations in the signals will last for several seconds thereby affecting the system's dynamic performance. In the worst case (having negative damping increasing oscillations will appear indicating that the system is unstable (fragment taken out from paper [1]).” \n",
    "<br><br>“sssi” python script has the following configuration:<br>\n",
    "![table3.png](attachment:table3.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sssi(signal,ts,stepsize_min,variance_min,Freq,Damp,Nm):\n",
    "    \n",
    "    '''\n",
    "    SSSI\n",
    "    DESCRIPTION:\n",
    "    Small-signal-stability index (SSSI) function \n",
    "    INPUTS:\n",
    "    -       `signal`: Active power flow of relevant lines.\n",
    "    -           `ts`: Time vector with or without variable step size.\n",
    "    - `stepsize_min`: Minimum varable step for signal analysis.\n",
    "    - `variance_min`: Filter signals with variance lower than var_min.\n",
    "    -         `Freq`: [fmin,fmax] 1x2 vector with frequencies of interest in Hz.\n",
    "    -         `Damp`: [d0,d1,d2] 1x3 vector with damping ratios to compute index in percentage.\n",
    "    -           `Nm`: Number of modes used in Prony to reconstruc input signals.\n",
    "    OUTPUTS:\n",
    "    - `INDEX`: Three layer sss index with smi, ami and gmi.\n",
    "    \n",
    "    LAST MODIFICATION DATE:\n",
    "    12/13/2021 BY FERNANDO FACHINI\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    import numpy as np\n",
    "    import matplotlib.pyplot as plt\n",
    "    import math\n",
    "    import pandas as pd\n",
    "    \n",
    "    yh,x,l,yhd = data_filter(signal,ts,variance_min,stepsize_min)\n",
    "    \n",
    "    if yhd.size == 0:\n",
    "        out = np.array([])\n",
    "        sss_smi = np.array([])\n",
    "        sss_ami = np.array([])\n",
    "        sss_gmi = np.array([])\n",
    "        modelb = np.array([])\n",
    "        detail_l = np.array([])\n",
    "        detail_f = np.array([])\n",
    "        detail_d = np.array([])\n",
    "    else:\n",
    "        t = yhd[:,0:1]\n",
    "        y = yhd[:, (x.astype(int))+1]\n",
    "        y = np.reshape(y,(len(y),len(y[0])))\n",
    "        \n",
    "        plt.figure(2)\n",
    "        plt.plot(t,y)\n",
    "        \n",
    "        out_t = t\n",
    "        out_y = y \n",
    "        \n",
    "    \n",
    "        fmin = Freq[0]\n",
    "        fmax = Freq[1]\n",
    "        dmax = 0.25\n",
    "        \n",
    "        sig = y[:,:]\n",
    "        nL = len(sig[0])\n",
    "        \n",
    "        lamda, model_Poles, model_Res, model_K, model_that, model_yhat=pronyitesla(t,sig,Nm,(t[0])*np.ones([1,nL]),\n",
    "                                                                                   t[len(t)-1]*np.ones([1,nL]), 0.1, 1, t[0], float(t[len(t)-1]), 1)\n",
    "    \n",
    "        \n",
    "        dmodes0 = pd.DataFrame()\n",
    "        \n",
    "        Poles_sys = np.diag(lamda)   # number of modes\n",
    "        n_modes1 = len(lamda)        \n",
    "        fmodes0 = np.abs(np.imag(lamda))/(2*np.pi)     # frequency of modes in Hz\n",
    "        \n",
    "        for i in range(n_modes1):\n",
    "            dmodes0.loc[i,0] = -np.cos(math.atan2(np.imag(lamda[i,0]),np.real(lamda[i,0])))\n",
    "        \n",
    "        dmodes0 = dmodes0.to_numpy()      # damping of modes\n",
    "        \n",
    "        fmodes1 = np.array(sorted(fmodes0))       # sorted modes  \n",
    "        mode_idx1 = np.argsort(np.transpose(fmodes0))        # sorted modes\n",
    "        \n",
    "        dmodes1 = dmodes0[mode_idx1]   # sorted dampings\n",
    "        \n",
    "        jj = 0; kk = 0;\n",
    "        \n",
    "        fmodes2 = pd.DataFrame()\n",
    "        dmodes2 = pd.DataFrame()\n",
    "        mode_idx2 = pd.DataFrame()\n",
    "        \n",
    "        fmodes3 = pd.DataFrame()\n",
    "        dmodes3 = pd.DataFrame()\n",
    "        mode_idx3 = pd.DataFrame()\n",
    "        \n",
    "        for i in range(n_modes1):\n",
    "            if (fmodes1[i,0] >= fmin) and (fmodes1[i,0] <= fmax):\n",
    "                fmodes2.loc[jj,0] = fmodes1[i,0]\n",
    "                dmodes2.loc[jj,0] = dmodes1[0,i,0]\n",
    "                mode_idx2.loc[jj,0] = mode_idx1[0,i]\n",
    "                if dmodes2.iloc[jj,0] <= dmax:\n",
    "                    fmodes3.loc[kk,0] = fmodes2.iloc[jj,0]\n",
    "                    dmodes3.loc[kk,0] = dmodes2.iloc[jj,0]\n",
    "                    mode_idx3.loc[kk,0] = mode_idx2.iloc[jj,0]\n",
    "                    kk = kk + 1\n",
    "                jj = jj + 1\n",
    "                \n",
    "        fmodes2 = fmodes2.to_numpy()\n",
    "        dmodes2 = dmodes2.to_numpy()\n",
    "        mode_idx2 = mode_idx2.to_numpy()\n",
    "        \n",
    "        fmodes3 = fmodes3.to_numpy()\n",
    "        dmodes3 = dmodes3.to_numpy()\n",
    "        mode_idx3 = mode_idx3.astype(int).to_numpy()\n",
    "        \n",
    "        n_modes2 = len(fmodes3)\n",
    "        jj = np.array(range(0,n_modes2,2)).reshape(-1,1)\n",
    "        fmodes = fmodes3[jj]\n",
    "        dmodes = dmodes3[jj,0]\n",
    "        \n",
    "        print('\\n  f(Hz)        d(%)')\n",
    "        print('---------------------')\n",
    "        for i in range(int(n_modes2/2)):\n",
    "            print('\\n {:.4f}       {:.4f}'.format(float(fmodes[i]), float(dmodes[i]*100)))\n",
    "        print('\\n\\n')\n",
    "        \n",
    "        detail_l = lamda[mode_idx3]\n",
    "        detail_f = fmodes0[mode_idx3]\n",
    "        detail_d = dmodes0[mode_idx3]\n",
    "        \n",
    "        \n",
    "        th = pd.DataFrame()\n",
    "        \n",
    "        for i in range(len(Damp)):\n",
    "            th.loc[0,i] = math.acos(Damp[i]/100)    \n",
    "        th = th.to_numpy()\n",
    "        ths = np.pi - th\n",
    "        \n",
    "        thl = pd.DataFrame()\n",
    "        for i in range(len(dmodes)):\n",
    "            thl.loc[i,0] = math.acos(dmodes[i,0])\n",
    "        thl = thl.to_numpy()\n",
    "        thls = np.pi - thl\n",
    "        \n",
    "        jj = 0;  kk = len(Damp);\n",
    "        \n",
    "        SMI = pd.DataFrame()\n",
    "        \n",
    "        for i in range(int(n_modes2/2)):\n",
    "            for m in range(kk):\n",
    "                SMI.loc[i,m] = thls[i,0] - ths[0,m]\n",
    "                \n",
    "        SMI = SMI.to_numpy()\n",
    "        AMI = pd.DataFrame()\n",
    "        \n",
    "        for m in range(kk):\n",
    "            AMI.loc[0,m] = np.min(SMI[:,m])\n",
    "    \n",
    "        AMI = AMI.to_numpy()\n",
    "        \n",
    "        GMI = np.min(AMI)\n",
    "        \n",
    "        sss_smi = SMI\n",
    "        sss_ami = AMI\n",
    "        sss_gmi = GMI\n",
    "        \n",
    "    return  sss_smi, sss_ami, sss_gmi, out_t, out_y, model_Poles, model_Res, model_K, model_that, model_yhat, detail_d, detail_f, detail_l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 04 - Putting everything together\n",
    "\n",
    "<br><i> Now that all the functions have been discussed, we can run an example. You can play around with two data sets to test the output of the sssi script.\n",
    "<br>They are:\n",
    "* 7bus_fault\n",
    "* KYH_Nordic32_fault\n",
    "    \n",
    "An easy way of installing the PyEMD package in a specific environment, where your Jupyter Notebook script is present, is by utilizing the following two lines of code:\n",
    "<br>`import sys`\n",
    "<br>`!{sys.executable} -m pip install EMD-signal`\n",
    "    \n",
    "<b>Once you run the previous lines of code once and you have everything installed, you can comment the aforementioned lines out.<b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyEMD runs with numpy version 1.12 - 1.19.5.\n",
      "The version of numpy stated through np.__version__ command is: 1.17.1\n",
      "The location of the peaks are [24 25 80 81] and the length of the array is 4.\n",
      "The last element of the loc_peak array is:  81\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.io\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from utils.small_signal_processing import *\n",
    "\n",
    "print(\"PyEMD runs with numpy version 1.12 - 1.19.5.\")\n",
    "print(\"The version of numpy stated through np.__version__ command is: \" + np.__version__)\n",
    "%matplotlib inline\n",
    "\n",
    "#import sys\n",
    "#!{sys.executable} -m pip install EMD-signal\n",
    "\n",
    "#mat = scipy.io.loadmat('data/sssi_data/KTH_Nordic32_fault.mat')\n",
    "mat = scipy.io.loadmat('data/sssi_data/7bus_fault.mat')\n",
    "\n",
    "time = mat['ts']\n",
    "signal = mat['P_line']\n",
    "\n",
    "step_min = 0.04\n",
    "var_min = 0.1\n",
    "f = [0.1,2.5]\n",
    "d = [0,5,10]\n",
    "Nm = 10\n",
    "\n",
    "# ----------------------------------------------------------------------------\n",
    "# --------------------------- Original Signal Plot ---------------------------\n",
    "# ----------------------------------------------------------------------------\n",
    "\n",
    "plt.figure(1)\n",
    "plt.plot(time,signal)\n",
    "plt.xlabel('Time (sec)')\n",
    "plt.ylabel('Intensity')\n",
    "plt.title('Original Signal')\n",
    "\n",
    "\n",
    "sss_smi, sss_ami, sss_gmi, out_t, out_y, model_Poles, model_Res, model_K, model_that, model_yhat, detail_d, detail_f, detail_l = sssi(signal, time, step_min, var_min, f, d, Nm)\n",
    "\n",
    "print(\"Small Signal Stability SMI index: {}\" .format(sss_smi))\n",
    "print(\"Small Signal Stability AMI index: {}\" .format(sss_ami))\n",
    "print(\"Small Signal Stability GMI index: {}\" .format(sss_gmi))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "OCT",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}