{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":35332,"databundleVersionId":3723648,"sourceType":"competition"},{"sourceId":3910635,"sourceType":"datasetVersion","datasetId":2314156}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport datetime\n\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.metrics import roc_auc_score, precision_score, recall_score, auc\n\nimport lightgbm as lgb","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-04T09:23:50.958320Z","iopub.execute_input":"2023-11-04T09:23:50.958612Z","iopub.status.idle":"2023-11-04T09:23:55.526466Z","shell.execute_reply.started":"2023-11-04T09:23:50.958587Z","shell.execute_reply":"2023-11-04T09:23:55.525488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## A. Evaluation dataset:\n\nOne of best practices when solving kaggle competions (and ML industry problems in general) is to have a good evaluation set - ideally one that replicates the same conditions as the test set.\n\nThe CV scores should be tracked of course, but there's a risk of over fitting (CV scores used to optimize model). Like so, we should create an additional data partition - **validation set (20% of train data)**.\n","metadata":{}},{"cell_type":"code","source":"# Read the single statement data including the mentioned lag features: \ndf_train_lag_pq = pd.read_parquet(\"/kaggle/input/amex-fe/train_fe_plus_plus.parquet\").rename(\n    {'D_63_count':'ntot_statement'}, axis = 1\n)\nprint(df_train_lag_pq.info())\ndf_train_lag_pq.head(3)\n\n# Create validation set:\ndf_valid_lag = df_train_lag_pq.sample(frac=0.2, random_state = 42)\ndf_train_lag_pq = df_train_lag_pq[~df_train_lag_pq.customer_ID.isin(df_valid_lag.customer_ID.tolist())]","metadata":{"execution":{"iopub.status.busy":"2023-11-02T19:30:57.956943Z","iopub.execute_input":"2023-11-02T19:30:57.957920Z","iopub.status.idle":"2023-11-02T19:31:05.273567Z","shell.execute_reply.started":"2023-11-02T19:30:57.957872Z","shell.execute_reply":"2023-11-02T19:31:05.272296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## B. LGBM Baseline\n\nThe baseline will be build with a classic CV 5 fold with target stratification (using `StratifiedKFold`) and simple hyperparameters (no optimization). Baseline performance in the evaluation set will be our reference when trying to improve the model. Besides the AMEX metric, I'll track the AUC, PRGAUC (my implementation), and the `y_pred`, `y_true` distributions. \n\n### B.1 GPU accelaration:\nIt is possible to speed up lightgm training using kaggle's P100 GPU as explained in this [notebook](https://www.kaggle.com/competitions/amex-default-prediction/discussion/328606). To do so we only to turn on the GPU on the session settings and add the following lightgbm parameters: \n```python\n'device': 'gpu', 'gpu_platform_id': 0, 'gpu_device_id': 0\n```\n\n\n### B.2 AMEX Metric:\n\nAn implementation of the AMEX metric (used to compute the leaderboard score) can be found [here](https://www.kaggle.com/code/inversion/amex-competition-metric-python). Additionally, [ambrosm](https://www.kaggle.com/ambrosm) made a [discussion post](https://www.kaggle.com/competitions/amex-default-prediction/discussion/327464) visually explaining the AMEX metric:\n\n![amex_metric_explan.png](attachment:7c549a14-5f9b-4047-89e8-1b00215fc427.png)\n\nThe **AMEX metric** is $\\frac{G + D}{2}$, where $G$ - normalized Gini Coefficient, $D$ - default rate captured at 4%.\n* The _normalized Gini coefficient_ ($G$) is simply a stretched AUC: AUC is the light red area under the curve, which has a value between 0 and 1. The normalized Gini coefficient is equal to 2*AUC-1 and is between -1 and 1. The larger the red area, the better is the score;\n\n* The default rate captured at 4 % is the true positive rate (recall) for a threshold set at 4 % of the total (weighted) sample count. It corresponds to the y coordinate of the intersection between the green line and the red roc curve (marked with a green dot) and is always between 0 and 1. The higher the intersection point, the better is the score.\n\nIt's a metric heavily focused on making sure the all the defaults are associated with the highest estimated probabilities. Like so, model calibration isn't a requirement (e.g. estimated probabilities ~ real probabilities). The focus is discrimination, not calibration.\n\n### B.3 Other metrics - PRG AUC:\n\nThe **PRGAUC** is the area under the curve of the precision - recall graph. It's a [better alternative](https://towardsdatascience.com/imbalanced-data-stop-using-roc-auc-and-use-auprc-instead-46af4910a494) to the ROC AUC for unbalanced data, and also contains the same properties of normal AUC metrics (0.5 -> random classifier, between 0-1, 1 corresponds to the ideal classifier). Like ROC AUC, it returns higher values if the 0's are associated with lower probabilities and 1's with higher probabilities (doesn't take into consideration model calibration). See [this paper](https://research-information.bris.ac.uk/ws/portalfiles/portal/72164009/5867_precision_recall_gain_curves_pr_analysis_done_right.pdf) for more details. I'll be working with an implementation of the metric made by myself (based of this paper).","metadata":{},"attachments":{"7c549a14-5f9b-4047-89e8-1b00215fc427.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAeYAAAG5CAIAAADUKQC6AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAEOoAABDqAYLTCpgAAJQnSURBVHhe7Z0FWBT5/8dPPc+4s671Uj2vvPLUU884u7u7G5MURaQlpbu7u2EDMDAJFVEURVARkJKu/f/fy+yP4wbwPF0U1s/7eT08u/P97u7Ezms+M8zOvPV/FAqFQukkIWVTKBRKpwkpm0KhUDpNSNlvUAQCQXV1dX19PR6IBklcMHW1tbX4K3reDmloaMBH1NXViZ53+LzAPEF/RIK/J503pOw3KLm5uRMnTgwNDS0pKRENkrhwudzjx497enqKnrdD0tLSDh8+rKenV15eLhrUsYN5Ii0tHRwcLHr+HNHQ0FBTU7t9+7boOaXDhJTdUXL+/PkFCxZMmjRp8uTJs2fPPnToUEJCAiodUbM4UlNTk5KSUlhY2IkqxGcH0+Lr67tr1668vDxmSHFxcVZWVn5+PvO0PZKamrpt27YTJ06UlZWJBrV/7t69a25urqOjg4UoGvTcwTzJzMx88uSJ6HmLXLlyRVtbu/l27n5jKisrRc8pHSak7I6SqKioESNGqKioYM1xdnbesGGDgoLCpUuXRM0dLA0NDaJH7Rzsmz/jsx4/fmxmZoZdh5ycHNGg9s/LK/vZE9Vqrl+/fuTIkd27d1dVVYkGPV/wQf96fIPP52Ozd+rUKdFzSgcOKbujBMqGerDyoARGdaOlpbVu3Tp3d3c0odZGkWVvb4+B+vr6sbGxBQUFGI61EdYIDQ01NDTEnqypqSlKdQxHIXbjxg0bGxtNTU00xcXFMUdC8FdXVxfGQcHl4+MTHR3NlGz19fWPHj0yMTG5c+cOVu/c3NyIiAh8EF6OD83OzmaK/fj4eHwEPh1NKPdYlSxUcvXqVQiUGZNz584xRkM3Zpytra0x/lZWVpcvX2aOq+KzsOvt4uJy8uRJAwODyMhIVM0YjumytLQMCwtDE17i7++P4RhbPT09VVVV/I2JiUFZjU88e/bs3Llzv/zyy8OHD6MJ1WJ6enpQUNCFCxfwEffu3UNlGhIS4uDgwEwLWoXj2ji2eBNIChOCDSTmoaOjY8sCFkMwT5g5z4wJ5gZm4NatW2VlZe3s7DDmeJPTp08zB0nwFx+HaVFTU8NffASzpPBxHA7HqzFGRkZubm6YV3gfzC70xLh5eHhgETMex9xGhYvxwZtjZmITji03OowfP37UqFHHjh3DJ6IzvidYyhh5jBgWK+ZeaWkpXo6tFxafk5MT9j/wzpiN+FZgcWMHC62Yb5gKvC1mF8bw4sWLN2/exAIaOXLkrFmzUDFgXqEPXoV5wuy74NuIryW+Hurq6uiJr0FFRQWGU15LSNkdJSxlQyWrV6/Gigf1YE8fNd2+xqDOgp4CAgJgk6dPn2JV3L59+549e/bv3w+JwFZY4TMyMuTk5Jj+O3fuRLWOdRg6wMr81VdfQQHwzoEDB9AfUsNH4+P8/PxmzJgBz0Ix3t7eeDdESkpq06ZN0ApTw2KUvvvuO7xKRkYGf2H2xhEXBh8KG+I9MTJ79+7dsmWLoqIij8fDh8IIQ4YMwU7DwYMH8YbYDkEWcA2mCzUyHIEPwniiypOWlsZH4yWYA2PHjl2+fDlegibICyMM0zFPMUqYKEwptkCYXZMmTRo4cCAG4n1gcIgGcwMbBswfyOXnn3/esWMHRgYjtnbtWviR8TJeiP4bN27EX4zw4sWLISxW1YzRgBmPHj2Kj0Mf9MSmiFEtpmLOnDmYsZjkJUuWoP6F+/ASbFrgtUOHDmE8Mcl4VWBgIKaouLgYSxCjiiGYTGxIYPnw8HCmJ6Z95cqV2DBgw4nO8DUsjK0COqOV0TFmwu+///7TTz/hC3D8+PFbt25hTOBuvBbjgBm7efNmzHB8JTAmeOFff/0FueOrgsWNLd/69eux+ccEYtuJkUd/vDNeiHHAdg4GHz58+IQJEzAQ7/nw4UOMGFrT0tLwEsxGjAzmKl6FgfjKYXKYWUR59SFld5RA2aihUBRjpUJ9tGrVKpiCKZAhaKgHZRrWn8zMTMZ98DjW23HjxsEXeFBdXY01DW5FVQvRo2g6c+YMij6s2MzKidWsSdlFRUWQGhwE8eGjIRooD25CRQmXoTMew8h4OYo7rPxcLheGhbIHDx4M4cIL+Ljmu9so211dXX/44Qes3mhNSEhgDuxg5KFsvArbA9R6GAeMG0QJ/2JaoDNspVCWohsmASbCVOOtGGVPmTIFo4fKEU/xQmxOmKbk5OQVK1ZgmwHBtTwwwlI2thbYwEA9Dx48gPswuxgtwkFwHzpg2uGgP/74o6WyMZfQ9O233wYHB2MEMDfwWsxeKBuGxRJJTEzEQAgRE4siFC9B+YnFhM0ePh3Ti7mKT8Ek4OVQNuYD9h4wzpgKfC4mGRtXPMZ7QpGQL7OIMQlQM+Y8+mBriklDtcs6MIJXoUaGpvHdwIdi0WMc0AEbTigbnwv/MksNY4I5zCgbWz6M5/z589Ef74BFjHdmimjWgZEmZWM3CF+2RYsWYbbjo9Efr8JLRP0orzyk7I4SKPvrr7/+5ptv4ALs6aOkwi4qjIC1GuUSXHbt2jUUfYi2tjYqR+gMVdW7776LDlgzRe/SeD4DWrG+QdZMf8gC6zOKqebKhlBQgqHIZY664KMZP2K9RUHK7Kcj2EKg/oJn0QRlT548OSkpSfRJzYKqGYLGC5ldZtgcZocWsTMOZQ8bNszY2Jg5kAJPoeaFtaE/9IcOUBozn4UaFjaBmmETKBsVMT698e2Fh26gdYw/dgvQE1OHQvLSpUv/qmzsFqBEZbYuMNf48eMvXLgA9fz666+2trZ4TwzHyKPybalsjCo2XZh12AiJBjUGyoYrMebMgQhMI0ZGSUkJj5mCumk8sV3BgsMCYpSNwhzia3wPYU/MK8wEpidUjveE/fEUm5bZs2czewNNYSkbWy8sZcxhjCRejr8WFhYLFizAxhLKxhzAvGVeiGCeMMqGbe3s7PDmYWFheAkjX8ycZygbmyVUEijG6WBIBwkpu6Ok+YER6BW76pqamtAczLt69epevXoNGDDg/cb069cP6oRx/Pz8Bg0axBxwbAqUilaovHl/+Ag1V3NlYw1EVbtmzRpU1v7+/qhGsebj5crKyp9//jlewrwWee+991AUY/WGsrG2o1JjPqh5YCJGW031F17C7CVA2ShUvby8oC0Mh2QhL1RteB+MVf/GMB+ED0W1C5swyjY3N4fRmHfDA0h/xIgRH3/8MXpibuAd8Ob/quzRo0djw8Y0wTvYY8BT1JjYIcCuDD4IwzHhWlpaLZWNOQ8vYwODDadoUGOg7Ob/foTX8Iny8vJ4jLoYG1Q47pNPPsH8x3hC01iajLJR/KJz43v8H3QfERExc+ZMLEH0xEz+/vvvra2tsf3AcsH7M/8/aApe2FzZEC42Bnh/ZtYhmI3YPqFOh7KxqZCTk2NeiDQpG4+xc6auro6dsE8//XTZsmX4vmGH6RnKZvZUUC6wNiGU1xVSdkdJc2Vj9bCxsYHy3NzcmCobJRjcBzExgapgB3SGmrHr2rLKRs2LfWRR78Y9a6zqeNCkbNRWrq6uqMXgVqyf+/btY6pgrLfwAna3mRcicCvKTObACATUqrIhGggFomleZW/ZsqVllY3JQekKecHC6A+5oFgWfVJODlP3McpG2cgoG++JPYCffvoJZTJejhHA/gEKTEw+5gPM/p+UDdE/Z5XNjCpk17LKblXZeLfg4GBMLDaBUCrGE0sBm97Lly8zysZCbFI2xm3t2rV4fyym+/fvY3GvXLkSo91Wlc1SdkFBAZYyPAsFM7MOwaYIrYyymU0Ik+bKxrzFeGLGYoOEqcCcRB2NefKMKhvzENtRqrI7SEjZHSXNlY2nKP1Q4mG1wSodGBiIAhNlMtyB4gsugE2gXfz9888/4QVokXUsG8o7c+YMVjOs+XgrpgNam5SNj7hw4QJei5J8zJgx+HRmncQIYHVFoYe1GuZFgYmVFpsHND1D2cyxbFR5MBFzLBuOwJujlsRHDx48eNq0acyxbEdHx6Zj2dgwTJkyhTmWjenC6KGag3RYykYr3hxvwrTifRYuXMhUiPhcFxcX1KfMuS7o/DzKxvvDVtAfOmDSfHx80K2lsjG22I/55ptvMIZ4jFmBOYzZ25ayMTLOzs6fffYZHI3OePOlS5fOmzevVWVHRkZiZurp6WE88bZY1lj6GG28EH4cPnw461g2NszYCi5ZsoRZTJgEGBbbLWzM0AcbVMyr69evY149Q9lY7vg6YV7h5VhM+FDsG2E8sYj379+PTYLoBf88lo2mBQsWYOuLmY8xwQc17UtRXn1I2R0lLGVDYdjlRzWKQhtVmIaGBlYhrGCwA1YnKAzuwFqHgg6VMlNwoVTECowXYp2Ec7GmoTNegkoKmsaqzlI2FGNqavr5559D2Xg3plSHJmBSvBVeCKlhrUYliG0Dmp6hbHwo1Hzw4EFU1hgZjLaioiKmpeF/Z4xgCMYH+w1r1qxRU1ODOPBxqJFPnjyJV+GzEEyXtbU1nMVSNqQGp8BW8A5GSUlJae7cufggvD+2Q0lJSRMmTMBTjCq2Us+jbHw0/mJ81q1bx0wm3hxuZSkb3VDz4uMgaEwU3hb7CihO21I2jIY9hpkzZ2I8MQSLYP78+cuXL29V2XiAd0YHzBMYFm+IyhqjjTmGOYxZjc5owvvgMd4WiwzLBeOPBY0lgrkKQWPxYaoRTALmHr4qmGPPUDa2DV5eXngH5guzatUqbEEfPHiALxh2VmbMmIG5YWBggG5NysY0nj59mpkDGB/sjdEZI683pOyOEpQz9vb2cARWWmYIVhispcxvILH/C03r6+tra2tbWloyJwIzVTBMZGJiguFY4bG64oXoD8mi1kYRhxUeO91QHgYy/11ErcpUSfgglE54Idbbpg9FUEnFxsZCT1qNp4FjFYVbMRxCxDgw/3NrGQgLOsO44Q1hW9TCjNEgl2HDhmEIgl1+Ozu7K1euMJsljD/2ALAJgSaY6YJeoVqMDMYZ09LkUAgLkwm/o5uvry+8g6eYV3gHbLcgI0wpmjA5t27dgqHgOAgXMw3bAHwE8yaYJw4ODswcxhzgcDiYb5ghCDYb0DdmJtOzKZhpmZmZmD+YjRh5bBGx2cPcCA4OxgYDo4o+2NrhE2NiYvAYI4NdIvREfzzAeKJOhwExc9AfQ9C58Y2FZ3BjPmAmY2OMkcTkeHt7Y7TRhA+FtbH3wCw+TC8mAZODj7a1tcVkYluOScM8xMRi7qEb5gxmL5Yylg4UHB0dzYwPE3TDph1bGrRi1wpTjQ/F+0DxzDtjy4cvG94W78NUA9gJiIiIYMYW8wobPzgdrRhhfA2oyn6NIWVT2jdQ9rfffgtbteX61xJsDFAqQlWwD7ZPa9euhcUgVlEzhdJR8yqUjZIBe68omlAaiAY1CwqoGzduoL5AFYAipXm5R5GAdExlY58GhW1UVFRYWNjx48fnzZvHnFwoaqZQOmpehbINDQ2nTp36448/MueuNg92LbGzuXr16gkTJsyfPx97c/SPaQnLnTt3/vrrL8ix6ShHRwg2JLt37/71119/+eWXlStXenl5Mf+9pFA6eF6FspOTk1HCKCoqtlT21atXVVRUdHR0mIOV48ePb/W/W5TOG6iwvuNdoZsZK5TVCB7Qvh2ls+RVKLu6uvrRo0dqamotlY3qRlNTk9d4MYqsrKyffvopKSmJ+a8OhUKhUFh5Rf9+LCgo0NDQaKlsMzMzLS2t69ev43FeXt7EiROZyywwrRQKhUJpntesbCMjIyj71q1beJyfnz9jxoyma0g2BUX3mTNnDA0NUY8rKCgcOXIED5pHRklmj8KeXfK7Dh49qKquKhpKeTVRV9dQUVFXVlZXUlI/dkxNUVFNQUFVXl5VTk5FWlrlwIET+/eD47t2Ke/de0JK6sTevcp79ijv2oUhBNFBOLZr99Fde47t3n1CSamD/xf6NSvbysoKyk7535V8x44dGxcXxzo9tra29tKlS7a2tqampgsXLly+fDkeNM+k3ZPeXfTuW3Pe+nb9t8rayqKhlP8YExMTIz29U6qq+keO6B46pL17t+bGjarz5ilPmnRi6tQT06efmDJFeeJEpbFjj/3+u+LPPysMHSr/xRdyAwfKfvKJ7EcfyXzwgcz770v37y/dr590377Sffoc7t37UI8eh955Bxx4662D3bodfPtt4d+uXQ926YIhrxfhaDSOD/MYY3v43Xel33tPSJ8+Mn37CunXT6Z/f9kBA2QGDDjUvbv0u+/KffSR3Mcfy33yiTz49FP5gQPlBw1SAJ99pvD553h6uGfPI198ceTrrxUHD1YcMgQcHTr06DffCBk2DBz79ttj330n5Pvvj/3wg9KPPyoNHy7ip5+O//yzkF9+EfHrr0D5t9+EjBgh5PffTzTCdDsxapTK6NFC/vhDxJgxIsaOFTFunGoTf/4pZPz4JtQmTBAxceLfTJqkDv76S33yZBFTpjBoTJ36N9Om/YPp05vQBDNmiJg58x/MmsWgBWbPFjFnzj+YPfvk3Lkn5837b8yf/y8sWNCSgxPn/PXF9AmfT981fvKY336rrq4WqadD5jUrOzg4WF1dPSQkBF5OT0//5Zdfrl69+oytHHMyv+jJ/7I/Yv+Huh++pfLWVOepN5/cFA2ltIxAIKitrSspqXn8uOr+/Yr09LKkpNKzZ4t5vMLIyHwfn1wHhxwjoyx19UwFhYx9+9I3bkydPv3yiBFJY8Yk/fnnlT/+uPzbbxd/+OH811+fGzjw9LvvxnfrFvfWWy9MfNeu+JvQowfeKuGddxK6dz/bt+/ZAQMwJPHjj88PHHh+0KALn32GB3h68csvLw0efGnIkEtDh1765hvhgyFDLn/77ZUffkj68cek4cOTf/45+ddfU377LWXEiJTff08dOTJ11KjU0aNT//jj6pgx1/7889r48dcnTLg+cWLa5MlpU6bcmD49fdas9BkzmMcZixZlLF58e+lSIcuX31mxQsiqVZmrV2euWYO/N+fMyViy5N6mTfc2b763ZUvW1q33t227v337/Z077+/alQ1278bTW/Pn5+zd++DgwQeHDj0Ehw8/lJZ+JCsrRE4uV14+V0Eh98iRXEXF3KNHHx87lqeklHf8eJ6yMsg/cSJfRUWIqmqBmpoQdXWsOUI0NZ9oaT05eVKItnahtjbTrVBPr0hfv8jAQMipUyIMDYuBkZEIY+MSExMRpqZCzMxEmJuXAgsLIZaW4KmVlQhr66c2NmV2dn9jby/EwaG8CUfHcienv3F2FuHiUsHg6irCza2yCXd3IR4ef+PpCaq8vP7G07Pax6fa1/cf+Pn9Z/z92QQENCfH3stur/WOmSY6W2z5WgYzJ0yoIWU/ePDg9OnTUlJSO3bsuHHjRn5+PgpqpKamJiMjQ1dXV1lZ+dy5c/b29itXrmR+o9xWWlX2ochDH+t9DGVPdpqcli/6QfCbHoGgoaqq9smTmkePKjMzy1NTn164UJqYCDs/dnd/YGaWpaZ2e//+tFWrUqZNuzxy5PkhQ+K7d2/u0+cH4j7ds+fp3r1h3jN9+pwbMODcBx8kfvhh4kcfMQjN+/nnF7/6SujcoUMvDxvGqBaSxd9r48ZBmkK3jhp1c+5cqBN6vbt+vdCJO3bAg3iQuW7dg/37hdaD8hQVYTqI75GCArQl9JeuLoQFN8E+cE2ZrS3MApsIldEoiCpv75qAgNrg4LqwsPrISAGXK+Dz/y8ujiAidPy3z7PW3e3yOCDmqoMDKVsYPT29KVOmDB06dMiQIQsWLHBzc3NqTG5uLsq+2NjYTZs2jRkzZuHChVeuXHn2XkmrypaJlvlE/xMo+y/Hv64+vioa+iZFeCfBmpqGysr6p09rCwtrCwpqHj5E+fzQyurusWPXly8/P3RoXJcuLNW2BYrf+LffRs0rrHx79MBfxshnYOT33hPSp4+wHO7X72z//gD1b8rIkcIa9q+/0mfPvrtuXda2bVBtzr59IHvPHqYYRB2H+ktoz8DAOniTpEm8PgT8uOoYXi2H76zoFXbSPz8oBgNJ2aLU1tZWVlZWNKaq8TptTJhzdevr6zEQHfD3X0+PbVXZcjFyA/UHQtmTHCcl54ouIf9GpfrRoyIe75GtbcaePVdGjbo4fPjpPn0SevWCcOPfeQf+jWs8BPGvCI9UdOmCujhlxIhrEybAvxlLlqTPnHln9eocKSnszkO+qGqLTUzKbGxQvdYEBdVFRDRERTVERzfExIiIjRVwOP8AVS2PJwSaZmixChHEq6Qymheg7nvJOqQolFPH5cPgGEjKFn9aVfYRzpHPDD6Dsic6TLz88LJoqARHIEAdXRgTc0dGJmn8eKGjmePLgwadaTwQDE2zXNyccwMGJH788YXPP08dPTpj2TJUxA8PH85TVi42MiqzsxMeT3B1rfLyEh5J9POrCQgQ/g0MrA0JqQsPr4+IqI+MFAoaXqbDC0TnpDA01lPZe+S3BrH6AWWR3KbhpGzxp1VlH+Me+/zU51D2ePvx53OEtxiXvAjq68vT0h5aWNzcvPn6okVX58y5MnIkHB3fvbvwf4CtHfQ48957wn/ZDR2aOmrUrfnz723alKuoWKClVaSvX3zqVLGx8VNUym5uULPwSEVYGCxM9S8h8eQFxrgpec8fZ2ot4/HQLxoldlMTKVv8aVXZyjzlLw2/hLL/tP/zbPZZ0VCJSH1ZWfnVq7lOTnePHk1bseLSzz+ffu+9uK5dmXMtmnO6d+9z779/5dtvb0yffnvp0uxdux7JyuarqBSoqxedOoXyGYVzfVRU07eTIN5McgNiAjX8Tkm55gfF1PP+UaOQssWfVpWtEqfytdHXUPY4u3EJWQmioZ02goaGqnv3imJiHru5Zevr39y27fKIEaf79mWV0qd79bo0ZEjK779fnzABRfTdDRug6cdHjz61tKx0d2+IjqaSmSCa88g/5opNaLJt2APfaNTarFZAyhZ/WlW2erz6YOPBUPZYu7H8e3zR0M6WhhrhHaWK4+Ly/f3vKSunTJt2fsgQoambafrMe++dHzjwyg8/XB0z5ubcuQ8PHy7U1X1qY1MbGirg8VjfP4IgmsgPinVW9NLa4RyrH8BqaoKULf60qmytBK2hxkOh7D9s/+BkckRDO0kEdXXVDx+WX71aHB+PmvrSzz+zfp+Cp2f69oWpUVPfmD79/s6dxYaGtSEhrG8bQRBtURnNc1fyXjPNQneX8PxrVmsTpGzxp1Vla5/RHmYyDMoebTM66naUaGgHj0DQUF1dX15eeffu7f37z374IfuXLF26xL/99umePS98/nnGkiVPtLQaoqLoWAdB/FcE/LgU+7CFf5qa7Hdjzr9uC1K2+NOqsvXO6X1n+h2UPcpmVERGhGhoR45AUFdamqOvf3XOnMQvvkjo2bPledPnP/kkfdasAjW1moCA+qgoAYdDviaI/0oDP+5pBKc6hpfpEfkkJBZPWR2aQ8oWf1pV9qnEU9+bfQ9lj7QeGXIrRDS0Y0YgqLp374GJSerMmYmffXb6vfeaHwY53bv31bFjc/btKzY2LndwqPL2ro+IIFMTxItRGBLrq+KzdJLZI//oWg6/4Z/nh7SElC3+tKps4/PGP5r/CGWPsB4RmB4oGtrBIqitfXr5crae3vUlSy5+911Cr15Npj7Tp0/qyJFZ27YVaGo+tbKq8fcX1tQtvk8EQTw/jwNj3I55L51oZifnURr+XCsUKVv8aVXZphdMh1sMh7J/s/rNL81PNLQjper+/Ue2tmlr1lz49tuEHj2aZH3lhx8yFi3KkZIq0tOr9vNjfYEIgngx6rh8vlHQ3sU2+ntcC4JiWedftwUpW/xpVdkWFy1+tvgZyv7V8lfva96ioR0jgrq6yoyMeyoql37+WXjMmjkA0qsXZJ0+a9bjo0eF51DHPOtfIgRB/CdKwjnFYZyLViG2cp6tnn/dFqRs8adVZVtftoasoexfLH/xuOohGvraIxA0VFaWp6ZmHjly9oMPYOr4rl3P9u9/+bvvbsyYkaesXBsSQidTE4R4KQiODdb0O2sWXB3D+9eD1yxI2eJPq8q2u2L3m9VvUPZPFj+5pLqIhr7WCH8Xk5dXzOff3LaNOXsvvlu3xI8/vjlnTtGpU3SomiDag6cRXI/j3iunmBvuc0OhzWr9V0jZ4k+rynZMdvzd+ncoe7j5cDwWDX19EdTVlael3VNWTvz8c+ZIiPCKph9/nCMlVRMYyPqWEAQhFgT8uAB130m/Gv7r+ddtQcoWf1pVtnOK8yibUVD2j+Y/2ifZi4a+pjRUVeU6OCRPnHi6T5+mE/gufvll3vHjtcHBdCSEINqJBh5fb7drkIbfI//o/3pIhIGULf60qmzXVNfRNqOh7B/MfrC+bC0a+jpSV1qapa5++bffEnr3Zorr071735w7t+jUKTpyTRDtRGFIbJxx0E2XiHteUSXhwlsWsDo8J6Rs8adVZXtc9RhjOwbK/t7se4uLFqKhrzaC2lrhfxrl5S9+/z1zGp/wbOvRox/s319mby+8rl6L7wdBEC/P48AYdyXvTbOsUu3DajkvKGsGUrb406qyva97j7MbB2V/a/qt6QVT0dBXmNqiooLAwBtr1wrvOfD22/D1+U8+uTV/fqGurvD6TfTzRYJoHx4HCH8vs32etc4ul/ygmBc7HtIEKVv8aVXZfml+f9r/CWUPMxlmlGgkGvpKIry2dePPZFKmThVdz7pLl4tffpm5evVTa2vWF4IgCPGS4RZpfshdY5vzM67P9/yQssWfVpUdmB44wWEClP2NyTcG5wxEQ19BBILq7OwcQ8MLw4YJZf3WWwk9e14eNix7z55KDw/Wt4EgCDFSGBJ7zyvqtntkpkdkrjh8DUjZ4k+ryg6+GTzJcRKUPdR4qM4ZHdHQ9k9daekDE5ML33zDFNene/dOGTGiSF+f7tdFEO1KeRTX+4SP1g7nBJMgVtPLQMoWf1pVdtitsMmOk6HsIUZDtBK0REPbOQKBIM/LK+nPP5n6+nSvXhkLF1b7+dFNxwmiXann8QPUfcf+eMrkwAuef90WpGzxp1VlR2ZETnGaAmV/bfS1ery6aGg7J9/fP2XyZObkkHPvv39n5craoCA6jY8g2pWqGN5Fy5Bxw0/Zyno88I1+zus9PSekbPGnVWXH3ImZ5jwNyv7K8CuVOBXR0PZM6blzV+fMOd2nD3x9tm/fG9OnV7q7sxY/QRDipZ7LL4vkFoXEhmv7P/CLfuHzr9uClC3+tKpszl3ODJcZUPaXhl8e5x0XDW2fCGprK27cSFu58syAAfD1mT59bkydWmpmRmfyEUS78jgwJkjDT2+3Sz2XVxHFfcnz+VqFlC3+tKps/j3+TNeZUPbnpz4/yj0qGtoOaaipqUhPv7Vnz5n+/eO6dEno0ePa+PGFOjoCLpe17AmCECO5ATGux7x3L7Qx3u/WHrJmIGWLP60qOyErYbbbbCj7s1OfycfKi4aKO4K6uoqbN+/Iy59+9134Ov7tt5N//jnv+HH6ZSNBtCv1PH64tv/+pbZqW53EdT5fq5CyxZ9WlX3m/pm5bnOh7EEGg2SiZURDxRuBoCorK0tLi/lxI5R9aciQXAUF4Y8bWyx4giDERXkUtzyS63rMy+SAW65/+5ZHpGzxp1VlJ+YkznOfB2UP1B94KOqQaKhYU/f06QNz89P9+gl9/dZb5wYMEPo6OJi1yAmCECOV0bw4o6D73lF1z3Gz3ZeHlC3+tKrsiw8uLvBYAGV/ov/J/oj9oqFiTUFgYPKECYyvwUNp6dogcZ7DTxAEizouP1DDd9zwU44Kns95v92XhJQt/rSq7CuPrizyXARlf6z38d7wvaKh4kvFzZu3du5k7ol+unfvjCVLqn186BRsgmg/nkZw/dWEv5exkRGef/0KSmxAyhZ/WlV2cm7yYq/FUPZHuh/tCt0lGiqm1OTmZsrLM/eXie/e/fI331S4udGdwAiiXckPitXe6Wwv5yk8KiLu86/bgpQt/rSq7KuPry71Xgplf6j74bbgbaKh4oigri5bT+/i99/Hde0a363bpcGDcxUUWIuZIAgxkhcYc9Uh7I575CWrkPygmFfma0DKFn9aVXZaftpyn+VQ9gc6H2wO2iwa+vIRCIq43CtjxsS/8w5K7AtffHF/+/a6sDDWYiYIQlw8DozxOO6tud053Tmc1fQKIGWLP60q+2bBzRU+K6DsAToDNgRsEA19uQjq62vz89NWrRL+auatt8598MHt5curvL1Zy5ggCHFREBwLX+9ZZKO62am9z+drFVK2+NOqsm8/ub3KdxWU3V+7/1r/taKhL5f6srI8H5+zH34IX4PrkyaVmpuzFjBBEGLkklWI9Eo74e9lXoevASlb/GlV2ZlFmWv81kDZ/bT7rfJbJRr6EhE0NFTeuZP0558JjYdEEnr2vL9jB/0qnSDaiepY4WVDsryibjiHi+X+Mi8GKVv8aVXZ94rvobiGsvue7Lvce7lo6EukJjc3x9Awvls3psS+u3FjlY8Pa+kSBCEWajn8IA0/1S1O1xzCG3h8AZ/d4ZVByhZ/WlV2dmn2+oD1UHYfrT5LvJaIhr5oBPX1xfHxF7/7jvF18i+/FOnrU4lNEO1BaTgnQN33z+GnPJS8C4JjWa2vGFK2+NOqsh8+fbgxcCOU/a7Wuws9F4qGvmjK09JuHzwY3707fB3fteuDgwdrAgNZi5YgiJenLJLLMwycM8bE7KDbqzz/ui1I2eJPq8p+XPZ4U+AmKLu3Zu+57nNFQ18o9RUVj11cRLff7dIldfToMltb+qEjQYgdCBoldppTuNE+t0f+4r9fwQtAyhZ/WlV2QXnBlqAtUHYvjV6zXGeJhr5QypKTb27fLvzhTNeuZ/v3z1NRoROxCULs5AXGxBoEROn613J4JWGc13j8ujmkbPGnVWUXVhZuC94GZffU6DndZbpo6H8PSmzhHdMbS+zTvXqlTZtG11YlCLHzJCTWU9ln/1JbQynXV3PxkOeElC3+tKrskuqS7SHboex31N+Z4jRFNPS/52lSUtrKlcJD2N26XfzyyxILi4bY1/z/EIKQMBr4ccGafptmWalucXr0ms6/bgtStvjTqrLLasp2hu6EsrurdZ/kOEk09D9G0NBw78QJ5vJPZ/v3v7VgAWtxEgTxktTz+JVRXOmV9qYH3F7j+ddtQcoWf1pVdmVt5e7Q3VD222pvj7cfLxr6H1N+7drV2bPju3aFsi8PG1ZiYsJanARBvAx1XP5Fq5Cn4ZzScI7wfrsd7w7XpGzxp1VlV9dV7wnbA2V3U+02xnaMaOh/DErs8199BV+f+/DDOytW1EdFsRYnQRAvTEk4J1Ddd8oIo0Tz4KroDnoKFilb/GlV2XUNdVLhUlB2F9Uuo2xGiYY+fwSCqnv3UqZMSejRA8pO+e23In191rIkCOKFKQiO9VHxWfinqfE+txzf6PoOcD5fq5CyxZ9Wld0gaNgfsR/KBr9b/y4QCEQNzxdBXV2OsXHil18KS+wBAzLXrKETRQhCjGS4RervddXc7vzIr0Ocf90WpGzxp1VlIwciD3RR6QJl/2b1W0NDg2jo86Shoa64OHny5ITevYUl9siRT06eZC1IgiBejPyg2Pvewos9cU8F5gZEd5Dzr9uClC3+tKXsQ1GHuql2Y5RdU18jGvocaaisLI6PP/vRR/B1wjvvZK5eXePvz1qQBEG8AIWhwuMhLopet90jWU0dE1K2+NOWsmWiZbqrdYeyf7X8taK2QvB/z3tspObRo5s7dzIl9sWvvnp87BhrKRIE8QKUR3H91HzXTbdQ3dzhzr9uC1K2+NOWsuVi5Hqo94Cyf7H85Wn10+dUtqChoezatTP9+sV16QJlZyxaVO7gwFqKBEG8AFF6ATNGGWvvdOmA51+3BSlb/GlL2QqxCr00ekHZP1v8XFRZ9Jz/gazJzX1gZBTXeC52Qs+euUeONER3jnKAIDo4N10izpgGZftE1Xekn6Q/G1K2+NOWshU5iu9pvccoO788/zmV/fTSpeRJk+BrcG3ChFJLy//reKf3E0QnoiScE6nrb7zf7aFfdEUUtyOfH9ISUrb405aylXhKfU/2hbJ/svgptyy3QfDvJ43UFhY+srU93bcvfB3frdtDGZnaoCDWIiQI4vkpCI71VfVdMdncV8XnaUTnuysIKVv8aUvZynzl/tr9oezh5sNzSnOeR9kosdM3bBCW2F26XB469KmVFd16hiBemLIIboSO/7Z51mpbhf9v7LC/l3kGpGzxpy1lq8SpvK/zPpT9o/mPWcVZ/6psQX19nqfnhaFDmRL73qZNdG4fQbwwdVz+k5DYcG1/xtcd/PzrtiBliz9tKVstXu1D3Q8ZZd8pvPOvyq7Jy7t34oTwIlBdupzt27fMxoaus0oQL0ZRGOeKTehtt4iqGG5JGIfV2okgZYs/bSlbI0HjY72PoewfzH64WXCzXlAvamgjRTExV+fORYmd0KNHyu+/0y/UCeLFqIjmBaj7bpptaXnYvROdHNIqpGzxpy1lnzx98lP9T6Hs782+T8tL+1dl5xgYnBs0CMo+N2DAg4MH6bp9BPFihJ70++tXo851/nVbkLLFn7aUrXNGZ5DBICj7O7PvUh+n1jc8S9mVmZk3t26N79YtrkuXi198Ue3jQzfkJYj/ioAfVxnNXTfdwlrGI8urM51/3RakbPGnLWXrndX74tQXQmWbfnfl0ZVnKzvP0zNp3DiU2Gf69EmfMYNOFCGI/0p5FPeGc3h1DC/ZNhT1dS1HEn7QQMr+Ozk5OR4eHoqKitra2nw+XzS0MZWVlRcuXMBwFRUVZWXlc+fOlZWVidpapC1lG5wz+MrwKyj7W9NvLzy4UNdQJ2poGYEg88iRcwMHQtkXv/46T1mZtdgIgng2BcGxAeq+exbZPPLv0BdT/a+QskURCASurq4KCgrq6uonTpw4cOBAQUFBfb2oEL527ZqhoaGsrKyJiQlaDx48mJ6ezjS1TFvKNkw0HGw0GMoeZjLsXPa5Zyi7Kjv7+pIl8e+8E9+tW+qoUdU+PqzFRhDEM8gPEl6fb8d8a5XNTk+CYyXgeEgTpGxRioqKpKWlDQwMSkpKkpOT161bFx0djeKaaQ0NDZWTk/P19a2trU1JSZk8efLZs2eZppZpS9kmF0yGmgyFsr8x+SYhK+EZys7z8rr8228osc/275+xdCkdxSaI/8QVm1C1rU7H1js88pO0C/KQskWBpuXl5T09PfE4OztbT09PRUWluLiYaU1ISDh+/Li2tva5c+fs7e23b9+Ouptpapm2lG120WyY6TBG2bx7vLaULaivv7VjB3N17Cs//EBHRQji+SkJ5xSHcS5YhoRr+z+UOF8DUrYofD4fUg4JCcHj3NxcW1tbKSmpJ0+eMK2owZ2cnEaNGvXbb78NGTLEw8OjsLCQaWIiEAgqKiowsKCgQFlZ2cjISNTQLBaXLL43+x7KHmo8NDYztrahVtTQPAJBzePHTIkd37XrjenTq7y9WcuMIIhWqYrhBWv6hZ30u+8tsWfEkrJFebayz58/b2Bg4Orq+vTp0/j4+GnTpp09e7b5pfgqKysdHR3nzJkzbty4H3/8UUtLS9TQLFaXrX40+xHKHmI8JPJ2ZKvKFtTWPnZxOf/118xRkawtW+i6fQTxnISd9J/8m5HBHtfHgZ3+/Ou2IGWLkpKSIi8vj/IZj+/fv6+jo6OmptZ0YMTBwUFaWvru3bvQNOw8ZcqUiIiI6mZzraGhAZ3xQvQ5cuSIoaGhqKFZbK/Y/mTxE5Q92Hhw2K2w2vpWlN1QXX378GHRUZHvv6cb0BDE81AZzYvSDYCvTQ+43fWMlKRTRFiQskUpKSk5fPiwnp4ezHvlypVVq1Zxudyqqiqm1dnZec+ePais6+vrs7Ozx44dy+Fwms4nYaWtY9n2Sfa/WP4iVLbR4KD0oFaVXV9WljRuXELPnlB22pQppRYWrAVGEERLqmN4l6xCXY96ZXpESsb5121ByhYF5bOXl5dcY2RlZWVkZAoLC2NiYi5evIgHSUlJEPG+ffvU1dXRino8IyND9MoWaUvZjsmOv1n9BmV/bfS13w2/lnfsFdTWVmZkCE/H7to1oXv3exs30nVFCOLZFATHJpgEBWn4FYdxCkNiJdvXgJT9dx4+fBgYGKitrW1ubs6cw5eQkICKG3V3eXk53G1iYnKqMTB4RUUF86qWaUvZzinOv1v/DmV/ZfiV13Wvlsquf/q0ICCAuaFB4ief5Coo0IFsgngG+UHC+xVIr7TzOP6m/JeelC3+tKVst1S3UTajoOwvDb90T3Vvqeyahw9v7dzJHBVJGTmyUFeXtbQIgmiiIoobouW3d7GN4joHiTyfr1VI2eJPW8r2uOoxxnYMo2xU3C2VXZ6Wdvb995k782YsXlzu6MhaWgRBMNTz+NC03h4XpQ2Ob46vASlb/GlL2V7XvMbZjYOyvzj1hX2SfXX9P2a6oL7+6cWLkDVDzr59dCCbIFqlhiP0dVkEB5SGd+L7FbwApGzxpy1l+6T5jLcfD2V/fupzm8s2LGXXFRU9dnVtUnaesrKADmQTRAsE/Lgwbf+F402DNf1Qa3fS+4G9MKRs8actZfvd8JvoMJFRtuUlS5ayqzIzM+XkGF8n/fRTiakpa1ERBFERxQ3X9p8ywsjsoFuOz5t40w9StvjTlrID0gP+cvwLyv7M4DPTC6bVdf+Y6U8vXkyZMoVRdsbixRUuLqxFRRBvOA18frZP1JppFsb73e64S/j5121ByhZ/2lJ28M3gKU5ToOxBBoMMEw2bK1tQW1sQHHz2ww/h6/iuXXMVFOhANkE0pzSck+Mb/SQ41uO49wPf6DfT14CULf60pezQW6HTnadD2QP1B+qf1a+qE/20EqnNz88xMhKW2F26nHv//VJzc7oNDUE0URAcG6LlZyntURzGecP/x0PKFn/aUnZ4RvhMl5lQ9qf6n2qf1m6u7PKrVzP27BGW2I33NKh0c2MtJ4J4YykKjfVX82XOvy4IimW1vmmQssWftpQddTtqtttsRtkaCRrNlV0YHZ08eTKUndC9+90NG2r8/VnLiSDeWOKNgxhfP/B9g86/bgtStvjTlrJj7sTMdZ8LZX+i94kKX+VvZQsEufb2zNX7Enr0KNDSqo+IYC0ngngDqefy63n8EC2/YE0/ybu/zItByhZ/2lI2J5Mz32M+lP2x3sfHeccr60R3Kat98uTeiRPCHz126XK2b98qb2+6cxhBQNYR2v5xRkEP/KKrY3kNEnT/xpeBlC3+tKVs3l3eQs+FUPZHeh8pchSblF2WknJzyxamxE7+7be60FDWQiKIN43yxvOvp/1u7K/mWxT6Zv2+8dmQssWftpQdlxW3xGuJUNm6H8nHyFfWipRdGBmZOmsWlH2mT5+MRYvqIyNZC4kg3iiKwzhhJ/3njTXR3+N6+009/7otSNniT1vKPn3/9DLvZVD2h7ofSkdJNyn7kY3NxR9/hLITP/zwkZxcQ4zE3gOJIJ6HotDYaL2AE5scs7yiyNcsSNniT1vKPpt9doXPCij7A50PDkQcqKgVXnG7obLy7tGjwmtkd+ly8csvyx0c6Ixs4o3lSXBsin3YZZvQvMCYR37RdPy6JaRs8actZSfmJK7yXQVlv6/zvlSYFKPsyszMG+vXM6f3Jf/6q7DEpqtBEW8kxWGcQA0/1S1Onspvyv0KXgBStvjTlrIvPLiwxm8No+xdobsYZRfFxKTOmAFln+3f/9bChazFQxBvCDUcfthJ/02zLI+spfOvnwUpW/xpS9mXH15e778eyh6gPWBb8DZG2Y+srC79+iuUfeHzzx8cPMhaPATxhnDLNWLVFIsTmxzp/OtnQ8oWf9pSdtKjpA0BG6Ds/tr9NwVtKq8tx8BMWdlzn3wCZV/54YdivKrFEiIIyUbAjysOja2M5mZ6RD7yj66n49fPhJQt/rSl7JTclM2Bm6Hsftr91gesh7Lry8vTli+P7949rkuX1NGjawIDWYuHICSbskhuhI7/uukW6S7htRwe+fpfIWWLP20p++rjq1uDtgqVfbLfar/V5TXl5TduMNfIPv3uuzfnzqVzRYg3iqJQTrCm36IJZgZ7XHP9Y960+8u8GKRs8actZaflp20P2Q5l9z3Zd7nPcii7ICjo8u+/Q9nnBw7M2rqVtWwIQoJp4PFT7MMOLrdV2tB4/jWX6uvngpQt/rSl7PSC9J2hO6HsPif7LPFaUlZTlqOvf+Gbb4QHsr/77rGiImvZEISkUhrBKQqNve4YZnLA7SGdf/1fIGWLP20p+9aTW7vDdguVrdVnocdCKDtjz56zH38MZaeOHl1saMhaNgQhkZSEc2L1A7inAiujuHRP6v8KKVv8aUvZtwtv7w3fC2W/p/XeXLc5pRXFyX/9Ff/OO3FdutyYPr3S05O1bAhC8qiK4YWe9F8/w/Loese8QLo2w3+GlC3+tKXszKLM/RH7oex3td6d6Ty94HbapeHDUWLHv/32nZUr6YKrxJtAlG7AzNHGJzY5PqDzr18IUrb405ay7xXfOxh5EMrurdl7it2k+2EBF4YOhbITP/00R0qKtWAIQiJxPertq+KT6RFZT/9vfCFI2eJPW8rOLsk+HHWYUfZEq7E3dDXODRok/N/j998/PnaMtWAIQpIoi+TyjQKTbELveEQWhsTS9fleGFK2+NOWsh+UPpCJloGye2n0Gmc64vKWdWfefx/KvjZ+PP3vkZBgikI5IVp+q6danDYJqoiiHx+8FKRs8actZT96+kguRg7K7qnRc5Th8ISxIxJ694ayby1cWOHqylowBCEZFAp97b9jvrXiWof73lF1dDzk5SBliz9tKftx2WOFWAUou4d6j990hkZ/0j++W7e4Ll3ubdpUFxbGWjAEIRlA0/bynjIr7R/4RjfQ7xtfGlK2+NOWsgsqCo5yjkLZ76h1H642KKJXV5TYZ95998GhQ3S6CCF5lEZwcnyjsryisn2i8IDVSrwYpGzxpy1lF1YWKvGUoOzuqm9/e3RARM8uUPbl774r0NBgLRWC6OxUx/LCdfxVNjsGqPuymoiXgZQt/rSl7OKqYmW+MpT9tkrXIQrvhTcqO23y5FIzM9ZSIYjOTrRewOw/TNS3OtH9CsSLJCs7Ozs7KirKysrKwMAAf6OjowsKCurr60XN7Za2lF1aXaoSpyJU9okuX8r0CO/5FpR9e+nScicn1lIhiM5LPY9/2z1iwZ+mwvuju0XQ/xvFiwQqG1KOi4s7fvy4jIzMiRMn9PT0DA0NdXV1lZWVpaSkNDU1L168WFEhvCNMO6UtZZfVlKnFq0HZ3U68NehwV0bZ2bt21fj7s5YKQXRS4OuScE5pOCdKNyDTI7KGQ/+kETMSqOza2lonJydnZ+egoKCEhITk5OTr16/jb3x8vK+vL4ZzOJzc3FxR73ZIW8qurK3UTNCEsrueeOuTw29B2fHduz9WUhLeorfFgiGITkdRKCdGP0B1i1N5JLcqmu5X0C5IZpV94cKF4uLihoYG0aBmyc/Pv3v3bmFhoeh5O6QtZVfXV588fZJR9ofSQmWf+/DDQj091iIhiM5IYUhsiJbfnkU2ShscK4SX6GN3IMSCBCobpkZZndJG2vWQCJO2lF3bUKtzRgfK7nLirQGyb4X1fOvKjz+W0P8eic5PA59/xixYZpXd4RV2Ob5R5Ov2QwKVXVNTo6qqqtBG7t+/L+rXbmlL2Q2CBr2zelA26CsvVPb1SZPK7OxYi4QgOheV0TyU1UEafto7XXJ86Pzr9kUClV1XV+fv7+/RRgoKCkT92i1tKVsgEBicM2CU/d6Rt0J7vpWxZAn9VJ3o1NRy+OfMQ265RtTE8uh+Ba8ACVQ2Ultbi1q71cCbok7tlraUjRgmGnZV7cooO6j3W3c3b6ry9mYtEoLoRETrBcwZY2Ik5VoYEstqItoDCVR2dXX10aNHD7aRe/fuifq1W56hbKPzRm+rvc0oO+Ddt7L3768JCGAtEoLoFFTH8qL1hb+X0d/jkkHnX78qJFDZKLFdXV1t20heXp6oX7vlGco2uWDSQ70HlP3ukbf8+rz1+IRyXXg4a5EQRKegIoprJe1hedg93UV4VITVSrQTEqhsgUDw8OHDnDZSU1Mj6tdueYayzS6Y9dboJVS24ls+A7oWmhgJYml3kuhkFIVyUuzDbrpEXLYOfeQfTb+XeZVI5rFsJhUVFampqYGBgW5ubs7/S35+vqi53dKWsqtr6/XPGPdSEyq7t2IX28/7xRnaPIjg1fDYS4UgOiyFobGhJ/2VNzletAqhi6m+eiRZ2YmJifr6+suXLx87diz+Tpw4ccKECdevXxc1t1taVfbTqtpzmU+2eav1UBYeGOl1tKvGj0P3HLE1tYtJDuKW032ViM5AaTgnTNt//1Lbg8tsc3zo/OvXgCQr+9ixY1paWhDoli1bbt68aWpqKisr+1r+/VhXL+DfzNtgf2HQ0b09jgmV3fNYN9lRv43aaztMMWy/UdTFQG4dffuJDs81h3DFdQ6Hlttl0/nXrwlJVra8vLyzs7OXl9fu3bvr6uoqKytnzJiRmpoqam63tFR2dmHFRvsLQxXDP1XY1VuxO5TdQ+ntnZMm/iLl8JVC2DeKYXKmUbmRdECQ6LjU8/g1sbz8oJg77pF0PdXXiCQr++jRoy4uLuHh4du3b4e7g4ODR40alZycLGput7RUtk3CnSl6/MFHwgfK7eqjIFR29+PvzFy58btD7lD24CNhC7QifNzo/5BExyVGP+DEJsfTpsH1XD5d7+k1IsnKPnfu3LVr17Kysnx8fPbu3Yta28HB4bX8+nGTw4Xvj0fCzp/L7uwv1w3Kflu5x+it+7+R8cJA8NOxMAUz2tMkOiJVMbxY/YC5Y00sDrvnUH39upFkZZeVlZWXl9fU1Dx+/DgxMfHMmTN4UFdXJ2put7RU9lzjhK+PhH0tH/rV4W0fSneBsrue6PnNQeWv5XwZZaPQ3nmKlE10OKpjeVftw+aPMz250/mGczidf/3akWRl83i8K1eu1NbW4rFAIIDBPTw8XstPaeabnB58JHyYtN+3+9Z/ckh4jZGuJ3oNkj/5pYI/o+whR8J2kbKJDkYdl18azrnnFXVyp8sdd7pfQYdAkpV9/PhxJycn5rczDQ0NxcXFixYtei0n+e11v/LziegRe51/37Zk4MFGZav0+uSISNlfK4SNPB6uakm7nEQHoiiUc8Y0OFjTr47DfxrBaaBLPnUMJFnZcnJydnZ2TVV2aWnptGnTXssZI14X7880jJ+4zXzS2hmf7Rcqu4tKr48VNRhlo8ReoR0R4cVhLRuCeF2UhHPCtf0PLrM9stahlg6GdCQkWdlqamq6urooq/Py8nJycuLj4xcvXpyWliZqbre0VHZeadVBz+QFO01nrPzri31Nylb7UsEPJfZvx8M1raKfRFMVQ3QIBPw4vmHgjvnWB5bZ0vnXHQ1JVnZycvLOnTvHjRu3YcOGpUuXDhs2zNfXt7i4WNTcbmmpbOTqgxITNdtli0Z/JcUou+fHiieGHPH7/mjYcYvojDAqZIgOAXyNslpzm/PJHc45vuTrDockK7umpoYprr29vQMDA1Ful5WVtXpDSPGmVWXX1jfcdHBVXT968F7mWHbPz5WUthoEBLjH5kTy6KePRAfhklVIXmBMYUhsSTinnq6n2vGQZGULBIL09HQ3NzcvL6+CggLo++bNm5WVlaLmdkurykYeWFrqLvnhm91CZfdS66XjpJISElYSy68nXxMdgKoYHscgcPEEsyjdgLJILquV6CBIsrLj4uIMDQ337Nlz/Pjxurq6hw8fHjly5O7du6Lmdktbyr6vpaUz5dNvdwmV/a56Tx+fE6Uculg20SEoDedE6QUs+NNUc5vzTZeIWrpOWUdFkpWtra2toaGhpKS0e/duVNwlJSXTp09/LWeMIIL6+jtyclojen23Q6RsD2+lElI20TF46B9tcchdbpV9hlsEnX/dkZFkZSsoKDg5OXl7e0tJSTEn+U2dOvV1KbuuuPjm9u2a33f5YbtI2S5eR4s5YazlQRCvmKIwzn3vqFuuEdxTgTk+UXQ/sA6OJCtbV1fXwsLCwMBg06ZNGRkZUVFRq1atSk9PFzW3W1pVdkV6etqKFSe/fWv4NpGyHTwVSNnE66U0ghuh428j65FkG8pqIjomkqzsCxcuqKmpzZo1a8qUKdLS0pMmTfLw8CgqKhI1t1taVXYxj5c6c6b2sLd+3iJUdm/1nrYeckUcWk+I14bwfrt6AVvmWB1cZpvtTefzdQ4kWdkCgeDx48fR0dEmJib29vbXr18vKSmpr68XNbdbWlV2nofHlTFjdL5569fNjcpW62nlIV0YS8omXhsJJkGLJ5gprLGHrwUtWomOicQqu7CwMCMj4+HDh0+fPs3Pz09OTtbQ0Jg8efK1a9dEPVrk5s2bqMqnT5++Zs0a1OOiof/LgwcPbG1tZ8+ePWfOHCkpqVu3bokaWqRVZT8wMbk4fLju0LdGbGKU3cPM/eCT2BDW8iCIV4OAz8/xjTprFnzLNYKOX3ciJFDZdXV1YWFhW7ZsWbt27fbt27W1ta2srJYtWyYrK+vl5QWVi/r9Mw0NDXp6ekpKSiEhIU5OTnhtVlZW05Va4f3AwMC9e/fGxMScPXs2ISHhGQdYWlX2PWXlxK++0h/61siNImUbue1/EhvMWh4E0d5URvO4hoGGUq6ZHpEVUdwaOp+vUyGByq6qqtq8eTNM7enp6ebmJi8vv3z5cn19/cTExJKSkrZ+/YgiWkZGxsLCorKyMi0tDXU0Cu2ysjKmFUW6gYGBkZERc11AfMQzfkXZqrIz9u498/77BkPeGr1eqOxeaj0MXPcWkLKJVwtz/vXSSWa2sh55gTGsVqLjI4HKrqiomDhx4uXLl+vr6+FfS0tLlNiPHz8WNbeRixcvKioq+vr64jH0DTvjaVMpjdL78OHDhoaGrq6uzs7OcDrrV5Sox1NTU7GRcHR0RIWOLYSooTGChoa0VasSevQwHPzWmHUiZeu47MqPDWItD4JoP6pieIkWITsXWEuvsLtN5193TiRT2WPHjoV8rzQG5faKFSuYx0h5ebmo3z/D5/OPHz8ONeNxbm6ura0tCu0nT54wrV5eXqtWrdq/f//Ro0cPHTqkoqJy69at5oU2qm8ej6eurn7kyJHZs2fr6uqKGhpTV1SUOmtW3FtvGQ3p+uf6royyNZ2355GyiVdFHZdfHMaJMw6SWSm8Pzrdv7GTIpnKnjp16p49e2BPZN26dTNnzmQeI/fv3xf1+2dOnz597NixwMBAPH706JGFhQXK6qYD397e3lC2pqYm1Pzw4cPJkyeHhoa2dbmSlgdGylJSkidMgLJNh3WftKVno7LfUXPakhcTyFoeBNEelEVyrzmGpzuH13J4peF0ZfZOjGT++xHmRV3capoKZ1bS09NlZWWdnZ1RO9+9exd1tIGBQUlJCdMaFRWFGjw4OBiP0WHt2rWenp5t3fm3pbILQkIujxgBZZt/32vqzr5Qdk+1d5QdNz4mZRPtTw2HF6MfsHWO1YlNjnS/gs6OBCobgbXbikAgEHX6Z2praw8dOqSmppaXl5eYmIjCPDU1lflnI5KRkWFubq6lpQVfQ/rjx49HlV1VVcW0stJS2Q+trC589x2Ubfl97+k7+jPKPuqwPjcmgLU8CELscE8FLhpvqrjWIcuLfi/T6ZFAZcOzOjo658+fbzrfoymQLI/HCwkJuX37tmjQ/wKV8/l8WBuyXrZsGexcXV3t5OQENT969AiPExISdu7cuXDhwgULFqAVlXhb9m+p7HsqKolffgll24/8eIHCEEbZCvZrH5GyiXamlsOTW21vtM+t8agIHb/u9EjmgZHw8HBY9eDBg4qKiqdOnUKBrK+vLy8vLyUlpaure/bs2aYjHs1TWlp67do1Lpd75syZe/fuYUhmZuaDBw+YY9bFxcXJycnQelxcHFoh8cYXtZKWys7Yu/fsxx9D2U7jvliq8hOjbBm7lY9i/FnLgyDERVUML8U+rCyCk2wXmu0TVUOHRCQCCVQ2it8nT57Ex8e7ublZWFgYGhrC1/gLcWPIuXPn0NpWgSyWtFT29RUrTvfrB2W7Thq66uRoKLuH2jsH7ZY/jPZjLQ+CEAul4ZxY/YCtc61uukSQrCUJyTyWzQTldl5eXlpaGqrjGzdu5OfnP+P3L2IMW9kCQfKkSfHvvANle0z7Yb3hX43K7r7PdskDUjbRDhSHcaJ0A3YvtDm4zPauZyQdD5EkJFnZryv/ULZAUF9efunnn+O6dInv2tVn7ogtFrMYZe+2WZgT7ctaHgTx8qCy1t3tsmeRTbY3nX8taZCyxZ/myhbU11dmZjKni5x+993A5RO2W89nlL3deh4pmxAv5VHckjDOVYewCG3/+3Q9VUmElC3+/EPZtbXF8fHnBw+Gss8PHBi+Ze4u6wVQ9jtq3Tdbzs6O9mEtD4J4Yeq4fM6pQF8Vn+uOdIM6iUXylV1XV/eMszvaI82V3VBdne/rm/jFF1D2paFDY/at3GuzSKhs1e4bLGfej/JmLQ+CeGG4hsL7o6tsdnrgG81qIiQGSVZ2VlaWgYHBokWLlJSUampq0tPTDQ0NHzx4IGput/xD2ZWVDywszg0cCGUn//wzX37LPtsljcp+e43FtCxSNiEO6nn8M6ZBSyaaaW53TnOi868lGUlWtpWVlba2trS09K5duxoaGvLz8xcvXnz9+nVRc7ulubLrKyqyVFWZk7JTR48+c2LvQbtljLJXmk/JivJiLQ+CeAGg7BtO4a5Hva85hFfTKX0SjSQrW1FR0dbW1sXFRUpKCk+rq6snT56ckpLCtLZf/qHssrLbhw6def99KPv6pEnntWWl7VZA2d1V315qNukeKZt4OUrCOYkWIV7KPoUhsYB8LfFIsrJVVVWhbOYyqpWVlRcvXpw/f/6rrrKfPk3ftOlM4+9o0mfOvGx0XNZ+FaPshabj70Z5spYHQTw/xWGcaL2AwyvsjPe7VUaTrN8IJFnZ3t7eZmZm0tLSS5cuDQgIOH78uJqaWk5Ojqi53dJc2XXFxcmTJyf07g1l31q4MNVKS8FhDaPseSZjM0nZxItSE8uLNw46uNx290Kb+95RDXT+9ZuBJCu7sLDQy8trw4YN48ePnz59uqKi4p07d2pra0XN7Zbmyq4tKDjTv39cly5Q9p2VK9OcDRUd1jHKnmXyR2akB2t5EMTzUM/jFwTFWEm7Sy22uU/X53uTkGRlp6WlZWdni540Xuf67NmzpaWlouftFpayT7/3HnwN7m3ceMvLSslxA5T9tmq3acYj75Cyif8OCuq8wJjScE5VDO9pBJfVSkg2kqzsEydOeHp6ip40XpR1yZIlr/hYdk1+/pkBA5gqO3vXrrvBLiccNzHKnmz02+1Id9byIIh/hWcYuGaaheMRz1oOT8BntxKSjQQqG9X0jRs3AgMDFyxYsGfPHlgbcXd3t7a2HjduXGpqqqhfu6VJ2YK6usrbt5kqO/7ttx8cPHg/0lvVaQuj7Amnfs6IdGMtD4J4BnVcPny9ZKKZ1nbnmy4R5Os3EAlUdn19fWJiIqQ5ceLEZcuWaf4vOjo6tra2eXl5on7tliZlN1RVlZw7x/zv8Wz//rlHjuRE+2o4b2OUPc7gx4wIUjbxvEDQTyM4B5baau90SbUPo/P53kwkUNkCgSAnJ+f8+fPm5uaurq7xjUlISLhw4UJ+fn5dXZ2oX7ulSdn1T58WhIQk9OoFZV/88st8FZWH0X4nnXcwyh5t8P2tCFfW8iCIVimL5GZ5RT0N53gq+9xxjyRfv7FI8rHswsLCGzduQNZRUVER/0txcbGoud3SpOy64uJcF5eEnj2h7Cvff/9ESys3JkDHZReU3U216wj9YTdJ2cRzUBLG4RoGGu5zyw+KpZP53nAkWdnXr1+3s7PbsWPH3LlzN27ciL8LFy68efOmqLnd0qTs2oKCB8bGCT16QNnJv/5aqKubFxOo77qHUfYvekNvRriwlgdBsHgawY3RDzi43HbnAusHvtENdPz6zUaSla2qqnr06FFZWdmVK1fGxsbKy8sfOHCAualju6ZJ2dUPH2bKyzP3o7n6xx/Fhob5sUGGblKMsn/U/TqdlE38G5esQw+vsNs535ruj04ASVY2HO3s7Ozl5bVnzx48ra6unjVr1tWrV5nW9svfys7OvrV7d3z37lD2tT//LDE2fhIbYuK2n1H2tzpf3IhwZi0PgmBx2jQ4WNOP7ldAMEiyshUVFV1dXYODg3fu3BkfH5+cnDx69OikpCRRc7ulubIz9uxhlH117NhiI6Oi2FAz94NQdlfVrkN0BqWFk7KJNuEbBUbq+N9yjaiI5tZx6RA2IUSSlR0ZGZmYmJienm5lZbVhw4Z169bp6Og8evRI1NxuaVJ21b17NzZsYJR9c968MlvbEk6YpfthRtlfaX+SFu7EWh4EAWo5/DijoGWTzBwVPB8HxLBaiTcZSVb2kydPioqKKioq7t69Gx4ejnKbw+GUlJSImtstTcquvH376oIF8W+/DWVnrl5d6epayomw9pARKluly2daH14Pd2QtD4KoiOKeMQ1e/pe56manFPuw6hg6n4/4G8lUtkAguHPnzqVLl27duvX06dPy8vKUlBR7e/sTJ040v+pIO6VJ2RXp6Ul//hnfrRuUnbV1a5WPTxk30s5TjlH2J5oDrpGyiRaURQqVfWStw3VHul8BwUYClV1fX3/+/HlDQ0N1dXUTExN/f/+goCAlJaUdO3YoKiq+SmWXX7t24dtv47p2hbJz9u+vDQqq4EY5eipA2V1Uunyo0e9quANreRBvMiXhnDSn8ItWIUWhsQ/9oul+YERLJFDZVVVVmzdv3r17t66uroqKyoYNGxYsWLBr164zZ87U1NSIOrVnmpRdlpR05oMPmGtC5R45Uh8ZWcWLdvZSZJTdX+O91DBSNiECxXWsQeCJzY5G++kyBkSbSKCyKyoqJkyYcPXqVYFAUFdXZ21tvWTJkldwzdWmiJQtEJReuMAcFYG189XUBFxuDS/WzfsYo+w+6r1Tw+xZy4N4M6nn8eOMg3bOt94hPP86ktVKEE1IprKnTJnSdJFVb29v5t6PryyMsgW1tcV8vtDXb72V0KNHoY4OZnctn+Ppc7xR2W/1UuuRQsomGnngFy212ObAMtu7npF0fT7iGUigssvLy4cMGTJ58uQFjRk5cuTgwYPxYP78+fibnp4u6tduYZRdX15eEBTEKPvsgAFFBgaY3fV8ro/PCSgbdFd9OznMrvnCIN5MngTHlkdx73lF3feOovOviWcjgcquq6vz9fX18/MLaAyqbA8PDzzw9/fH38LCQlG/dguj7NonTx7Z2THKPj9wYLGREWa3gM/381VllN1FpUtSqG3zhUG8aTD3b9y5wPqCZQge19Mln4h/QwKVLRAIysrKnraR+vp6Ub92C6PsmkePck6dElXZ/foV6eszczzAVx2yZqx9OdRaEEdr6RtKWSSXbxS0coq5ymbHDLeIBj59E4h/RwKV/drDKLsqK+ueisrfVbahITPHg/w0uql0Y5R9MdSqgZT9RiLgx93zilJYa39wmS2df008P6Rs8YdRduWdO5kKCoyyrwwfXmphwczxYD/Nd1S7M8pODDGv59O6+saB+rowJDbTI1Jvt8s9z0g6/5p4fiRc2cx5fq/mdOym/K1seXmRsn/4odTcnJnjIX5avdR6MMo+E2Rax6c7ZL9ZlEdxT5sGh530Lw3nCOhgCPEfkWRlQ9bFxcUZGRm3bt1i3P0qj2VX3LyZsX8/o+z0mTPLnURXgAr1P/meWi9G2fFBxrU8TtPCICQe4fWejIO2zROef53jE81qJYh/RZKVnZSUdODAgV9++WXnzp3V1dV4unXrVuhb1NxuYZRdlpSUtnIlo+xbixZVuIruGRbmr91P/V1G2dzAU9W82KaFQUg88SZBq6aYH1pud9s9gu4vQ7wAkqxsTU1NQ0NDXV3dPXv2oMouKSmZPn16amqqqLndwii75PTppPHjGWVnrllT5enJzPGIAJ0PNPsyyo4N0K/i0aU13yAidf29T/hccwij49fEiyHJyj5y5IiTk5Onpyfz60cU2pMnT05JSWFa2y8iZZ89m/zXXyJlr13bpOzIAN2PNfszysbjSh7tHUs+1bG8BJOgc+bBt1wj8oNi6fwQ4oWRZGWjvjY1NdXT09u+fXtWVpa3t/eqVate2e16SxISkv78k1F29u7d1X5+zByPCtAbqPUBo+ww/5MVXLpBlIRTFsmNMw5aPdUiVMuvOIz+dUG8FJKs7Pj4eBsbm717986fP19fX19WVtbZ2bmgoEDU3G5hlF3M5V765RdG2blHjtSGhDBzPDpA//OTHzHKDvbTLOfSNYAkGfg63jho90KbvYts6Pxr4uWRZGWXlJScPn3a0NBQRkbm2LFjHh4epaWlAoFA1NxuYZRdFB19YehQRtl5Kip14eHMHI8N1P9K+1NG2QG+ak+5EU0Lg5A8HgfGeCr7bJ1jddczkq4fQrw8kqzsK1euJCYm3rt378mTJ9WvcAoZZReGh58bNIhR9pOTJ+ujRAdAOIGnhmoPYpTt66tSyhGpnJAwKqK5uQExOb7RjwNisryjBC06EMQLIMnK1tLSmjhx4pIlS0xNTa9evVpeXg5xv7Iq+0lw8Jn+/RllFxsZNcSKTubjBp4apvM5o2xP7+MlnLCmhUFIDPVcfoJJkNJGR/ND7qwmgngZJFnZFRUVOTk5MTExR48enT59+l9//aWpqfnK7rCe7+fH3FsdPLWyEvBEBzF5QYY/6H7FKNvN61gxJ7RpYRASA3y9eqqF3Cr7TA/6XwUhTiRZ2Uh9fX1ZWdmtW7dcXV0XL178zTffXLlyRdTWboGyjQwMHru5MbcQS+jevdzh7xuG8YOMftIdzCjb2UuxMJaULWk88o/eNtfqxGbHVPuwGg79v5EQJ5Ks7MLCwrNnz5qamkpLS8vIyKirq9vY2OTl5Yma2y1Q9iktrUe2tsISu0uXM/36VTg7N83x+CDjX/WHMsp28FR4Eis6k4SQABp4/CfBsU8jOLH6AdedwqtiyNeEmJFkZfv6+mpoaOjp6VlYWPj7+6ekpJSWljY0NIia2y1QtoGq6gNTUyg7vmvXxE8/bfq1OkgINvnd4FtG2bYesgWxwU1NRKemLJJ7xixYca1DXmBMdQzdr4BoFyRZ2cHBwe7u7sw/HkWDXkmgbD0lpWw9PaGyu3W7+NVXlW5/3zP7dLDpaIPvGWVbekjnxwQ1NRGdl9JwDs8wUGqJza4F1o/8o+n6IUQ7IYHKRh197dq1ioqK9PT0My3y9OlTUb92C5StKy+fpa7OKPvKd99Veng0zfGzwWbjTg1nlG3mfvBxTGBTE9FJEfCxIoUd2+CwebbVPS86/5poRyRQ2dXV1QoKCllZWRoaGrObZdasWfh748YNUb92C5Stc+jQ3WPHhMp+++2UESOaLjACzgWbTzD8mVG2sdv+3JiApiaiM1Idy6uK5vINA1U2O93zpPNDiPZFApUtEAgeP35cW1ubl5d37969+81y9+7dV/CbGij75N69dw4fFiq7e/er48ZVeXs3zfHEEIu/DH9jlH3KTepRjH9TE9HpaODHXbENTbUPq4iiW1UQrwLJPJbd0NAAcfv6+p4+fbr+f4GsDQwMcnJyRJ3aLTC21rZtGbt3Q9kJ77xzY8aMal/fpjl+PsRyqtHvjLL1XPc8jBZdLorojCSYBK2dZqG+1elxAF1El3gVSPK/H1VUVDw9PUVP/u//ampqFi9efP36ddHzdguUrbF+ffrGjYyyMxYtqvH/u5S+EGI1w3gUo2xtl50Pov+2OdGJEPBFv5c5sckx2TaUzr8mXg0SqGwU1PHx8fDm2LFjFyxYoNiYI0eOHDp0aOHChenp6aJ+7RZ8tNqKFdeXLxcqu0ePzLVrawL//h/jpVCr2cZ/MMrWdN6eHe3T1ER0IqBsXxUf0wNul61D6fxr4pUhmcq+cOGCmZkZBL1lyxad/0VXV9fX17ewsFDUr90CZasuXnx13jwo+3TPnlnbt9cG/33y9eVQ63kmYxllqzltuR/192FuolNQFsm96hB2zTE8xS4sxzeafE28SiTz349Pnz598uRJYGAgj8d72JhHjx5B1lVVVa/mpzQn5s1LnT5dqOxevXKkpJoulg2uhNosNB3PKPuE46YsUnanojSCwzcKOrreIVY/gGRNvHoks8q+dOlSWVlZSkpKdIsUFxeL+rVboGzlWbOYu4hB2Q8PH64L/ftCIklhtkvNJjLKPua44V6UV1MT0cGpiOLGmwRJr7TbMNPyHl3/mngdSKCyq6urDxw4cPfu3ePHj0/4ZyZOnPhq/v14fNq0pHHjhMru3fuRgkJd2N9XWE0Os1th/hej7CMOazMj/z5lm+jgZLhFqG11gq/v0vnXxGtCApWNMBfFxt+WYTq0a6BspUmTLv/+O6PsPGXlplvSgNQw+9UWUxlly9mvvhP59w8jiY6MgM8vi+Rke0dledHtOonXhmQqm0lmZmZ2dnZNTc3jx4+9vLxg0jt37uCpqLndgg86NmbMxe+/h7LPvPvuEx2dplvSCOd4uMM6y+mMsqXtVtyOpEvgdwJOmwapbnEK1/av4/LpeAjxGpFkZRsYGAQEBBQXF4eFha1evVpXV3ffvn13794VNbdboGzF335L/OILobLfe6/EzKwh5u/fWVwLd9xoOYtR9kG7ZRmRf18xiuiANPD4p02D106z0NnpkuFGN+okXjOSrGxZWVlHR8f09HRTU1M1NbXU1NQZM2ZcvXpV1NxugbLlhwxJ6N1bqOw+fcodHZvuIgauhzttsZrDKFvKZvHNiL+vy0p0NOq5/Id+0ZtmWSpvdLxE518THQBJVvbRo0fV1dUtLCxkZGQuXrxYXV09depUiFvU3G6BsmU/+wy+Fiq7b98qLy8B9+8LUNwId95hPZ9R9m7rhekRLk1NRIeijssvCePkB8Xq73FNtQ8jXxMdAUlWtoeHB4preXl5LS2tkpKSoqKi7du337p1S9TcIsXFxZcuXQoODo6JicnIyBANbRZI//bt2+7u7tnZ2bW1taKhLdJc2Wf79q0JChLw/z76CUfD1Iyy4W4YvKmJ6Dg8jeBetg71UfGp4/Aro3kNdL8ComMgycouKys7c+ZMQEDA5cuXnzx5cv/+fbi4rfOyBQJBbGystLT0kiVL1q1bp6GhUVFR0fwMEzzOycnB8H79+gUGBpaWlooaWuQfyh4woC7iHwdAb0a4StksZpS91WpuWrhT81aiI1AexY03DpJZZbdljlVlNFdA9ysgOgySrOyampoHDx6cPn3a19cXtfONGzfKy8vb+vUjOsPXMHJ+fv758+fnzp175cqV5qeXVFZWxsXFzZo1a+zYsREREc9StpaWzKBB8HV8167nBw5sfroIuBXpdsBuGaPsTZazr4U7Nm8lOgIXrUL2L7VdP8PyLt0fnehgSLKyExMTt23bNnHixM2bN69cuXL48OHOzs5tXWMkLS1NVlbW1dUV1fS9e/fgbm1t7ZKSElHz//0fSnU1NTU3N7cpU6ZERkY+S9kaGjKffipU9ttvXxw8mKXs25Fuh+1WMMpeZznjavjfN18nOgL1XL69vKfqFqfbbhFUXxMdDUlWtqqqqq6u7vXr1x8+fHj//v2zZ88uWrQIahY1/zMJCQlKSkpBQUF4/OjRIysrq4MHDz558oRpRent6emJDrD51KlTWyq7qqrKz89vy5Yty5cvHzNqlFT//lB2wjvvXPn+e5ay70R6yNmvZpS9xmJaaph981bi9ZJkG3rfOyovMKYgOLaWQ8eviQ6HJCtbWlra1ta2rq4Oj1E7Q7KTJ09OSUlhWlnh8/nHjx8PCQnB49zcXLxQSkqKUTZeGxYWhrr7/PnzcHeryq6vr8/Ozr5w4cKZM2ekdu489MEHQmX36JH8668N0dHN53hmpMcRh7WMsleYT04Js2veSrwu6nn8s2bBm2ZZ+ar4lIRxWK0E0UGQZGWfOnVKQUHBy8srMTGRx+MZGxtv2LCh1VNBkEuXLh05csTHxwePc3JyDA0Njx07VlRUhKc1NTV4Ct1jG4B8/fXXa9asgeLLysoaX8qO5okT0h9+CGWf7tXr6tixLGXfjfI85rCeUfZS80lJYbbNW4nXQmU0D75eP8Py2HoHuv410ZGRZGVfuXLFxMREWVlZV1dX+ItERcXAwEDGwi3z6NEjGRkZMzOzioqK69ev7969G/ouLy9HU21tLYfDUVVV1dTUPHr06JAhQzZv3pyQkMC0toymsrL0++8zyk4dPZql7KwoL2XHTYyyF5mOvxJq07yVeC0Uh3GcFb2kV9gl2ZKviQ6NBCpbIBBApklJSVwu183NzdTU1Nra2sPD4+LFi884mbqhocHIyAhGRlVuZWW1cePG7OzsGzdu3L9/HxIXdWo8qN3qgZHm0VRSOjxggFDZPXum/P47S9n3o7xVnbYwyp5vMu5yqHXzVuIV8zSSm+0Tdcc9Mt446J5nZE0s+Zro0Eigsuvr6+Pi4g4dOrRs2bI1a9bo6+s/583Dbt26paOjs3jx4q1bt3p7e2MIjB8REZGbm8t0QFCkowA/e/ZsWyU2onH06OF+/aDsM++9d2PatOYXGAHZUT4aztsYZc8xGXOJlP36qIjmJpgEmx9y5576+05vBNGRkUBlV1ZWzp07F6Y+f/58WFiYgoKCqqqqqO2VRENR8XDfvlD22b59MxYtYik7J9r3pPMORtkzjUddDLFq3kq8Muq4/HPmwXsW2aydbpFJ518TnQQJVHZFRcXkyZOTk5NRbtfW1np4eOzfv1/U9krSpOwzffrcnDOHpewH0X66LrsYZU8z+v18iGXzVuKVkWQbun2e9e6FNjddIuj36ERnQTKVPW7cuNDQ0PT09Bs3bhgaGq5bt455jKAGF/Vrt6gfOXKoTx8o+9wHH2Rt29b8Mn7gUYz/Kde9jLInG/2WGGLRvJV4NcDR+YExiRYhKXZhNRw6fk10GiTzwMjs2bO3bt26rzGrVq2aNWsW8xjl9iu4Xra6gsKh996DshM//vjBwYMCzj9O8s2NCTB228coe4Lhz+dCzJu3Eu0NczzEaJ/bdcfwimhuNZ0fQnQqJFDZdXV1wcHBvm2k6QeN7Rd1eflD774rVPYnnzySl2cp+3FMoJn7QUbZ404NPxts1ryVaFcqonnw9cZZlnp7XO560v3AiM6HBCr7tUddTu5Q4/0Nzg8c+FhJqfnFskFeTJClx2FG2X8Y/HA62LR5K9F+oL5Od4nYu9hm90IbOv+a6KSQssUfdVnZQ716CZU9aFCBujpL2QWxwTYesoyyRxp8mxBs0ryVaCfg66cRnCSb0O3zrG66RND510QnhZQt/qhJSx/s2RPKvvDZZ4U6OixlP4kNsfeUZ5T9m/438UHGzVuJ9qAympfhFnHVIQziLo+k618TnRgJV3ZDQ0NtbW3TlaHwFH+ZpvaL2sGDB7t3Fyr7iy+KDA0FvH8UdIWxIU5eioyyf9IdHBdk1LyVEDv1PH6iRcjuhdY75llVRP1j80kQnQ5JVjbsnJqa6uDg4O/vD2vn5OSkpKQ841eL4orK3r37G29Jc/HLL0vMzVnKLuKEunkdY5T9o+5X/CDD5q2E2DlnHrxptpXUYpsbTuFUXxOdHUlWdmBgoLKy8oYNG2RkZOrr63Nzc7du3fqMez+KK6r79x/o1k2o7K++emplxVJ2CSfM0/s4o+xvdb7gBpKy25EGHv/UXlfd3S5XbELp+DUhAUiystXU1AwMDHR1dffs2YOKu6ysbNq0aa/gDusqUlL7u3SBsi8NHlzh4sJSdikn3MfnBKPsodqfcQJPNW8lxEUdl59sG/okODbFPuyORySdH0JIBpKsbAUFBScnJ29vbykpqYaGhqKioqlTp75SZQ8dWu3r+3/Nbq8OnnIj/H3VGGV/rf1pbKBB81ZCLFREcxMtQnbOt75oFVIZTbImJAdJVraFhQVzp4KVK1cmJCRYW1vv2rXrzp07ouZ2C5R9oGtXKPvyN9/UBgezlF3GjQzy02CU/cXJj2MC9Ju3Ei9PWST3jFnwviW22+dZJ9uF0e8bCUlCkpV97do1c3Pz9evXL168+OjRo2vWrHn2da7FFeG/Hxur7Mvfflsfyb5EXAU3Ksz/JKPsQVofRAXosToQL8l9nyjTA26rppgLz7+m64cQkoUkK7u+vv7x48fx8fGOjo6enp63bt2qqqp6BSf5NZ0xcuW77xr++Wt1UMmLjvTXYZT9sdaAyABdVgfihamK4T2N4N7xiIzWC7jjHtlA54cQEockKzsvL+9ui1S3/6T+rezvvxf886gIqOLFxAToM8r+QLNfRIAOqwPxYgj4cRcsQ9yOeSWYBLGaCEJikGRlq6ioTGnMpEmTRo8e/cUXX4wfP/769eui5nZLc2WzZjeo5sVyA08xyu6v/l64vzarA/FiJFoEb55jJbvK/q4n3a+AkFgkWdmZmZlXGnP58uXTp0/b2tquX7/+FZyXzSg7vlu35F9+Yc1uUMvjxAUZMcruo94r1P8kqwPxAiTbhW2dY6W4zuGiVUg1nX9NSC6SrOzmYX5Ks3jx4rS0NNGgdotI2d27p4wcyZrdoI7POR1kwii7t1qPED8tVgfiv9LA49/zjHRX8j5vQaf0ERKOJCsbBfW5xpw9ezY+Pt7Ly2v69OnXrl0TNbdbGGUnvPNO6h9/sGY3qOdzzwWbMcruodo92E+T1YF4fiqiuEm2oW7HvPMCY4pCY+n3MoTEI8nKDgoK0v5fNDU15eXlT548+fDhQ1Fzu+VvZY8ezZrdoIHPuxBiySj7bdVuQX4arA7Ec1IWyT1rHiy90u74Rsf8oH/crY0gJBWJVbZAIEBZbWho6Ojo6OTk5OHhweFwysrKXtlJfm0dyxbE8S+HWjPK7qLSJcBPndWBeB6Ev0e3C5NfY7/8L/NbrhG1HLrfLvFGILHKbmhosLOzi42NFT1/hWGUHdelS9Lw4azZzZAUatulUdnA31et5YmAxLNp4PFLwjiex73XTDW/4x5B1+cj3hwkU9l1dXWopo8dO+bm5iYa9ArDKPvMe+/dmDGDNbsZksPsuqu+zSjbx0elnk9HYP8bBcGxRaGxNRxeWST7l0oEIdlIoLKrqqrWrVuXkZFhYWGxfft2DQ2N4ODgqP+luLhY1K/dwij7bN++N+fMYc1uhpQw+55q7zCFtqfP8Vo+eec/kGgRvGuhteE+18poHtXXxJuGBCq7oqJiypQp169fh7I3btwIa8vIyMjJycnKyuJvVlaWqF+7RaTs/v1vLVzImt0MqWH276n36qLSBcp291aq4dG/zp6XRIuQzbOtjq5zuGwT2sCjA0rEG4dkKnv06NG2trZGRkbH/hklJaXs7GxRv3YLo+xz779/e+lS1uxmSA1z6Kf+HqNsFy/FKl4MqwPRKnC01nZn9a1O58yC6fxr4s1EApVdXl7+008/HT58+OT/ovW/4PGDBw9E/dotImV/+OGdlStZs5vharjDBxp9GWU7eipUcKNYHQgWEPRdj8iySG6gum+aU3hlNN3CkXhDkcwq+/fff0eJ7dci/v7+hYWFon7tFkbZiR99lLl2LWt2M1wLd/xEc0DXRmXbe8qXc+maGM+iPIp73jJEe6dLtk9UHZcOhhBvNJJ8LFv0/JVHpOxPPrm3cSNrdjNcD3ccpPUho2wbD9mnnAhWB6KJimheokWI/Br7pZPM7nhE0vnXxBuOBCq7srJy8eLFN2/eFD1/5WGUfX7gwKxt21izm+F6uNOX2p90Ve0KZVt6SJdwwlkdiCauOYYrrnOAr2+7RQhatBLEm4YEKlsgEDx+/Li2tlb0/JVHpOxBg+7v2MGa3Qxp4U5DtAcyyjZzP1TECWV1IJq47hgedtLvpks43a+AIIAEKhupr69/BT9MbyuMsi989ln27t2s2c1wI8L5W53PuzUq28TtQGFsCKsDARItQoI1/ZLtwsoiuXQ/MIJgkExlv96IlP3559l797JmN0N6hMsPul8xyjZ0kyqIpbuosDlvGbJ1rpXJfrf73nQ6DUH8DSlb/BEp+4svcvbvZ81uBij7Z70hjLINXPfkxQSyOrzJ1HL4153Ct8+3PrLW4Yzw/Gs6n48g/oaULf4wyr745ZcPDh1izW6GmxGuI/SHMcrWddmVGxPA6vAmUxXDS7YLO7reofF+BeRrgvgHpGzxR6Tsr79+JCvLmt0MtyJcR+l/1021G5R90nnHw2g/Voc3k/IoboZbxHmL4Ioo7iP/aLofGEG0hJQt/jDKvjRkyOOjR1mzmyEjwm2swY9vNypb03nbg2hfVoc3kMpo3nnLEJUtTnKr7VlNBEE0QcoWfxhlXx42LP/ECdbsZsiIdBt/6idG2WpOW7OjfFgd3jQE/LgrNqEHl9kumSg8/5rVShBEE6Rs8YdR9pXvvivQbP2+jrcj3f8y+o1R9gnHzVlR3qwObxpFobHqW522zbVKcwqn6/MRxDMgZYs/ImX/8EOhri5rdjPciXSfZvw7o+zjjhvvRXmxOrxR5AXGlIZz7ntHob6uoePXBPFMSNniD6PspOHDiw0NWbOb4U6kx0yT0cyNaY46rM+M9GB1eEMQ8OMuWIYcXmkXaxBQHcOj+pog/hVStvgjUvZPPxUbG7NmN0NmpOdck7GMshUc1t6OdGd1eBOojuVdtgndMd9afo39ZetQukQfQTwPpGzxh1F28s8/l5iasmY3w90oz4WmfzLKlrNfnRHhxurwJvAkJFZtq9PeRTaJdP41QTw3pGzxR6TsX38ttbJizW6Ge1FeS8wmMsqWtl95M8KF1UGyqYzmwde5/tGa251vOIVXx9Dxa4J4XkjZ4g+j7JQRI8rsWz/FOCvKa4X55HcalX3IbvmNCGdWBwmmKoZ3xTY0QN03P4hun0YQ/xlStvgjUvbvv5c7OrJmN0NWlPdqi2mMsvfbLrke3no3yaOBx79sHbp/qfD86zsedC8egvjPkLLFH0bZqSNHVri0fsTjfpT3eosZ76h2h7L32iy6FubA6iCpXLQK2bnAets866v2YfV0fghB/HdI2eKPSNmjRlW4tf5/xewon82Ws99REyp7l/XC1LA35SfaiRbBPid8zluE0PVDCOLFIGWLP38r2731s/dyon22W8/r0ajsHdbzk8PsWB0kjAa+sL6OMwpKcwrPC4yh80MI4oUhZYs/ImX/8UeVV+s/a8yJ9t1ls4BR9laruVdCbVkdJInqGN4Vm9BdC6zdjnnn+tO/HAnipSBliz+Msq+OGVPt1/pVVR9E+0nZLmaUvclq9uVQa1YHiaEqhpdiF3Zgqe22uVaJFnS/AoJ4WUjZ4o9I2WPH1gS2fruZh9F+B22X9VB7B8reYDnzYmjrp29LAEWhsSFafovGm6Xah9H51wTx8pCyxR9G2df+/LM2OJg1uxkexvhL263s2ajstRbTzodYsDpIANWxvCfBwt/LFIdx7npE1tPv0QlCHJCyxZ8mZdeFhrJmN8OjGH95+zWMsldZTD0XYs7q0NkR8OOSbEOPb3Q8tv5NOeWcIF4NpGzxh1H29fHj68PDWbObITcmQNFhPaPs5eZ/nQk2Y3Xo7FyyCtm10Gb3QpvrjmGsJoIgXgZStvgjUvaECfWRrf/A73FMwHHHjYyyl5pNPB1swurQqSkO5RxZ6yC/2v6cWXAVHb8mCLFCyhZ/RMqeNKkhNpY1uxnyYgJVnTb3alT2ItMJ8UGtX6O10yHgx+UHxTyN4PIMgy5bh1bQ+SEEIW5I2eKPSNl//SXgtV5j5sUEaTpv66XWA8peYPonP8iI1aEzgoI62S5MbYtTpkdkTSyPfo9OEO0BKVv8ESq7S5e0KVME/Na1lR8bpO2yk1H2XJOxvMDWb17TiUBBfck69PAKu02zLNOcwun8EIJoJ0jZ4g+UfaBr1xvTprHmdRMFscH6rnsYZc82+YMTeIrVodNx1zNSc7vz/HGmV+3D6PohBNF+kLLFH6Gy3377xowZrHndBJRt6CbVu1HZM4xHxQTqszp0Imo5/OoYXpJt6LH1jhmuEXQ8hCDaFVK2+CNUdvfu6bNmseZ1E09iQ0zdDzLKnmr0e3SAHqtDJ+K6Y/h5i5CikNgGHl/AZ7cSBCFeSNniT9OxbNa8bqIwNsTC/XBvtZ5Q9l9Gv0UG6LA6dBYuWYfuWWSjsMYhxzea1UQQRHtAyv47d+/eNTc337Jli7S0dGhoqGhoY27cuGFiYrJ9+/bdu3fLycndvHmzpqZG1NYiQmV37Zo2dSprXjdRxAm18ZDtrS5U9kTDX8L9tVkdOgWXrIW/l5FbZX/GlM6/JohXBClbFIFAYGVlpaioaG1tferUqR07djx69Ki+vp5pTUtL8/Lycm8MlK2rq5uVlcU0tQyUfbBnz5vz5rHmdRNFnDAHT4V3G5X956nhof4nWR06BbEGAaYH3OKMgiqi6PxrgnhFkLJFycvLk5GRMTY2Lisru3r16rZt24KCgioqKpjWoqKi3Nzcurq6ysrK4ODg1atXJyUlMU0tA2Uffu+920uWsOZ1E8WcMBcvRUbZYwx+CPHTZHXoyKCgTnMKT7ELu+YYft87inxNEK8SUrYoV65cUVBQ8Pb2xuOcnBwDA4Pjx4/D1ExrU6DsgICAzZs3p6amiga1iFDZffrcXrqUNa+bKOGEu3srMcoebfB9kJ8Gq0OHpSKad9kmVHGdg/cJn+IwDquVIIj2hpQtCp/Ph6NDQkLwGAW1ra2tlJTUkydPmFYmtbW1d+7c2bp1q6mp6YMHD0RDG9PQ0JCfn5+enn79+vUDq1dL9+17Z/ly1rxuopQT7u2jzCj7d/1vA3zVWR06JtWxwt83Hl3vMGeMyVUHOv+aIF4DpGxR4uPjlZSUgoOD8RjKtra23r9/f3NlQ8rZ2dmGhoazZ8++d+8enooaGoPq293dfdmyZTNmzBj5+efS/frdWbmSNa+bKOVG+PuqvafeC8r+Tf8bP181VoeOyX3vqFNSrnPHmNyi868J4jVByhbl2rVrcnJy0C4e379//+TJk5qamsXFxUwrwpTeU6ZMga9RbouGNktdXV1VVRXcrbRzp3T//s9Q9lNuRJCfBqPsn/WG+PiqsDp0QODoymjuI7/omy4RDXT+NUG8JkjZolRUVBw+fBimRmV98eLFRYsWJSYmVv9vvhQUFDg4OKxbty4hIQG+FggEzPBWo7J374Fu3W5Mn86a102UcSPD/E8yyv5R92tvnxOsDh2NS9ahWjucPY5713F5tRyqrwnitUHKFgUWDg0NlZaW3rRp044dO5SUlJ4+fRoQEMDn8/Pz8yMjIyHx4cOHb9++ff/+/SoqKpmZmaJXtsi/njFSzo2MCtBllP2d7peePsqsDh2Ky42/l1Ha4Jhk2/pNdgiCeGWQsv8OSum4uDhnZ2dfX1/mhJCkpKSbN2/C3enp6Z6envb29h6NCQ4OhseZV7UMlC3z/vuZa9ey5nUT5dyo2AADRtnDdD5391ZideggCPhxpeEcmVV28mvs+YaBdD4fQbx2SNnij1DZH3xwd9061rxuooIbzQ806qPeG8oeoj3QzfsYq0NHoI7LLwrllIRzjPe7JZoHl5OvCaIDQMoWf4TK/vDDu+vXs+Z1E5W86NNBJoyyvzr5iYvXUVaH105lNO+6U7jbMW9U1tWxvAY6P4QgOgakbPFHpOwNG1jzuokqXkxisDmj7C9OfuTkdYTV4fVSFcO7YhN6dL3D3LEmjwNi6Hw+gug4kLLFH6GyP/ro7saNrHndBJR9McSSUfYgrQ8dPBVYHV4v1xzDFdbYz/7D5JZLBNXXBNGhIGWLP1C27Ecf3Wtb2dW82KRQm76Nyv5E6317T3lWh9eIgM8P1fI7sckx1T6M6muC6GiQssUfobI//vjepk2sed1EDS82NcyeUfZHmv1tPeVYHV4XVx3Cb7lGPPSLfuQfTb9HJ4gOCClb/BEpe/Nm1rxuopbPSQt36qvxLpT9vkZfaw8ZVofXwhWb0H1LbR0UPAuCY1lNBEF0EEjZ4o9Q2Z98krVlC2teNwFl34pwZZTdX+M9Kw9pVodXTC2Xn2Iftm+JrcxKO55hIJ3PRxAdFlK2+POvyq7jc+9EuvdrVHYf9d4WHodZHV4xVTE8fzXfw8vt4o2DyNcE0ZEhZYs/QmV/+mnW1q2sed1EPZ+XFeXFKLu3ek8z94OsDq+MymjeQ7/ou55R5y1C7rhH4imrA0EQHQpStvgDZctB2du2seZ1Ew18Xk6UL6PsnmrvmLgfYHV4NTDXvzY94Oal7EMnhxBEp4CULf48h7L5udEB/TTeg7K7q3YzdtvP6vAKEPDjUh3CZFfZzxxtfNMlHE9ZHQiC6ICQssUfobIHDry/fTtrXjchiOMXxAT1b1R2V5Wuhm77WB1eAbdcI+TX2K+fYZlkE1rHpRKbIDoHpGzxR6jsQYPu79jBmtfNKYwNYZQNDFz3slrbm/rGSz5dsg49Zx5cFUPHrwmi00DKFn9Eyt65kzWvm1PMCXtfsy+jbD3XPQ38V+fNJNtQs4NuF61CKqO59P9GguhckLLFn389MAKg7I81+3dR6QJl67jsquO/ilPrajn8qw5hB5baMr9HZ7USBNHxIWWLP1C2zIABmWvWsOZ1c0o4YQO1PujaqOyTLjuqee3+g8MGHj83IEZ2ld22uVbxRnT+NUF0SkjZ4g+UffjddzMWL2bN6+aUcMK/OPlxN9WuULam8/ZKXjSrg3ip4/LLIriZHpFrp1lcEh4SoeMhBNEpIWWLP1C2/OefZ+/Zw5rXzYGyv9Ye2E21G5St7ry1nBvJ6iBGaji8LK+oKzah9Tx+VQyXzucjiM4LKVv8ESr7iy/+VdnDdD5/u1HZqk5bnnIjWB3EyDXHcJmVdvPHmRaFclhNBEF0LkjZ4s/zKLuUE/697pfdVd+GspUdNxVz2uufgcl2YQeW2W6aZXnJKoTuV0AQnR1StvjznMr+SXfwO43KVnLcUBgbwuogFgR8vruSt/ZO53jjIDp+TRASAClb/BEpe++zfiBTyon4VW/oO6rdoeyjDusLYoNZHV6ea47hD/2iU+3D0p0j6PwQgpAMSNniz3MpmxsxUv/bHmpCZR9xWJsXE8jq8DLUcvjC49er7LmnAksj6Pg1QUgOpGzxh1F2zjOV/ZQbMebUDz3V3oGy5e1X58b4szq8MFUxPFTW0ivsNs2yjDcOqqD6miAkCFK2+ANlK3z5ZY6UFGteNwfK/tPwp16NypaxX/Ug2o/V4YXJC4qxl/ecPtL4oiWdf00QkgYpW/x5TmVPMvy1l1oPKPuw3YrsaB9WhxeglsNHTf3IPzrOKOimSziesjoQBNHZIWWLP8+p7KlGI3qr9YSyD9ktz4ryZnV4Aa45hjsregVq+DbwhZfDZrUSBCEBkLLFn+dRdhk3cobxqHfVhcreb7v0bpQnq8N/Jdku7NByuz2LbNKdw1lNBEFIDKRs8ec5lT3HZMx76r2gbCnbxXci3Vkd/hO33CJkVtlD2XzDwMpo+n8jQUgspGzxR6Tsfc+61wyUPd/0T0bZe2wW3YpwZXV4fuq5/Ae+0V7KPhyDwPJI8jVBSDKkbPFHqOyvvvpXZS82ndBHvTeUvct6QXq4C6vD81DD4aU5hUPWD/yiS8I4dD4fQUg8pGzxh1H2g/3PuglvOTdyudlffRuVvcN6Xlq4E6vDv1IVw7vqECa/Rng85I5HO14IkCCIjgMpW/x5PmVHrTKf0k/jXSh7q9Xca2GOrA7PpoHPh6ZPbHacMsLoklUI3b+RIN4QSNniz3Mqe53l9H6Nd+zdbDU7Ncye1eEZNPDjKqK4UboBc8aYpDuH0/3RCeLNgZQt/jynsjdazWJusr7RclZSqC2rwzPID4rN9Y+pieUVhcbS9VQJ4o2ClC3+PI+yK7hR26zmDtDoA2Wvt5xxOdSG1aEtku3C5NbYq211Kovg0O9lCOJNg5Qt/giV/fXXDw4cYM3r5kDZO60XvN+o7LUW0y6GWLE6tArze5nDK+zijALp9+gE8QZCyhZ/nk/Z0XttFr+v2RfKXm0x9XyIBatDSwR8vo2sh8pmp2i9gDI6/5og3khI2eLPcymbF73fdtkHjcpeYT75bLAZq0NzajjC80OKw2IjdPxRaJOvCeKNhZQt/kDZR6DsgwdZ87o5lbzow3YrPtTsB2UvM5t0OsiE1aGJqhjedcdw9a1O+Esn8xHEGw4pW/x5TmXL2a/+qFHZS8wmxgcZszowVMfyrjmGn9js+NevRldsQknZBPGGQ8oWf4TKHjz4X5V9xH7tR5r9oexFpuN5gYasDgy33SNVtzhN+tWQzr8mCAKQssWf51S2kuOGj7UGQNkLTP/kBJ5idWC47x0Vqx+A+rqefE0QBCm7PfI8yq7ixag4bf6kUdlzTcbGBOizOqTYhwWq+541C34awaXjIQRBMJCyxR+Rsg8dYs3r5kDZGs7bPtV6H8qebfJHVIBu81b4+vAKO60dzjddIpoPJwjiDYeULf48p7JPuuxklD3TeHSEvw4zvIEfd9czSnaVPZQdpUvnXxME8Q9I2eIPo+yH/6ZsPZfdA7U+gLKnG48M9T+JgYLGW+6mO0cob3KMNQgkXxMEwYKULf4IlT1kyMPDh1nzujnVvBhDV6lBjcqeajQi2E+zKoaX5RV11iy4IoqbFxhTGU3HrwmCYEPKFn+eT9mxJm4HPjv5IZQ92eg3b0+1607hmtud10yzKA3nsDoTBEEwkLLFn+dUtrn7IUbZEw1/0TQ4cnS9w8SfDW84hdP5fARBtAUpW/x5TmVbe8h8fvIjKHuM3k8Ltu9ZO93iomUI/V6GIIhnQMoWf55T2fYecl80KnuswXBDiyPXHcPp+DVBEM+GlC3+iJQtLc2a182p4cU6eR358uTHUPYovR9cPJTo/jIEQfwrpGzxR6jsoUOfrewqbqyOuexAdWGVPVLvO2fPY6wOBEEQLSFliz9QtuLQoY9kZFjzujnlsTFb5fcNOC789+NIg29dvI6yOhAEQbSElC3+PFvZ1bG8JyGxT8KitskeZJQ9Qn+Yk9cRVjeCIIiWkLLFn2coW/jjRpcIT2WfDI9Qfx/1r7U/hbJ/0R9q7ynP6kkQBNESUrb48wxlw9eK6xwm/Gx41SkkyEdrsPZAKPsnvcG2nnKsngRBEC0hZYs/bSn7mmO4/Br71VMtzpkHV8TGRvjrDNEeBGX/qPuVtcezDnwTBEEwkLLFH6Gyv/nmkawsa16nOYX7q/lyDAIro7m1fG5MgP5QHaGyv9f90tL9WSdxEwRBMJCyxZ+Wyr7qEMY9FZhqH5YXGMNcn6+Oz+UFnhqq8xmU/a3O52buz7rsH0EQBAMpW/xprux6Lv+Gc7jcanvLw+533COb5juUnRBk/E2jsvHXxO1AUxNBEERbkLLFH5Gy5eRqOfy7nlGK6xw2zbKM0fvH/Qqg7LPBZsN0PoeyB2sPNHLb19REEATRFqTsv1NRUXH37t2kpKRr167l5uaKhjZGIBCUlJTcuHEDrampqU+fPm1oaBC1tYhQ2cOGocoui+DGGQXNHWuSYBJUEfWP+xVA2RdCLL9tVPZX2p+ccpVq3koQBNEqpOy/c/78+YMHD44dO3bu3LmnTp2qq6uDqZmmqqqqqKioxYsXjxs3bsqUKVwut6ysjGlqGShb/ptvb+6TzQuIqYrmZXlF1nDY13uq53OvhNp8p/sFlP3FyY/1XfewOhAEQbSElC0KqmZZWVkVFZXs7Oz4+PiZM2eipq6trWVaL1++rKWlpaOjU1xc7O7uvmPHDlTiTFPLQNmrPv1JbuTBnfOtBXy+gM+e6aCez0sJs2OU/dnJD3VcdrE6EARBtISULUpGRgaU7eDggOL6zp07x44dMzU1LS0tZVp9fHwUFBRQhsPsRUVFEyZMSEhIYJpaRnHvgcE9/lwzVue8RQhrdjcBZV8Pd/pe90soe6DWBydddrI6EARBtISULcrZs2eh6YCAADx++PChubm5jIxMYWEh02pjY3Po0KGcnBw8FggEkyZNioqKqqmpYVpZObRXdvonY3w2qrKOXzcHyk4Pd/lB9yso+2PNAZrO21kdCIIgWkLKFoXP5x8/fjwkJASPc3NzbW1tpaSknjx5wrRaWloeOHAgPz+feTp16tTw8PCqqirmKQJ9Q/rGxsYnT54c/+dfIwcOkZk+9+SuXW2htWvnkW3rfl3yzZfzP/lu0ZfL105mdeiwaO7YobRx48mdO1nDJQnVrVtVt2xhDZQYtHbuxBLEX9ZwiUFlyxa1rVtZAyWGgytWjBg+vJqUff78eUVFRT8/Pzx+8OAB5KugoFBUVMS02tvbHz58OCsrC48bGhomTJgQHR3ddKQbweMLFy5YW1sbGRnNmzdvycKFRrq6EomSgsKvw4framiwhksSq5ctW7ZoEWugxKCqpPTjt99qnjjBGi4xLJ43b92qVayBEoOCtPTvI0aQsv/v/v37MjIyVlZWmBfp6emHDh1ycnJqOi0kODgYQudwOPX19RD3tGnTUFMzTS2DQhvGFz2RuNy8eXPy5MlPnz4VPZfEMHtLoicSl+zs7D/++CMvL0/0XOKiqqqK3WLRE4nL1atXZ86c2dZR2Q6SV6Fs1M7q6urHjx+Pj4/38fFBmZWZmYlyu7CwEHPn+vXrhoaGysrKly9fxvosJycHc4le2SKk7M4eUnanDin7tedVKBu5du2avLz8yJEj58yZY2triyEWFhbQd05OTl1dHVS+YsWKUaNGocS+cOFCZWUl86qW0dfXNzc3Fz2RuNy6dQvfmGecli4BweLDQhQ9kbjg+zxx4kQJVraWlhaz/kpkoKkFCxaQsoWpr6+HiFE/lpeXM3Okurq6traW+aEjrI3haIWt0LPpVzYtExcX94zDJp09+fn5NjY2HfxQ2ksGiw8LUfRE4lJcXIxtkgRvdDkczvnz50VPJC65ubn29vZQkOh5h8wrUra4ArNXVFSInkhcsOkqLCx8xhZLAoLFh4UoeiJxwdr+5MkTphCRyGBrJNkrYFFRUQdfATuZsikUCuVNTsdVNqqV2NhYTU1NHR2dmJiY5qUZNoOPHz92cHBAq4GBQWpqavPzuDtFUI7dunXLzMwMk+Dm5paRkSFqaMzly5etra0x4YaGhkFBQU2/FO1EycvLCwsLw9SdOnWKz+c3P2uTSXV1dVJS0tGjRy9cuNDpCjfm3+b47mloaHh7e9+7d0/U8L9ggTo7O2tra5uYmODbKxraqVJQUBAVFYUlqKenx+Fwmv+HCZOPb6+pqam+vr6urm58fHxxcbGorTMkPz/fxcXl+PHjJ0+ebLnsoBcsMiMjI7SGhoY2/YKkg6SDKhtzjcvlysvLKygoKCkpbd26tfllSfD9iIiI2Lx5s4qKipSUFNaKZ5xk0jGDHUx8Y+Tk5DB1Bw4ccHR0bH4AFLLDegIdKDYmPDxc1NBJgsWH7zoW35EjRyDlXbt2ZWdnY69T1Nx4EtHDhw8xaYMHD3Z1dW36KWxnCSoGfDOx+LAQ8Q308fFpXlJgE4vv5LFjx9TU1CC14OBgUUPnCZZgdHQ0swJiQrZs2XL79u2mFfD+/fu2trZ79uzBNgmLGPOhcx3gzsnJ0dLS2rFjx8iRI8+dOyca2hjUUllZWdLS0phwpOk3gB0nHVTZ8JexsTFchk0cNokLFy709fVt2tyhwDlx4gRqHGztL126tHv3bubX8J0lkFdmZuYff/yBrwsKTHNzc4j76tWroubGGg0Ty1TiWOf3798vaugkefr0KWyFrztcjFUdiw8boeb7CiUlJShksNGdOHEilmznUjbMdeXKFaztKCNQe2LlR9LT05lWyO706dN79+5lrruAdMZj91gBUTTIysrie4iN6+zZs7GKNf38LTEx8dChQ9jNxaYXE75ixQrsajBNnSJY6eAQrHHTpk1jKRsL1MnJCWscWrFhxkyQkZHBQhQ1d4B0UGXDWTo6OtjnwmN847Fzra6ujpWfaY2Li9u4cWNycjK+MfDaunXrrKysOtH/fGA0rM9LlixhdsqwP4FdsFa/9CgHrK2tUaWKnneSpKWlYQcICwWPsX+N5QiDY81nWrHFSkpKQnWGCZ83b56fn1/nUjbGFv7C4mMUhv0JKBt/mVZ8D5WVlbGqQ9zYMGOqn3HSaocNdls1NTWNjIzwGNtXTA5WwKZjCFj1UHpjCDZamBUotzvdwR98CTE5M2bMYCkbE7tz507s9TKX9be3t8fax7rE/+tNB1X25cuXoWzYCo/xjXd3d8cGHyJgWrHLhhUGKwMqGjzdvn07CtVOVMtgVccU4Zvx4MEDPL148SKkZmNjw7Q2BSs/h8PBPhq+N6JBnSRYDbA+e3h44HFxcTFWgIMHDzLXJEBQvHh6ejLFC5Zjp1M2VmBMEfbtsHrjaUJCAuzGTCyCGmLr1q1SUlLLly+fPHkyyu0zZ84wX9ROFHwnsR3CZOIxKm4XFxd8D7HPx7SiSkXN8euvvw4bNmzQoEEWFhZYpkxTZ0lbysa6ie9kUFAQsz328vLCVxdbJqa1I4SU/RrynMpOSUlBLXP48OFO9+/HZygbiwwFKdZ/7EhJqrI3bdqE6cJk3rlzB4Xq+vXrO90vWp+tbNTgxsbGhoaG+fn5PB5vzZo1qLU712aJlC3m4MvR/MAIvhzYs25+YGTDhg0wWvMDI3jAtHb8YAWOjIxcunQpYzHmwAi+HEwrExgN+9dHjx7FhqrT1WjXr19XUVFhtrgFBQXYrVZVVWUOjFRXV2Pj9MMPP8yZM2f+/Pmffvrp6NGjnZ2d0a3xpZ0gzIERLD5mrQ4LC2MdGIHdsMXCgoMXAgMDt2zZ0nSku7MEI4ztEHNxCFQM+vr6WAGbDoygZsJEMd9MVNxYAbEj2LmO/zzjwMiOHTucnJyYgyEODg7YNtOBkX8PpIbyBNs3rB5YmbHdg9Ga1upr164pKSnh+4QyLSkpadeuXcxlAjtLamtrYeQ//vjj/Pnz+KJbWloeP348OTlZ1Nx4VXE4DsrGPnXL0+M6fvC9Z8YfJTYqTSy+4OBgZl8BW1ZselHFwHHQ2bhx4xQVFbFT1YlO08S37tKlSyNHjkSxidFGbQG7Ne0CQtmenp5YoJgobJ+wO7V69Wpmd6oTBQsLBZOcnBw2S/g2YuPa/L/EUPbatWtRX0PZjx8/xtYLjoMEmdZOkbaUjS0QNj+HDh2CZPLy8rCtwga4Q/0guYMqG1+FmJgYWVlZrPYoWFBTo3CD4LC/hvmIr05ISMjmzZu1tbUxc0+dOtW0wnSWYJskLy+PIhpTd+DAAVtbW4iMz+fDApCanZ3dxIkTly1bhurb3Ny8050lhsUHHWOFP9GYrVu3Yn8Cmx9YrPlZrp30wAim7tGjR/jiHTt2TENDQ0pKCo5GWRofH4/tE1oxsRiO/QzYHBskuK/T/W4AUxEeHo4VEIsPW9/169fjmwm7YVuVn5+Px5g0fG+x+4vJhNQ612UkMAnYscMK+M0336CIxnf1woUL0Askg7UPCxHFIopCbHfx19/fX/SyjpEOqmwEsxVfGswyfCew71lWVpaYmIjVHmsLChn8ZYpTrBsQQaf7pzy+GdjMYBuO1R6CvnHjBqYoMjIS23Y0ubm5YcIRCB2bJRhB9LLOE0xOQEAAJgEVKIoy2BlGwzqPxSrq0VjpYCFeuXKl050Gh7IrJSUFk4bFh5UfKzk0jSKDOQAC36ECZXSGkg2CY17VuYLaCIURlqCamhrWRKxi8DJWwNzcXDzG5OObiSZ8RVFqMMeIOkswCfjiQdnY3MrIyGBPKCEhAd9PZk8Xiw/CQbWEehG+bv6N7QjpuMqmUCgUCiukbAqFQuk0IWVTKBRKpwkpm0KhUDpNSNkUCoXSaULKplAolE4TUjblZXPu3LmTJ09a/C93795t9ec/9vb2L/bD36ysLA8PDyMjI7y5paVlcnLy858UyOVyb926VVlZmZOT0/zT/f3979+/L3ryH3P79m11dXWMjJmZmbOzc0ZGRlu/IikrK8vOzsZHi55TKC8dUjblZQOZDh8+/Pj/AjPWtHaxykmTJgW80DVyExIS1q9fv2bNGrz5rl271NTUrl27Jmr7t0DNKSkpUHxISIihoaFo6P/9n7W1NVQrevIfExUVNXDgQIzMsWPHMGKOjo5t/bgxMzMzMDCw6bfsFMrLh5RNedmYmJisXr1a9KTxF9sotJOSki5fvpyWltZ0vxJG2fX19Y8fP2ZakcLCQgyBUmG3K1euoIJGDYshzEuYQNnS0tIRERF4jA7MxdNLS0vh3EuXLuGtUMaizsXn4gE6YODVq1efPHkiEAjS09PxcY8ePYJe586di854FXriAT46NzcXJTxzzSYMRBOGVFdX47XY8GB8oHsMQVPjiIgSHR39+++/M48NDAwUFBTi4+OxY9E0XfiLd8AoocY/ePCgrKws3ic/Px99UNqnpqbinTFizFWlKJT/FFI25WUDZS9evBhqQyAmGFBbW3vRokUzZ87cuHFjUFAQ041RNlxmb2+PplmzZs2ePTs2NhbyTUxMPHr06Pz586HjAwcOwLDNrd1c2dAcXojClsfjSUlJTZs2DSJG3Q1dYtugr6+Pp3hzFL+ohaFajIBrYzB80KBBGCtFRUVI+Y8//sAburu747XMj62x2di6dauLiws2Hn5+ftu3b8dLlixZYmxsjCmC/RvHRZjmyvby8lJSUuJwOAUFBdiQYPwxUYiPjw9GCa/97rvvfvjhhxUrVqDMv3fvHubMqlWr8M7YXcAINH9bCuV5QsqmvGxOnTo1YMCAESNGjBw5Eu5GmVxWVlZVVQUJwtHz5s1jxMQoG/7V0dHx8PCAlGtqamBDVLvMr/YrKyuhy5MnT2poaDDXkGLSXNn4u2zZMlNTU2VlZTk5OXTDx/3yyy+RkZFoUlVV9ff3h6lRzzLvwCg7Ly8PI7lu3To0IRjOKBvFr66uLjY5GEM4d/jw4aiR8VaonfEXZfL169cnTJiAurj5hYGalI1JwJbGzMzszp07mBZMNf5iIIruffv24S+qaby/oaEhPhQfoaenZ2VlhVq+oqKCuVwUHjDvSaE8Z0jZlJcNlIfqMqcxKLRhXvhoy5Ytc+bMgRnhcdgQzmKUDWFBrCiobW1tmaPeFy5cQDELCULuKJAnTpy4Z8+e5tesgLKXL1+Ot0IH9EQhDONramriU/C26CkjIwOJnz59GqJcuXKlk5MT3plRc5Oy4U2U3swbIoyyYUwIFx+HzUZoaCgcim2GkZHRjBkzMBr4OFTxP/74Iwpk5uAJE9i8X79+aMWboFt4eDimAtsnjOe2bdtQYk+ePBkbMLxhUlISNM1cwhTBuGEDgDdHlf3nn3/iHe7evUuFNuU/hZRNedmwjmWfPXtWUVHR3NwcD2BSSAoqb1I2hJ6Wlubt7Y3yc8OGDcHBwWFhYVJSUnh6/vz5xMRE/IVwm59zAhXCqnAumlC35ufnw7Yoxpkr7qKaRsWNIhobg+TkZE9PTxTphw4d8vX1RcH7bGVDlxgllPAYB3jf0tLyyZMn2AlA/e7l5cWMz7lz5/Dy5ueEoMr+7rvv0BQbG4sJ19LSunnzJkYMmyJMBfpD8XA93pml7AULFujr63M4HLwzumFsqcqm/NeQsikvG5ayfXx8oD94DaWrjY0N6s2srCzYk1E2JFVcXIy6G85CUQzNwch4ALmjkoXZ4fTHjx+jv+jt/nlghElqaqqamhpsiPL24cOHqPFRdzP/SESZz+Vyjx8/fvjwYXi/SdmonbGFYEpvhFE2HsCqcPTevXv/+uuvq1evVldX462OHj2KbQk6w9R4bVVVVfNauOnACFqxYdixYwf8HhMTgzcJDAxET/zFu2Fi8YYGBgYYT+aFmAptbW1skPAYVTkmkxlOoTx/SNmUlw1L2SghZWVljx07ZmFhgfJ5zJgxzZV99+7doKAgqNzMzAytDg4O2dnZ1tbW6I/3wQM7OzsUsHCx6O1aUzak7+zsjFIatTzqa5TPkDiKd39/f7wDCmq8G5qg1CZlo/qGsm1tbZmD1E3KRlnt4uIyatSo2bNnl5SUQLio1uF3WNuqMXgf1PXNNyHN//2IbQ/GDYU26nRsfrCtwvjjwejRozGx2JxgSjdt2oRxgKnxQiUlJewfYCBmDrZtzJtQKM8fUjblZRMVFQXbip401o+QoKKiooaGhr6+PupK5lg2RHbx4kUom7m6PDqoq6vD5qiFc3Jy3N3dMQTDUT5fuXKlubLT09OhvORmd+1BIFZ7e3s5OTlIMDw8HB+amZkJFSooKBw5csTY2BjvDP/CjKdPny4tLYUxIV+IGErFJ+JVzBuiz+XLlw8cOICNB/OhGILSG1OEN8e7Yfyh9abyHEHtjFEVPWm8wzqqbLwJj8eDvtHk5uaGbQYmFm949uxZfCimHXsVZWVl6IzZwow2ynnRW1Aozx1SNoVCoXSakLIpFAql04SUTaFQKJ0mpGwKhULpNCFlUygUSqcJKZtCoVA6TUjZFAqF0mlCyqZQKJROE1I2hUKhdJqQsikUCqXThJRNoVAonST/93//Dx33XGvbRzjpAAAAAElFTkSuQmCC"}}},{"cell_type":"code","source":"def prgauc_metric(y_true, y_prob, sample_weight = None, plot = False):\n    '''Function to obtain the Precision Recall Gain Area under the Curve (PRG AUC). The baseline is computed from\n    the mean of y_true. Points with non-infinite recall or precision are plotted (TP != 0) (our decision-making).\n    Only points with the recall gain in the 1st quadrant (recall gain >= 0) are plotted (according to the reference).\n    Additionally, we do the manual inclusion of the [recall gain = 0, precision gain = max(precision gain)], to make\n    the curve artificially touch the y-axis to obtain a better PRG AUC approximation (our decision-making).\n    References:\n    https://research-information.bris.ac.uk/ws/portalfiles/portal/72164009/5867_precision_recall_gain_curves_pr_analysis_done_right.pdf\n    \n    The function computes the precision gain (_precG_) and the recall gain (_recG_), and computes the area under the \n    precision_gain(recall_gain). The baseline of this metric (random model) is 0.50 (similar interpretation to the \n    ROC AUC). According to the reference, the PRGAUC is prefered over the typical PR AUC since it has better\n    interpretability (PR AUC has reference the y_true avg., making the metric harder to evaluate different datasets\n    with different baselines). Like the PR AUC, the PRG AUC is better than the regular ROC AUC for imbalanced\n    datasets (ROC AUC gives equal weight to the positive and negative class, given the impression the model is\n    performing well even if the identification of the positive class is not so good). \n    \n    _precG_ = (precision - y_true avg.)/[(1 - y_true avg.)*precision] = 1 - (y_true avg.*FP)/[(1 - y_true avg.)*TP]\n    \n    _recG_ = (recall - y_true avg.)/[(1 - y_true avg.)*recall] = 1 - (y_true avg.*FN)/[(1 - y_true avg.)*TP]\n    \n    Arguments of the function:\n        y_true: (list, pd.Series) vector with the true labels;\n        y_prob: (list, pd.Series) vector with the predicted probabilidades of the positive class;\n        sample_weight: (pd.DataFrame) df with sample weights, with index in concordance with the y_true (case when this func\n        is called as the scoring of sklearn.model_selection.cross_val_score;\n        plot: (boolean) if True, plot the Precision-Recall Gain curve, along with the PRG AUC;\n    Package dependencies:\n        np (numpy)\n        auc (sklearn.metrics.auc)\n        plt (matplotlib.pyplot)\n        precision_score (sklearn-metrics.precision_score)\n        recall_score (sklearn-metrics.recall_score)\n    '''\n    \n    # Build threshold vectors (more granular between 0-5% and 95-100%, 5-95% the step is +1%):\n    \n    main_thr = 0.002\n    # main_thr was obtained by weighting in the runtime and the PRG AUC approximation (assuming main_thr -> 0)\n    # gives the true value of PRG AUC. Additionally, below we do the manual inclusion of the \n    # [recall gain = 0, precision gain = max(precision gain)], to make the curve artificially touch the y-axis\n    # and obtain a better PRG AUC approximation.\n    thr_vec = [0] + list(np.arange(0.0005, 0.01, 0.0015)) + list(np.arange(0.01, 0.05, 0.005)) + list(np.arange(0.05, 0.95, main_thr)) + list(np.arange(0.95, 0.99, 0.005)) + list(np.arange(0.99, 1, 0.0015)) + [1]\n    thr_vec = [round(i, 6) for i in thr_vec]\n    # (copy this code and print thr_vec to visualize the threshold values)\n    \n    # Baseline:\n    baseline_mean = np.mean(y_true) #(between 0 - 1)\n    \n    # Lists do save the precision gain and recall gain values:\n    precG_vec, recG_vec, thr_used_vec = [], [], []\n    \n    # Prepare sample_weight:\n    if sample_weight is not None:\n        if isinstance(sample_weight, pd.DataFrame) and (isinstance(y_true, pd.DataFrame) or isinstance(y_true, pd.Series)):\n            sample_weight = sample_weight=sample_weight.loc[y_true.index.values].values.reshape(-1)\n        else: #list case (just in case)\n            sample_weight = sample_weight\n        \n    for thr in thr_vec:\n        \n        # Build predictions vector for the thr:\n        pred_vec = [(0 if prob < thr else 1) for prob in y_prob]\n        \n        FP, TP, FN = 0, 0, 0\n        for pred, y in zip(pred_vec, y_true):\n            FP += 1 if pred == 1 and y == 0 else 0 # Build FP\n            TP += 1 if pred == 1 and y == 1 else 0 # Build TP\n            FN += 1 if pred == 0 and y == 1 else 0 # Build FN\n        \n        # Exclude infinite points:\n        if (TP) > 0:\n            # Use precision_score and recall_score to take advantage of the sample_weight argument:\n            precision = precision_score(y_true, pred_vec, sample_weight=sample_weight)\n            recall = recall_score(y_true, pred_vec, sample_weight=sample_weight)\n\n            #_precG_ = 1 - (baseline_mean*FP)/((1 - baseline_mean)*TP)\n            #_recG_ = 1 - (baseline_mean*FN)/((1 - baseline_mean)*TP)\n            _precG_ = (precision - baseline_mean)/((1 - baseline_mean)*precision)\n            _recG_ = (recall - baseline_mean)/((1 - baseline_mean)*recall)\n            \n            # Append to lists:\n            precG_vec.append(_precG_)\n            recG_vec.append(_recG_)\n            thr_used_vec.append(thr)\n    \n    # Create dataframe with accepted thr, precision gain and recall gain:\n    df_metrics = pd.DataFrame({'thr':thr_used_vec, 'precG':precG_vec, 'recG':recG_vec})\n    \n    # According to the reference (see the function's help), we should only consider the points with \n    df_metrics = df_metrics[df_metrics.recG >= 0]\n    \n    # Add the point for recall = 0 (according to the used threshold step, there might be a pice of the \n    # graph missing):\n    df_metrics = pd.concat([df_metrics, pd.DataFrame({'thr':[df_metrics.thr.max()], 'precG':[df_metrics.precG.max()], 'recG':[0]})])\n    df_metrics = df_metrics.sort_values(['recG'])\n    #display(df_metrics)\n       \n    # Compute area under the curve:\n    prg_auc = auc(x = df_metrics.recG, y = df_metrics.precG)\n    \n    # Graph plot:\n    if plot == True:\n        \n        plt.plot(df_metrics.recG, df_metrics.precG, label='PRG ({}%)'.format(round(prg_auc*100, 2)))\n        plt.plot([0,1], [1,0], linestyle = '--', color='gray', label='Rand. model (50%)')\n        plt.plot([0,1,1], [1,1,0], linestyle = '--', color='green', label='Perf. model (100%)')\n        plt.plot(df_metrics.recG, df_metrics.thr, label='Thresh. used')\n        \n        plt.title(\"PRG curve + PRG AUC\")\n        plt.xlabel(\"Recall Gain\")\n        plt.ylabel(\"Precision Gain\")\n        plt.legend()\n        plt.grid()\n        plt.show()\n    \n    return prg_auc","metadata":{"execution":{"iopub.status.busy":"2023-11-02T19:31:06.138318Z","iopub.status.idle":"2023-11-02T19:31:06.138707Z","shell.execute_reply.started":"2023-11-02T19:31:06.138518Z","shell.execute_reply":"2023-11-02T19:31:06.138535Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def amex_metric(y_true: pd.DataFrame, y_pred: pd.DataFrame) -> float:\n\n    def top_four_percent_captured(y_true: pd.DataFrame, y_pred: pd.DataFrame) -> float:\n        df = (pd.concat([y_true, y_pred], axis='columns')\n              .sort_values('prediction', ascending=False))\n        df['weight'] = df['target'].apply(lambda x: 20 if x==0 else 1)\n        four_pct_cutoff = int(0.04 * df['weight'].sum())\n        df['weight_cumsum'] = df['weight'].cumsum()\n        df_cutoff = df.loc[df['weight_cumsum'] <= four_pct_cutoff]\n        return (df_cutoff['target'] == 1).sum() / (df['target'] == 1).sum()\n        \n    def weighted_gini(y_true: pd.DataFrame, y_pred: pd.DataFrame) -> float:\n        df = (pd.concat([y_true, y_pred], axis='columns')\n              .sort_values('prediction', ascending=False))\n        df['weight'] = df['target'].apply(lambda x: 20 if x==0 else 1)\n        df['random'] = (df['weight'] / df['weight'].sum()).cumsum()\n        total_pos = (df['target'] * df['weight']).sum()\n        df['cum_pos_found'] = (df['target'] * df['weight']).cumsum()\n        df['lorentz'] = df['cum_pos_found'] / total_pos\n        df['gini'] = (df['lorentz'] - df['random']) * df['weight']\n        return df['gini'].sum()\n\n    def normalized_weighted_gini(y_true: pd.DataFrame, y_pred: pd.DataFrame) -> float:\n        y_true_pred = y_true.rename(columns={'target': 'prediction'})\n        return weighted_gini(y_true, y_pred) / weighted_gini(y_true, y_true_pred)\n\n    g = normalized_weighted_gini(y_true, y_pred)\n    d = top_four_percent_captured(y_true, y_pred)\n\n    return 0.5 * (g + d)","metadata":{"execution":{"iopub.status.busy":"2023-11-02T19:31:06.140650Z","iopub.status.idle":"2023-11-02T19:31:06.141236Z","shell.execute_reply.started":"2023-11-02T19:31:06.140914Z","shell.execute_reply":"2023-11-02T19:31:06.140940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def lgbm_train(df, params, weights_01 = [20,1]):\n    \n    start = datetime.datetime.now()\n    \n    X = df.drop(columns=['target'])  # Features\n    y = df['target']  # Target variable\n    \n    # Number of splits for stratified k-fold cross-validation\n    n_splits = 5\n    auc_scores = []\n    prgauc_scores = []\n    amex_scores = []\n    \n    # Initialize StratifiedKFold with the desired number of splits\n    skf = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=42)\n\n    # Perform stratified k-fold cross-validation\n    for train_idx, valid_idx in skf.split(X, y):\n        X_train, X_valid = X.iloc[train_idx], X.iloc[valid_idx]\n        y_train, y_valid = y.iloc[train_idx], y.iloc[valid_idx]\n\n        # Create LightGBM datasets\n        train_data = lgb.Dataset(X_train, label=y_train)\n        valid_data = lgb.Dataset(X_valid, label=y_valid, reference=train_data)\n\n        # Train the LightGBM model\n        model = lgb.train(params, train_data, valid_sets=[train_data, valid_data],\n                          callbacks=[lgb.early_stopping(100, first_metric_only=True, verbose=False),\n                                    lgb.log_evaluation(period=10, show_stdv=True)])\n\n        # Calculate the AUC score (you can use other metrics)\n        prob_valid = model.predict(X_valid)\n        auc = roc_auc_score(y_valid, prob_valid, sample_weight = y_valid.apply(lambda x: weights_01[0] if x == 0 else weights_01[1]))\n        prgauc = prgauc_metric(y_valid, prob_valid, sample_weight = y_valid.apply(lambda x: weights_01[0] if x == 0 else weights_01[1]))\n        amex = amex_metric(y_true = pd.DataFrame({'target':y_valid}), y_pred = pd.DataFrame({'prediction':prob_valid}))\n        auc_scores.append(auc)\n        prgauc_scores.append(prgauc)\n        amex_scores.append(amex)\n\n    # Calculate the mean AUC score across folds\n    mean_auc = np.mean(auc_scores)\n    mean_prgauc = np.mean(prgauc_scores)\n    mean_amex = np.mean(amex_scores)\n    print(f\"Mean AUC across {n_splits} folds: {mean_auc:.4f}\")\n    print(f\"Mean PRGAUC across {n_splits} folds: {mean_prgauc:.4f}\")\n    print(f\"Mean AMEX across {n_splits} folds: {mean_amex:.4f}\")\n    print(\"Runtime:\", datetime.datetime.now() - start)\n    return {'auc':mean_auc, 'prgauc':mean_prgauc, 'amex':mean_amex}\n\n# Define hyperparameters for LightGBM\nparams = {\n    'objective': 'binary',\n    'metric': 'binary_logloss',\n    'boosting_type': 'gbdt',\n    'device': 'gpu', 'gpu_platform_id': 0, 'gpu_device_id': 0\n}\n\nbaseline_lgbm_metrics = lgbm_train(\n    df = df_train_lag_pq.drop(['customer_ID'], axis = 1),#.sample(50000, random_state = 42),\n    params = params,\n    weights_01 = [20,1]\n)","metadata":{"execution":{"iopub.status.busy":"2023-11-02T19:30:54.912994Z","iopub.status.idle":"2023-11-02T19:30:54.913912Z","shell.execute_reply.started":"2023-11-02T19:30:54.913607Z","shell.execute_reply":"2023-11-02T19:30:54.913639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### B.4 Baseline results:\n\n* AUC (CV folds): **0.9598**  \n* PRGAUC (CV folds): **0.3398**  \n* AMEX (CV folds): **-0.0016**  \n\nRuntime: 0:31:52.797085 (_no GPU_) / 0:18:12.447228 (_with P100 GPU_)  \n\n(training with GPU is **43% faster**!)\n\n**Comment:** `AUC` doesn´t seem to be useful at all, given the extreme unbalance of the target. `PRGAUC` is in line with the low `amex_metric`. The `amex_metric` is far from the leaderboard values (around ~0.8).","metadata":{}}]}