{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Overview\n\nThanks for reading. I would like to share time series RNN approach.<br>\nThis notebook is almost implementaion of [Recurrent Neural Networks For Accurate RSSI Indoor Localization](https://arxiv.org/pdf/1903.11703.pdf).\n![image.png](attachment:image.png)\n### dataset\nMy [time sequence unified  wifi dataset](https://www.kaggle.com/ebinan92/time-sequence-unified-wifi/metadata) is made from two exellent notebooks.\n* [LSTM by Keras with Unified Wi-Fi Feats](https://www.kaggle.com/kokitanisaka/lstm-by-keras-with-unified-wi-fi-feats) by [@Kouki](https://www.kaggle.com/kokitanisaka)<br>\n* [Indoor GBM+postprocessing XY prediction](https://www.kaggle.com/oxzplvifi/indoor-gbm-postprocessing-xy-prediction) by [@Oscar Villarreal Escamilla](https://www.kaggle.com/oxzplvifi)<br>\n\nThe main difference between my dataset and [@kouki's dataset](https://www.kaggle.com/kokitanisaka/indoorunifiedwifids) is t1_wifi column.\nCan be treated as time series data by this column. I'll upload dataset notebook if really needed.\n\n### model\n* One hot encoding for bssi features of all test's building is intractable. So I used entity embedding approch followed by @kouki.\n* Many models are possible, but here for simplicity. I used many to many rnn model.\n","attachments":{"image.png":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAA0AAAAFGCAYAAABDpf92AAAgAElEQVR4Aey9h58WRbY3vn/Auze997f33fted3+797rBXdc1giIqJsziImJEBQMKiigoYgLFhKCIiiggWTKSo2RmGILkDENOAsMwpCHpeT/fejg956np7qefME/oPv359FQ/09XVVd86dep861RV/4L0UAQUAUVAEVAEFAFFQBFQBBQBRSAiCPwiIuXUYioCioAioAgoAoqAIqAIKAKKgCJASoBUCBQBRUARUAQUAUVAEVAEFAFFIDIIKAGKTFVrQRUBRUARUAQUAUVAEVAEFAFFQAmQyoAioAgoAoqAIqAIKAKKgCKgCEQGASVAkalqLagioAgoAoqAIqAIKAKKgCKgCCgBUhlQBBQBRUARUAQUAUVAEVAEFIHIIKAEKDJVrQVVBBQBRUARUAQUAUVAEVAEFAElQCoDioAioAgoAoqAIqAIKAKKgCIQGQSUAEWmqlMv6KlTp+jw4cP0888/p55Ilp88c+aMyfNPP/2U5Tfr6xQBRSAZBI4ePUonTpxI5pGcxy1EnZhz0DQDioAioAjkEQJKgPKoMvIxK5s2baJLL72UfvGLX9B7772Xj1mslqdDhw7R7bffbvL87LPPUmVlZbU4+g9FQBHIPQJDhw417fSXv/wlLViwIPcZCpCDQtSJAYqlURQBRUARiBQCSoCyVN3woMydO5dmzZplTnSiqXpUdu3a5aSDNA8ePFhjpVi/fj1ddNFFxkh5++23037P8OHD6d/+7d/oj3/8I82fPz/t9NwSAB7169c3eX766adzMroMz9OKFSuceuJ65xBlRz6DygDibd26lbp37073338/3XzzzSbE77Vr11IiTxdG2KdOnUpt27Y1z+J5XE+YMIEqKircYDT/W7lyJX300UfUr18/OnLkiGc8txt458CBA83zixYtcoui/8sQAvv27YuTNeiIVA7IGXSTlNNjx46lklSgZwYNGmTaKQZY0tUHZ8+epa5duxLIVJ06dUw5AmUiyUiZ1olJvt5E17aVGLVsyUPinORfDAwSQve/9tprZrAQ/QEGDfE7UZ+Qf6XRHCkCqSGgBCg13JJ+6uOPP3Y6enT26KRhkCZ7wAhG5440+KxJz0wmO3t02iAkNZ3vfCBAcpSYy+sW1q1bl6ZPn+5LYMrKygieLLfn+X933HEHlZaWVhMnGLSTJ082hJPj2iEIac+ePQkGg32A9CI+SDBkIZkjk7KTzHujFhfTsVq3bh0nH3/+858NAU8Wi8WLF9Nvf/vbuLTgpampI5MESLZ7yGxN5Tsf5Dof8lBTMpGpdLMlD5nKbzbS2b9/P7Vv394MQtr9gPyNPuGNN95IetArG2XQdygCmUJACVCmkEyQDhuSUsk8//zzdPLkyQRPxt/u3bu3MU5+9atfGRKF9DLhmYl/S9WvIB0tPB0YPcJZXFxc9bB1BWP8nXfeMfkHAcRIU00csuPLlQdI4gZvF0bYcN54443VDExg8dlnnxHWLdkHvDPNmjVzDFKkBc/NBx98YMKLL77YuQdj0j6mTZvmdHZ4T4MGDejdd981Z6NGjZx78JgBN/tguVUCZCPj/3vKlCmmvu+9994a80RwDuyBBdYx8IagzQU9QIBh9OB5qV/c5CpomoniBSFAQfWLxAEkbtmyZYlen9J92bZrUvf6ZS4f8uCXv5q8F7RtZUsearKsmUx79erVcYOn6A9uuOEG4/VBf/LYY4+R7E+gB9Cnb9u2LZPZqJZW0PZd7cE8+Ech5z0P4Mt5FpQAZakK2JC88MIL6b777jNGBkZpoZSCHnBb33PPPebZl156yRjT+UCAMHWFja5ExhKMfJSjJjdVyDcC5GYkgdjAK4iRNmAHgw2j7/YhDUTUuT0NDVPf8Bw8QPaI948//kh33XWXI2vz5s2rZhBzPu6++25TL/b7WW6VANnI+P+W9Zbu1C7/N5GZ4sme1Vq1ahmSC5kC2d6zZ0+ix537W7ZsITyPZ9966y1n6muiNu0kkMJFEJyS0S/whqH9Y2OFmjrygXzkQx5qCt9E6QaRGU4jG/LA78rncOPGjU7bRvvGjAI33cD9CQbrEA8nBsrQZ9fUkUz7rqk8pJpuIec91TKH6TklQFmqTWlIYh0MyA+USzKjtN9//73x+uBZzNHndS5uBnamihWko803JVAIBAj1g9H5bt26OR0NvGNyxB6bN/DUNximMFC9DjxnrwPConKM8gWRM7fpb3iXlFvIQjJHENlJJr1CipuMkZZuueRIN3TCkCFDnHofNWpU4OTZuwziNHPmzIIkQIELm0bEfJDrfMhDGhCm9Wg221ZaGc2ThzHLBINnTGgwZT7RzBOsKcSgGD8D3VBTR77ZD8mUs5Dznkw5wxpXCVCWalYakpiawSO2QUdpobAwZQ4KCeHu3bsjRYBg3B8/fjxQbaVCgOCZyuRi76AGClzo559/vqlXe7qeLIfXFDU/QKShgOtUDim3KFMyR1AMvNJEfdukziuu3/9Rt9neCVBij06yJg+bAAH3W265xcgUprUEkWvpXcagjKy7VGUnSJmD4JQNIyMZ/SKxCTr4lCv9EqQOEAf5g1c+E+0t6DtlPMgo8hDkCCIzQdLxi5MpPDAwBT0mB7b83mvfQ9uGFyudY+HChfTrX//a6INkvDlz5sxxnsPMk5ryAmWqfSfThtPBUz6bqbzLNPU6ewgoAcoS1rYhiZFZHl0JMkqLqXLsNUJ8aRx7dcI7duyghx56yKxFwLxprwMLI5966ikTr1evXnHRvDp7ELIOHTqYZ2rXru2UBXOIeb0LQqxtkd4Fnr+NtSjYvUwefK9FixZmZzIoNIxEY3oXY4UpY61atTI7osln5bXExiYVMh7K/f7778dtEMDrZPDedIwBL9zk+3Et49l5hUGCNSQo+3XXXWdIr/283++xY8c6uGGTg1QOW26TSUOWzU1GMf0OdQ05YflEnWAnIrkQH/U/e/ZsXyMCcot08B4YDJDPkSNHEjaZYNlBmlj/hHe4HZwGy59bHPwP6eJdaFtoY3zAQHjmmWfMPTmXHu0D8fnE85k8bAJ04MAB41lGuWH4wABKdLB3mePLuvMiQNxebRzsd2FdIK8RBOGXh5cxm6p+4Tps2rRptfbC+eX6TVW/SGzc5JrLlw/6hfPiFm7fvt2sycS6Qm4jCCG70ItuU6RkOsAPdfv44487U3nxPLzV0PvQw14Hpud+9dVXhqizlxrPoo1izePy5cvj2nuqbctPHuy8pYoH97OQcZZvbFb05JNPOrigjE2aNCHsqpnowLPY1ETqQNTRiy++6NvvuaUL4sXrboFvEFuD0wEpxQAKnkP+7W3qpf4Gzl4H9BM2XrD1czr2Q7ptOBd598JH/587BJQAZQl725BE5wLvD5RLos0QoMQwKou47DGSRr5XJyw7ai8jBsWX8ey0vO6hQ4LCR578Ttuo9zJ4kA++B28HOoF27dp5po0OASNUbofExn4/4gPPRDujoUxQ2kFGz93y4IWbHVd6gNDpydE+EMdXX33VwQALVYOOkuI9Mm0YJTbhtPPi9tuWW7c4Xv9LhIGsJ2yXjWmdtjHGsoUOePDgwXFGkXwv5xP1jd3wHnnkEQc3ToND7KKIvNkHp5HI28bx7HVRQXf+w/OZPGwCBFzlgEmiabbSu8weI1l3XrqD26uNg102jgf8bW+Y171U9YtX3SBP/K509YvExq0u80m/2HWB39Ah2Nae1x9yu7BDtEXoSZTHPiBzMKwlebGfd2tnSAtTNBO9G/cRj9+datvykwcuU7p4SHnAZymgp7zKB1KDwRyvA/ck8bExxYAYBsaCHnIQ7aqrrkp6QwMMnHEe7EE0qb/d2gHnUcbj/jjV9p2pNizzlK28Mx4a5g8CSoCyVBe2IoZiZ1KTaDMESZbYmAnSgKVihuLwOmQ8Wxl43YNxjql8MFp79OjhKEksnuZviCC0jW5WYFCqXsYQDCr2+lx++eU0YsQI2rlzpxnJxgg7K+R69eqR2/dOJDascGXZZSeD92CEHMQDI5ownuWWwp07d47zYMl0/K69cJPPoONNNDoHjGSH2KZNG08Phkwb1zBSWrZs6eAFXLErXDKeLVtu7Xf4/U6EgawnbAyCcsJw6NixI61atcrITqdOnRwjy8/Q5nxCJq699lpTZtQtDDhMF4UsYvoHyw5GYzEKKA9OI1UCBCKBUVK8C+2A34X2IdtEpndVciNAktTwoIksq7yWZIlHiGXdeekObst+9YL3cDy/Nm/fS1W/cB265YnzgXvp6BeJDd5nH/miX+x84Tf6HRjoTFxAcrC+A/oV7RHtDjsB8n20SVtPo26gF1m+gSU8iPB4YZ0iyBX0Nu7bssOyivQfeOAB8ykIPId3Y91Jnz59HH0HAoWBMBypti0/ecgUHlIeHn74YYMdcMUmN8BjyZIlcX0KNqaBB8w+ZD+PwUWkC13N/RL0otdmNXZa/FsSRwwK2ZvocDyvEIOMXM/QaUxIEV/qb7d2wGnKeNwfp9q+M9WGZZ6ylXfGQ8P8QUAJUJbqwk0Ry7m5TGzcsgOjBEpIEqUgDVgqZrsjku+R8Wxl4HeP00AHyUrS7z2IzwoM8e2OVd7D/eeee66awkZHKL0iboszJTascDmvGHliQxidCTpd+wAxgfcHeUi0+YD9LP/2ww2dCIjbK6+84hga2OzAzduEvMDYYIMEeYJRgqlciaaoIC/Ih/3dKDZYkHaiw01uEz3D9/0wQBxZTygXjAYQBdnJ4vrrr7925MseheR3cT6RDrD68MMPqy30LS8vdzxDiAOjTR6cRqoESKYlZdmWcxkvE9dsVKLsMu+sN/B/Jjb2+4AvD8RIoiTrzqtNcxndyIZ8D8dDPmws/O5xGsnoF65DtzzJdyEvqeoXiQ3eJ4980C8yP/Y1CAXrAy+vMGRCkiT2CnJakjDjo8xuU0phaMODb+9MCcN39OjR1aa4cdp2e8c0XvuQ9WjLkx3XTx4QNxN4SHmAXAFfYCQP2acgDgai7AP/wz2ehmrfx2+sZQRGQQ+QL17/Y88wCJKGnEVg96VSf9vtQKYt49lpIF4y7VvWfTptWOYpW3mXmOh1fiCgBChL9eCmiOUcW2l8yCzJOHKqXJAGLBUzFIfXIePZysDvHqeXqgKzOy+p3LwIAd4pO2C3US2Jja1wea0DDGC/D9HKjsMPO8bADiVuMOx5/QfWpEgyg2sQIXzs1OtA5wmDRHqCoPzhLQFRczNAZFqYo47RVjwjzyAfYXWTW5m237XEwJYrPCfrCWXzmhoit2f26sQ5nyif31RBrn/EA5GGQcYHpyFJBN+TIcdzM7I5npRlW845TqZCLwIkR5RtI5bfLePIQRhZd17yz2X0wwHv4XjA3MbC7x7nMRn94lc38l3p6BeJDd4nD5avbOoXOw8yP/Y1T2lC/saMGWPfdn5LDyIMaOhDPvij3n6GOsdNJVy3bh3hcxGQF7zLPmQ92vJkx/WTB8TNBB5SHrxIJd4lBzz9yoWyA4NMHLLtJCMn/G5ZNrsvlfrbL20Zz04D75F5RN36HbLu02nDMk/ZyrtfufRebhBQApQl3L0UcaJRWlaa6LDkiHWQBiyVl59ikfFsZeB3j6FLVYHZnZdUbvY9fhdCjC7yGg+3ec0SG6lwMXLGU84SbSoAUgFSik7YxkTmxeta4oY03E4sNl60aFHgET3kCYTHnl8OgoVpgn5T23Bv+vTpcZsCcJ6wYHzv3r2uRfGSW9fI1j8lBm4YetWTlUxcfcv6lPGC5hPT4VD3KLs9n57TCAsBgryzd8fLWGX9I73LwFXWnZfu4PZaiAQoHf0isZFynSv9IvMg24R9jd3InnjiCSP7ifQfnpWbqPTv398kJ3UvpmnB45Xpwwtffg/LHdqwXz0iPrdpNxnNBB54R6L8cr6l7nGrM/nZgiBbVXO6fiHwYT3v9k6/Z3FPls3WvVJ/+6Ut49lp4B0yj166hvMZtO6lnCayEbKVdy6DhvmDgBKgLNUFGhkUka2I5ei29PAgWxid5i+z29tQSqXi1YCl8vJTLDKenZbfPYYuUwosqHKDocFrLNxGyyQ2UuHKkXKM1I0bNy5ubYZcp4F7iIM6k2lwmROFEjfsBsRpw5PDU/CQNjxD+EhdMgfmj3fp0iXOIwSCjLSBjd8BbxIW6vIHdblz9JoO6CW3fu/gexIDW64Qx6ue+HkOZb151UXQfEqjBwRXes84jbAQIODHAyioZ+nhwT0v7zLuybrz0h3cXm2dxvXGIcdDHmyD1e8eP5+MfuE6dMtTkHfhnYn0i8RGyrWU02zqF5kHxswtlIM6Xu1IPienP/E7pBGPEfhMbC+PwRkQKegDnBgUYg8Qv1fmK2g94hk8D7lzk4dM4IF3eMmDzDOupb5zKxfWJGJtIutkTFWGrg4yVdl+F/+WbcftnRzPK5Rls2UmUXk4TRnPTgNxZB5Rt35H0LpP1IZlnvxwkfHSzbtfufRebhBQApQl3L0UsSQ59iisJEf2WhfZML0asFRefopFxrPT8rvH0GVKgQVVbngv4+lmVElspNJCJxtk5zrugDiUaXCZE4V+uKHDx65nIC14BwiS2/qfRO9ABw4jhPPp1sl7pYE8TJgwIW7XNYw6yilheJZxTiZtfqcfBojjVU/8PIfSsPSqi6D5lGnZZeI0wkSAJMmxp9kyObK9y8Bd1p2X7uD2auPI9cYhx4OcQlfIw+8ex0tGv3AduuUpyLv4nZyOW54lNojHR77oF86PHXrl247Hv2V8Jjvyf7Ls/EzQEMY++jR8r4r1IOsxGbq9I5V6dJOHZMsi4zMeKK/8v1t+GROp77ziYW2oHCADFiCDibYV53fYoZzK7aU77Wfkb9n2MrkJgtc7ULd+Ryp179aGg9QF8iHjueEn8UmUd79y6b3cIKAEKEu4Q+GhIbopYjZEcF+O0qKD8HpGNkwvZSoVs1/jlPHstPzuMXTJKAE/BeZ3j9/FIeOZSLlJpSUxA9mEQsd6kUSnnP/O708UJsJNLop1M0ATpc/3YeBKEuS1SQDHt0O5Y5XbtBjG2U1u7bTs34kwkPUh68lOR5IWr3hB8ynTssvEaYSJAAFLnuaGtsKbIciBFxih+HaQPGTdeekObq82jjIdXHM8t7bqd4/TSUa/cB265SnIu/idnI5bniU2iMeHlOdc6xfOkwy98i3jyGsZn9feyf/JssvnEl1jqhfvEgd8MaUXuze6rZN0e0cq9egmD8mWRcZnPFBW+X+3/DIeUj784mH9FXav5J0KgRFOTHX22pac32GHche4VKYsYhMLfj9Pg+R3BC2PjOemv5Np36nUPfKPd8hD5smvLmS8dPMu36/X+YGAEqAs1QMaGRqimyLGyCFPSeJRWvk/e2ocsiwbplcDlooZisPrkPHstPzucXqZUmBBlZt0b7vhKbGRSktOf0pk5HLZUg2D4IYpJjCUIBdei9SDvF9uVWrXX6LnJSZuWPrJbaK0E2HgVU92upK0yPqU8YLmU5bXlgFOw/6/fA+uOZ4bXhw3qCxz/HRCiY9b3qUnmeVM/k8OunA+ZN156Q4uox8OSI/juRkifvc4L8noF7+6CfIuvDORfpHY4H18+MkWx8lU6JUHv/SlMSw9GF7PyM0IuJzyvfYmIl7pyP/DwwGyw30hNqKxp3jJd/B7ZRpB6xHP+MlDJvDAOxLll/Mu9Z1buTgeh5BD7FKHTx+wlwz9BfqNoIf8DhCetXen80sHgyS846q9EQaeC1oeGc9NfyfTvoPWfaI2LPPkVxcyXrp598Na7+UGASVAWcLdTxEjC+ztQceAUVq5m5Dc/ICzKxumVwOWitlt1xlOS8az0/K7x89nSoEFVW5ygaOb10JiI5UWvvWDkTtg7LYwksuTiTAIbnKnpfPPPz+pjk3mUeJv15+M53YtjWc3QzaR3Lqlyf9LhIFXPfHzHMo8yvrk+wiD5lOuYbB3EOQ0EskGx3PDi/MUVJY5fjqhxMeNAElvD2+GwPrGyyiSdYeyuB1cRjfjSMbneIVCgBLpF4kNZIGPfNMvnC8O5RQ9ewMQjiND3pYZ9cYyINus3X7ks17X7FGAQe+1C6cXvpymnzxxHA792mom8MB7EuWX8yKxk3LD971CGPPYOpxJEKYqBz3wLG/8g3rEJwXwvyCHLJe9BhnPy/JIj5idtoznpr9l/8VyZqfBv4PWfaI2LPOUrbxzGTTMHwSUAGWpLvwUMbIgR2SbNWtGOKGw3BQP4ssG7KVM5QiXVyOHMpTffLDTkkrQvsfQZUqBSeXmt5Wx3AbbTaFKbOz7cttTN2LJZUo3DIIb3iGnJ7mNxAfJh0wj2Slw0tPoRiYTya1f/hJh4FdPMl1p4Nv1yfE4nzAS/DaD4IEFtC0YBtIYgGGB//sZ9PhuFDaMQLxCIUDASE6zbdu2LTVo0MCUwc27jPiy7ryMEjZmgYXb91qQjj1FE7pCHrLN2/c4XjL6heXArW7ku9LRLxIbWyfmm35hDBFKgsZEWN6X13KKriTJMg35f/ms37Vf/fBzfvgijqxHL5nhtPzeJ8uSKh54T6L8cl6kvrPlhuN4hfJZfL8Kg2dBD9n2sTlHkE13kP5LL71kdATaNwZM7CMogZTTrN30dzLtW9Z9Om04F3m38dPfuUdACVCW6sBPESML0hsAhcOnm+JBfKkQvZSpdH+7GQQw/iZNmhS3m5idVhDlLncLSmTESwVmd17yHgxZlN2eHiE7ZmDEaxpkNUpsbIUryRMMWbcPocq0Ur0OghvSlt9h4emP/M4ZM2YYIsxfQ+f/y1Aa5NIggeGJjvLbb7/17Cxt8osdBwt9EwTIBL4pBLmW5AaYyV2W3Aweue2v24YQ+FYTNqzgtunWprhu5Oi5m4xyPISQ6SFDhphpkDCgQfiSOSRBdPMAIS1JdDn/aGNegwBSftEu3Y5EC6yRrw4dOjh44b1+bd6+x+9MRr/46dlM6ReJja0v802/MIYcygEATIODnnA7pNFqk2Q54IKBNTd5RX/2+eefU9++feOST/QNIbQFLPhnb4eNLxJLpm35yQPSygQefvIgCy/7JbdyoS7s/o6fl57rZDxAeN4mM+j3MBXR67DbLQZhseGOfdgEElOx7QM6gndTRfu3+2PET6Z9Z6oN5yLvNjb6O/cIKAHKUh0kUsTIhlTGUBZQHPAMuR2JlCmegQHYrVs3xwCBIoOywQ5g2Eq5c+fOZgEqjEGceKetmIMod6mckWcscoUixwc4S0pK4gxRqcBsg0feQ15wooMtLS01hjm+VfPmm286naOXYpbY2AoXBj7KzeljoSm2XZUdD5TjqlWrDBb4WCg8ackeQXDjOgJp5PxIY5nxgDHQvHlzs5U2OiKUD8QH0yLkYmI5IiYxQJwePXoYWcL/caJ8cm456m3t2rXVihlEbqs9dO4fiTCQebTrSaYpDXyveJxPxhGYAQ/IDOocMvTUU085OLds2bKa4Sa/Co+F2V9++aUhTZANbEfLi5IxXRHv8SNAslPH4mOMup4+fZo2bNhAK1eulMWLG8zwSzPuIfFD4uNFgBCdp70xRl7eZcSVdQc5dDtAqnjHKuCNtgm9Av0C+eItfX//+987bdavzdv3+J3J6BeWAzccuT1x+RGmol8kNra+zIV+kdvs83b7Mly2bJkzsAEjW5J4fCAZ8sj6D3X6zTffOLtDuukFWe+M4ebNm00akEUYwtxWbHxQx0xuMIVOfn8Mg0HQSbJ+7OchE8m0LT95QFqZwMNPHliGEUp9Z5eL2zDaL77Xhj6ID+Ddrl07gwuw8xq04PhuoRwoA74YJEIe0E65T9i5cycNGDAg7ltxbvUv08fHdLk+69SpY/ooyBIGm1iOoEvPO+88k383/Z1M+85UG0YZsp13iZte5wcCSoCyVA+JFDGyAUUHo4Q7ALcRec6unzLlOAhtxcdpcwhDAQoVign/sxVzEOVud/qcNkJb4UkFZhs88h5ICpSvTEte4x46drdDYmO/H/HR6bVq1cozbfmeSy+9tEYJEPIjpyjwInX8H0azHwYynygPysWHX53I53CN3YXcRu+QVhC55XfaYSLZSVRPnB4bB8irW33KfGLL2A8//DDOq2mXF4QEnb192B4x+zn8xodop0yZYmTHzcjmNJFnkCy3NOw2JnFIRd4kPn4ESE6zRb68vMsog6w7tEuvQ3oK3MoKIxfkEVjhvl+bt+/xO/1k2ZYHP3nNlH6R2Nh1iTxnW7+44S7/Z8sERv95GqeMZ1/76QUMliTSTTCy7Tq12wWMZwwy1a5d22kr8ExhwAf5ccPXTkPm247vJw8sX+nikUge+D2yndv5lB4JLg8wwcey+TdCOcjF6QYNt23bFmdfyHTdrkFisRGG3wFZl7uQ2ulABjCo17FjR1MOu70i7WTad6baMN6b7bz74aj3coOAEqAs4c7rCxIZOHD/Q4lAcSxevNgzdyBL/E2bRC5xjLJh1I9HapA+rh9//HHjpcHH7FiJ2WnJdUT2PZk5LDqEcWi/A+WRU5F43QDiwVMkD6nc0HHCg4Q82mnif37TwiQ2KJfbx/owSoVpUjfccENcB8MKHPhjdBgdvcy/zK/fdVDckAYUMYgP3m3LBzpNYOhlbKCTmjlzphl1t/ODkfji4mKDIUbhuGwcoozvvvuu6/QGTiuo3HJ8GSbCIEg9IT0pn24dKOLYhg7IsRxMQJmBAWTUbToH5xuYocPmDzEyVviNaWqQG/5iu11XnAaHaHdNmzaNwx15wEd25YH8YPoj3pXKwnKJj99WtzCyXn75ZfMejNb6tSFZd2izXgfaRlFRUdyoMcoB2YL3GcYqp+XW5v30gXxnUP3iJ6+Z0i9cHpTTSyfWtH6BzDRu3DhOtlhW7RDxMCIvD+gceIVtOee6e+2113zbCdKCfsbAi61b8Bv/95IvyARkAzIi84q88JRdbs9e+AZtW37ykCk8gsgD3iX1nU2AcB91Cn1s49GwdU0AACAASURBVAKMJDYy38leYzqc2xbbsh7QpyAO4gY50DYxWGvLAfQvD1JyPXj1x0Hbd6baMJcrm3nnd2qYPwgoAcqfuqjxnBw9etRMg4JRj+uaONC5wWjHietkDlu58bPIKxQpRrCTTZPT8AuRPjxlMGzhGcDaqVRIj987MnEP+UT+kE/kNxksYJChAwaOqH9cw9gPy8EGk/TKoA5RTsxDB25yWkmicgMvfBsHeAFr/E71QL1xm3DLA3bDAjnAiWkZhXhAloA1picli3Uy5VX9kgxaieNyvUG3QtZRh8nKOmQabQRTiRG6ybhbTuznkn0v0kzUttze6/e/TODhl36Qe5wH7vNqqj8C/tBLeA/OVPpsWR60TbR96IBU85yofdeUjZCNvEus9Do/EFAClB/1oLlIcncfBUwRkAi4ESB5P1+vpfcPuy4FHXXN1/Lkc768jKd8zrPmTRFQBKoQ0DZchYVepY+AEqD0MdQUMoSAKrcMARnBZAqVAPHGJ4mmpEWwSjNeZNUvGYdUE1QEsoqAtuGswh36lykBCn0VF04BVbkVTl3lW04LkQBh2gXWNBXy1Ld8kwO//Kh+8UNH7ykC+Y+AtuH8r6NCyqESoEKqrZDnVZVbyCu4BotXiAQIu0dhoTA+yAoypEfNIqD6pWbx1dQVgZpGQNtwTSMcrfSVAEWrvvO6tCNGjDCj4dgFB4sy9VAEgiKAra+xk9G1115L27dvD/qYxosQAqpfIlTZWtRQIqBtOJTVmrNCKQHKGfT6YkVAEVAEFAFFQBFQBBQBRUARyDYCSoCyjbi+TxFQBBQBRUARUAQUAUVAEVAEcoaAEqCcQa8vVgQUAUVAEVAEFAFFQBFQBBSBbCOgBCjbiOv7FAFFQBFQBBQBRUARUAQUAUUgZwgoAcoZ9PpiRUARUAQUAUVAEVAEFAFFQBHINgJKgLKNuL5PEVAEFAFFQBFQBBQBRUARUARyhoASoJxBry9WBBQBRUARUAQUAUVAEVAEFIFsI6AEKNuI6/sUAUVAEVAEFAFFQBFQBBQBRSBnCCgByhn0+mJFQBFQBBQBRUARUAQUAUVAEcg2AkqAso24vk8RUAQUAUVAEVAEFAFFQBFQBHKGgBKgnEGvL1YEFAFFQBFQBBQBRUARUAQUgWwjoAQo24jr+xQBRUARUAQUAUVAEVAEFAFFIGcIKAHKGfT6YkVAEVAEFAFFQBFQBBQBRUARyDYCSoCyjbi+TxFQBBQBRUARUAQUAUVAEVAEcoaAEqCcQa8vVgQUAUVAEVAEFAFFQBFQBBSBbCOQFgH6+eef6ezJk+Y8U1lJekYbA5YFyIUeioAioAgoAoqAIqAIKAKKQD4ikDIBgpF75swZWv3CyzT7TxfT3Auv0FMxoMW3/oMqf9xP6ZKgn376iSr37KUTu/foGXEMIAfpylM+Kl/NkyKgCCgCioAioAjkBoGUCBCMkbNnz9KJEydo6YOPU/m4iUTbthNt2xELt58L9Xfk8Ci+sh4d3b3HyEcqIs2ytbVnb5p7wWW04Oqb9Iw4BnMvuJS29+5HIMV6KAKKgCKgCCgCioAikC4CKROgU6dO0aFDh2hx4yZU/t14oq3biLZs0zDiOBTXrkcHt2yl06dPpzRqbzw/lZW0oXsP2vLmO+eINci1nlHFoPT1jlTas7fxOGfLEwQ5PDh3Ph2YM4/2z56rZ4QxgAxAFrIle+l26vq8IpALBNA+9FQMciF7qb4zZQJUWVlJ+/bto4WNHo4RoNItRFu2EiEsPRfq78jhAQK0d8NGOnnyZEoGA6ZVVlRU0IrOH9OWN96JkWoQa0mu9Xek8Ch9rSOt696DoHOy4QXCO06Ul1PR5dfQykee0DPiGCy79yEque6WrBLwVDt0fU4RyCYCIDzQl3vGjKfd343TUzEwsnBwXlFK9l82ZRfvSpkAYfrb7t27qeTeh6h81FiizaVEm8SpvyOJR1Gt62j3uvUpESAoU3iOysrKaHGnD2jL62+fI5Ag1XpGFYPS9h1o1cfd6fjx4zVOgFgGy/fto3l/v1I9j+p5pcqSRVRS79asEfBsGwH6PkUgVQSwFOL4/gM076+X07oWrfVUDGjFA4/Rwvp3FcSAUVoEaNeuXbSg4YNUPmoM0cZNRBs3x8JN50L9HTk8impdS7vWrkuZAGFq5YEDB2jh2+/Rltc7xkikIdObqwil/iZCG+MBh5DjUfrqW7SiSzc6duxYjRMg4/3B4M6WLTECBC+2mda69Zw3Un9HDY/KBQtpwXW3GM80PNR6KAKKAJkRfvTXZTt2UtFldWODRVvPrQU34XYi/R1bGx8hPComTaMfGj1cEANG6ROgfzxI5SO+I1q/MXZuOBfq70jiAQK0c83atAjQ/v37qaTju7TltQ7niDXItZ5RxaD01Tdp+UefZI0AgWhtXb8+RoDgedwsTv0dOTwqi0oMAcKaV3io9VAEFAEyg1GYCbRrw0aaf+nVsSUQZno6Bov4xPR1vhaDSM7/9H7Y8KmYMJkWN3yQjh49mvJmWNlqXxkgQA9Q+fDRROvWE63fEAvXnQv1d+TwyCgBat8hJlOGTG8gMuR6wzmyrb+jgkdpu+wSICjuzWvW0NyLagtP2+Zz03zZ86a/Y9Oew49H5fwFVHxNfTp48KASoGxZJvqevEcA09+gK7esWh0jQFj7jcEiHiTS35HEAwRo4T33U3l5uZkGl8+CnBkCNGwk0Zp1eioGlFEC9OpbMQK5dr2GGGCIKA6lr7yRVQ8QOvVNq1fHCNCGTUR6RhqDyrlFVHzNzUqA8tmS0bxlHQEmQJuWr4gRIDMtGwMiekYZg4pxEw0Bgsc836cMZ44ArV5DxOeatVXX+J/+jgweRVdkcArcq2/GZGctyPXacwR7LZH+jhQepa+8niMCVCvmgXQ8j/A+Ss+j/o4CHpVz5ysByrp5rS/MdwSqESAeKMJ0db5GqL8jhUfF2Am0sEFj85mcaBCgb4cRrVylp2JAGSVA7d6IEcdVq4n4BKHma4T6O/R4lLZ9LTcE6G+1Yl43eN70jCwGlbPmKgHKd2tc85d1BKoRIAwO6Rl5DCrGjI8YARo8jGjFyhgBQsgnCAFfF/D9M0uXU3nRAirDOb+YTv2wNNTldeoshfrLKAF65Y14Ug3CI0lmAf/+acVKKl+wMCZTRQuocskPMSITkvKZespQ/eSUABnPI7yPekYVg8qZc5QAZd281hfmOwLxBKjOuWnq684NlKzT32baevTwqBg9LmIEaNAQomXLQ3eWTppG9a7/B/3rby91zj/9/XpaNnJs6MqaqfrLKAF6+XWi5StCdZbPL6bGDR935Illa9LX/UNVzkzWW2mbV3PnAQKJ0zPSGFTOmKUEKN+tcc1f1hGoRoB0kEgHytaspYpRYyJIgJYuI+ITZIivERbw722Tp9NddzzkGKwN7nqEDsyZF5rymXrKYP0UXXFN5rbBfvm1mOyABEmCXeC/y+YVUfNHWzoyVeeau2nj+MlVZSzw8hnyk8H6yikBkp5r9oiGJDyzzPJsQ1eHuLzGK5lC+Sq/n6kEKOvmtb4w3xGII0CX1KkaJNFp6VVYYPAsYnhUjPouQgTongeofMBgIkzhCeFZMXc+PXr/U46x2r7lK3SqZGEoy5qJ+ssoAWrbngjTDUN2nlq4iNq3bOfIFOSrYl5R6MqZqXorfbFd7jxAIfNAgpyWTp5O9a5v6MgfvJDGsz1qrHohXeq7ctr3SoDy3RrX/GUdgWoESA4urDi3/IEHiwr8d/VlENZgUYGXjzKY/4rho6JEgO6n8v6DiBYvIVr8w7kQ1+H4vWPSVKp/832OsdD3vY9FGcNX3li9pV5/RZdn0APUpn0MaxBryFNIwhipftKRqbdbv05nFy4KTfkyXU85JUDSg82e0hCEvp7tEJQvzrOdZnkqp05XApR181pfmO8IuBIgOYAAQlTgv30Hi0JQvpqon4phIyNGgPoNJFq0OJRnyYChdN6f6hhjFSF+h7WsmShXRgnQS68SLVoSunPTmIlU99q7HQI0vOsXoStjJuuttPUrufMAhcz7yF45eBzjPdvtCJ5Jvq9hlee5cvJUKq6r3wHKd4Nc85ddBKoRILnsIUTXroNFs+fGL4MIUXnjlq+kUK6KoSNo4d1R2QYbX3ztO4BoQQkRpobxGZLffTt1cQzV+jc1oh1Yq4EyhqR8ma6vzBKgdkQwykCuEfJZ4L9n9R7gyJSZejR0VFXZQljedOuvtPXLuSNAIfI8Ss+cq2c7JB5WWc5MeI5jBOgm/RBqdu1rfVueI1CNAPFgEYxmvkZY4L8r5hfHDxY9Zw0WFXj5Ml0/0SNAffoRFS8I3Xlq7nxq/0xbx1ht/vAzdHTGrNCVM5N1l1EC9OIrMaIJshmis8cb7zkydeet99PeSVNDVb5M11XpC21zRICuCK23Vz3bwWcsVE6cTMV1lQDluT2u2csyAtUIkFn2wNPnwxO6DhaFtKzxSyBSq8OKwUMj5gHq3ZeoqDh2zi+qusb/Cvj3gcnTqMHtDzrG6scvdyjo8pg6quH6KLq8buZ2gWv9SkyWQK4hSyEIT8ycTc8/9pwjU22eeIFOzpkbmvLVRD2VtmqTOwLEXseQeeb6durqyGD9mxvRjglTQudpTdfzyM9XTpikBCjLxrW+Lv8RiCdAV8XW6cqlECAJIfjtOVgUkvI5dZSh8sQI0H106NAhOnPmTF4L8i9Syd3PP/9MJ06coF27dtECTIHr9Q3RvPmhO5cNGEJ/uqieYyhM+rRn6MqY6XrLKAF64eUY4QRpC8m5Y8x4qn/TvY5M9XitU2jKVlN1VPp8jgjQhVeE0jN3al5Rdc/2zNmhLGsmvJGV4yYoAUrFUNBnQo1ANQIkB4tCdO06WBSi8jlLCzJUpopBQ2jh3VEiQF/3IcIo9tx5oQpHd+7uGKp1rr6DNg4dGary1UR9ZZYAtY3JFMg1ZIvPAv69+JuBdN4fr3LkqujrfjFSzWVDWMDlM3WU4fyXPv9SbjxAIEDS88hebvZEFujvA1Omx3u2X+kYK2eBlsfxDNdQ/ivHjlcCFGpTXguXCgJxBOjiq2IDKNCNmK4ektB1GQQGi0JSvpooR8WAwREiQA3up/KevYlmzwnVeXbmLHq7xcuOodr4rkeofOKUUJWxJuosowSoVZsY4QS5DsnZt8MHjkzVu7YBlQ4fHZqy1VQdlT73Yu4IUEg8j9I7t2zg0HjPdvee6oX0qefK78YpAUrFQtZnQo1ANQLEAxAhCl2XQYSofGaAL8Plqeg/KEoEqDGV9+xFNHMW0azZoQkrJk6mR+993DFW3362LZ2dMTM05aup+iq6LINrgJ5vE5MpyBVOkGy+LsDfp76fSe2fau3IVPMHn6ajk6dWlanAy1dT9VPasnXuCFAIPXOunm31PHp6XitHj1ECFGpTXguXCgLVCBAGEeD9D1HougwiROWrifqq6DcwYgSox1dEIAchOjcNHkZ1r77TMVaHv9vFt3wcv6TnN77xwoSRW1kySoCeeylGOEGuQ3AeHDeB7r3jIUemPmr9Gv3MpDoE5aupOiptkUMCBM8jiGlIwrOzZsd7tu8+59kOSflqop4qR45WApSKhazPhBqBagRIDhaF5NpzsCgk5XOWFWSwPBXf9I8YAfqiJ9H074m+nxGacNanXzqG6p/+Vo+W9e7vWr6y78ZRw9sfcOKWfNnbNV7Y8PEqT2YJ0IsxLEESIFt8Fujv1f0G0WVX1HdkZVKXzzzLV9KzjxPvX397qbnu3OpVz/gONsCoQPFxymDlv7TFC7nzABWgp9HPE1cxaQo9em9TR7YwzRfTfQvZs+pX3kx4jg0BuvpG/Q5QqM15LVyyCFQjQCEbRNHBotQG/yq+6UcL72oUkV3gGjSm8s+/JJo2nWjqtFiI6wL/3aPtW46RcOfNjWjv8FFx5aucOIlaP/IsNbztASobPZaGv9PZxC/p0bsKgxDhEbQ+iy67OnPbYLdoXYUly1QBh5M++tSRqcsur0+rvxngWr7Oz7cz8SBTBvdp04n/h5D/F5Ww9JlWuSNAIAcgZHwW+O9N3w63PNtdq8qGMp4rX9nYCdRQfAIAJBwecTzvYCHiO/8rcHy4/LI8lcNGULESoGTtY40fcgSqESA5PT0E17HBIrEMAoNF0HkhKFtNlqGid9+IEaDuXxBNnkI0ZWoowpPjJ1KbJs86xurzDzWnE2PH+5ZveMcPYwTo869844UJJ7f6ziwBeiGGJeQKJwglXxfg749faO/IVINbGtOBkd+5lqfzc69Q60eeocoJE537ZaO+M2S74W33E67DgEfQ+ix95vncESD2OobEszare09HBo1nu885z7YoX9mY8carDQKEa3jmYv97kOpefQdtGjQ0Up7IyqHDlQCF3JjX4iWPQDwButIZPJGDB851AQ6WBB0swvIInqXhF3Z+oX0kBo8qvu4TMQLU7TOiSZOJJk4KRVg5Ziy1fri5I9RPNWpK5SNH05Ivvqa2TZ6l/UOGVyvn8HO7e5V81jM0OKRSnxklQM+2iuEM2QrB2bll1a6Ct93QkHYNHkqb+w0ysrap70DfMpaNGEUNb73fxK0cO843bhiwkmUoffq5HBGgy0PnbXP1bFteVXi04dmO82ZPm+7u5baeDaNXsnLIMCVAydvH+kTIEahGgORgUQiufQeLRPkwNb11kxZUOWmKGRga3ilGiDCNHYNHGDDCwBH+70zzFs+H7X8VX/eOEgG6j8o/6U40fgIRRqxDEP48fgJ91KKtQ4CY1eP7Lfh/5XdjqpVz+Jvvmfgln/YIDQ6p1GfRpRmcAvfM8zEsQawhW3wW6O/hb71fTaYgW20eeYYODRvpW76S7rE1aQgNMWUsEBYoHkHrs/TplrkhQH+9PFSeNuPZfrSFI4PGsz1uQiDPauX4CcYraTyQo8c4nskoeCIrvx2qBCjkxrwWL3kE4gjQ36+MHywKwfR/38Gic+UzSyGaPEslX/Qy5effvDQCA0IYUGp2b1PaNODbKoxCgE/cYJcoT0XPXhEjQF27EWFUGue48VXXBfz7wOCh1OqBJ+m8P1xFv//z1dSycTNa1aMX/czl5PBceYe/8W6MAHX7IhTlT7U+M0qAmj8XkydgHILz+MhR9EHzl4w8Qa4evvMhmtu1O/3E7cSjjGVDhlHDWxrT8DffDQUOydZl6VMtckeAQuB5ZG8aPIfunu1eMc/20BGunkX2PsIDiWtOLyph5cDBSoCSt4/1iZAjUI0AyenpBX7tOViUoFyxqer3x6awwxmQIH4Y71f0+CpCBOju+6i8yydE8IqMGRvZcPjr78QI0MefRRqHokvrZG4ThKdbxrAEQYBsRTCsHDGKWj/4FA1/vVMky496L33y2dwRIPZos8etgH/7erZbvhzv2Z4wkdjryB7wzi3aVnm+Q4BHnIfbpzyV/QcqARLG/JkzZ6n8cAUdO3Zc/DfYZWXlSSo7dJgQJnvgfXgv3q9H7hGoRoAivgwCA0Il59ZYYraHGSAKybKQZJZDVHzRkxbeGZVd4ECAOnclGv0d0ajRVWfEfg9/7RwB6tI9hkHEys/1n1EC9GSLmFwBywie2H2q9QNPEmQriuXnMpc+8UzuCJCHVy5ZL1a+xD/wLTzbTwnP9hO06stznm2fspZ0+8IM8MATCY9kvpQnG/mo7DeAiuvoNtiHK47Qqx270HkX1HWmUd7W6AlavHRlQmt885bt9MhTbZznQKofb9GOdu/5MeGzSP/6ux5xnsX73+3Sg44fP5HwWY1QcwhUI0A+gwjOdGfEKYBBJDNYJNbs8iBQ3DIIl/LyNHczVd3lftBBl0LDS+a3ovsXESNAH3YhGjGSaOSoyIbDX+1oFHTJR90ijUNGCdATz8awlHKF6wj8xrqD1vc/QcPbd4xEeR394VK/pc2a544ARdyzLT378ELCEDDeyAjhUvlNPyquc0PkvwM0/LtJDglhgxBh89Zv0FEfb9DZsz/Re12rvqsnn/2oe2/6+eefPa30iiNH6dFnqjaPkc9OmjbH8zm9UfMIxBOg2rFBEZ7OHYIwNlgkl0E8Eb8MAgNGVjk7P9OG6l51O23q1Td0eJjBJqu8dvnxO0aA7o3Id4DgAXq/M9Gw4bFz+Iiqa/wvAr8rB31LrRs3ixkH7TpErvxO3Q8bTkWXZHAKXLNnYvIDGYrY2fmpF2g4ZEmUu6zfAGp298O06Yuv4v4v44TxurTp07kjQNLzCA93hH+XdOludFzn5i/GcIgIHpW9v1ECRESdP+3lSoDq/+Nx2rFrj6fF7UdiWr/6ru90uE2l26juLVUfGpcECPnRI3cIVCNAclAE1xH7XTZoCDWsf5+ZtVE5fGTkys/1XfFJd1p4Z2QIUCMqf+9DoiFDiYYOi1RY8sHHrh0ClHTDm+6lsm/6RQoP1H9GCVDT5pHDD+1I5Spej5Q+/lTuCFAEPdubenxNze5+hMr6D4rz6Jec+5Cv8UpGCJfKr3tHngCdPHmK2rz+gWt/B4ICouJ1YM1Pw0dauj6rBMgLtfz/fxwBuqh2/KAID45EKNzUszfVvfI2MgNEESq3GRQU5a34+NOIEaBO7xMN/lZPxYCKLrkqc5sgPP400bdDYnKFkE/gzNchu48F160bVX19Wo54xoh1Qyrr1Se05Xer79LHnswRAbpMeHMtz/aw8P7e9HlPqlv7Vur85AtO+cv6DjCDOmZgp29/ohCXP+bRrqrfyp5fR54AYZraB5985UpiEnmAsHFB48dbuT6rBCj/iY5XDqsRIAyKRPiMWwYRYRwqunwSIQJ0VyMqf+c9ooGDqs5Bg6uu8X/9HRk8MkqAHnsqhhvLj4aRxKP00RwRoL9cFkkPJDy5nZ+obrB2frJVJPGo/PKryBMgGMF9B49yJTFYo4Npbl7HqdOnqf3bXV2fTbQG6GBZOd3b5DnXZ3UNkBfi2fl/NQKkyyCMnJZ0/sQZODKDKRFYBiKXQVR89HHECFDHTkT9B+ipGGTWA/ToEypTKlNU2qRZbjxAIEDq1Y08BpWf96Diq3QThP0HyuixZ1+JIyOX1/sHFS38IaHFvXTFmrid3ODNvqVhU9+pc0gUnqdR46bS7y+qF/fe515+22ynnfDFGqHGEKhGgLAMIoJn5YBB1Pq+pnHyWfJB10higfqv+LALLbwjKmuA4AF6822ib/oS9e0XC3GtvyOJR9HFGZwC16RZvByxfGkYKVxKH26aOwIkPdt6He/Jjggeld0/VwJ0zoz+6aefzLd4sK4HJ7w7QQ/+fhCeS/Z7Pvz9IDwLb5PfznFB86Px0kOgGgHCYBFmachBI/0dOTwq3u9MC+9oGJFd4ECA3uhA1LuPnooBZZQAPdJUZUplikoffjx3BEg9cJH3wlZ++hkVX3V95LfBTs9c1qfDhkA8AaoVycGRuKUfERkQSlTmivc+jBgBev0toq97VZ29eldd4//6OzJ4FF18Je1au45OnjyZ9CgdRvVOnTpF+/fvp5KO79KWhx+P4abyExn5MXrEqu/Shx7LHQFST3YkPdlmBsO5mQyVn3yqBChs1ruWJ20E4gjQ32rFBkr69dcQg2YRxqGi0/vhJ0CVlZW0Z88eKr7nfpp/US1acMU1585rxTX+F+7fxVdcS+/+7Soqjmj5q9fvNbR73fq0CNCBAwdo0Tvv09w//t2SJZaxaITFl9clnFVtKxrltss75w8X0cqun9KxY8cIU3Bq8kD6R48epU2rV9Pcv1yqHkj1QFJll4+VANVko9O0CxKBagSIlz5oGL8UJGJ4VLz9bvgJEEb4f/zxR1o+YybNHDGSJg0cROP79adxfftF5hz7TV8a0+cb+uPf61GLp1oRfkep/FxW1Dvqf/ao0bRq9hwCgYEnJ5V52qdPnzZzRzds2EDzxoylyYMG04T+AyKHK+TqvvoNzIlrxjoqIWRq4oCBNH3oMFo8ZSrt3r2bTpw4kV0CdMGl9NMXPeinHl9qGGEcjr//IRVdWU+nwBWkma6ZrikEqhGgPt/EDxhF8PdPmL3Ag0YRLD/KXtHhHVp4e8jXAMFQPXz4MG3fvp1WrVpFS5YsoUWLFkXmXLhwIRUXF1OXbl+a3T/Ou+BqKioqopKSkshgIOsb9b9mzRrauXMnHTlyhM6cOZMSAYJSxUg/DN5169bRDz/8ECk8IVcLFiygESNH0y9//Sdzjho9xvwP9yTmYb9evHgxrVixgkpLS6msrCxlUp2MAcAeoC3r1lHxrQ1ozp8vieQ54w8XRbLcbvU9+w9/p8WPNDMyiH5PD0VAESCSBGjeXy6jk527VJ0fda26xv9D/ruycxc68eFH1LFhE6r88KPQl9evPg+8+DKV3HxneDdBQOOHoQAvEEgQPEEwWHft2hWZE8Rv48aNdEW9e5ztD9u89i5t3brVkIAoYYGyYjok1u5UVFQYQzXVqUrwGsHIAImCJwnpRglLyNXmzZvpgcda0j//19/M+eDjz9GWLVtox44dkcICOmXfvn1GkcL7gw63pg/ILd6Fdy9dupRmzpxJkydPpokTJ0binDBhAo0fP56uqHuLCaNSbrdyTpo0iaZPn24GtUDC0ddhYEcPRUARqCJApStXUfFdjWh+rWsjec6rdS3NrXUNtfpLbWMLTrn8GsL/oorH/EuuoqVt21N5eXne68tfpNqQYajCIIGxivVAMBqOHz8eiRNrBA4dOkTDR493yA++a4Bz585dZg1BVLDgcqL+QYghDzAiU5n+xrKI52FoIL0oyRY8XzCyFi1ZSv/ym4vjzqXLVxpyiTiMedhDyBTqHzIFXZOOTLFsJQrxDkzfRPvetm2b8WouX76cli1bFvoThA+e3EFDhtM//frP1LVbD+OBjULZvcq4cuVKMyCBQT7IY6oDO4nkTu8rAoWGgBwsQvuZNWsWTZkyxQwYYdAoCicGSTB4MnbsWPqP82MEqHGTZ83gEe5FAQMuoDnY9gAAIABJREFUI+p+xowZhJkb6DsxiJ2NQct02k3KBIhfCoMBDSFKJwwyGKo33d2kGgF6rm1HY7ij4qOECZc1U0ZqVOUKSuOx5m2ryVXTZ18xxAfEkLGOSpgpmWKdlSgEriBeaOPwQsL4hScq7Ce8XvBA3nlfU+N5vPSaO4zXEV7YsJfdq3zwamMkE4MNqU7rTSRvel8RKEQEoJcxWIT2Ab2BKesYMMC05aicGByDwf9ahw/j+uwpU6ebGQS4HxUsUE4sg8BsFfSb6EPRl+bzkTYByufC1UTeUKHwTEybMSdO4NkDhHDL1u2m4rNtuNVEeTXN7CDAncmKVWs85WrNug1qhGWhOlAXaOcweNHBo72H/URnBfK9YOGSOM/j51/1Mx5t3A87Bm7lQ/2zBzLfO/MsNA19hSIQhwDbQ9AdWKcJwxeDBlE4MTCGwSEY/P/99+vj+u0mT79kpqsjThSw4DIePHjQDBzCW57v3h8IshKguObs/wOGESoVU+AaPPh0nMBLAvRo87bOVDD/FPWuIhBDAHKF6W2PP/Oyp1zBC1QIoyphqVO096icIHuQvyZPt4mTvz9cepPxfrBHOyp42OUMi0xrORSBTCKAdgISBP2AgQIMGETl5EGjnn0GxelMtgWXr1xjdCcGVqKCCWSAZ6lANvL9UAKURA2hoUOQ5xaVuAo8Cz5CCD+UQiEIQRIQaNQaQAAyArny8/6wbKkXqAYqIOJJQv7Qca1eu95Vr33wcU9zX3VZxAVFi68IeCAA3RClE7YgdCY2ffr71Xe66s0Hn2htBizZDowSPh5iknf/VgIUsEogvBBkjJI+0LSVq8CzkYrwnoeeMQ0EDUUPRcAPAZYrP+8Py1azFu3UC+QHpt5LGgHIH9a4NH/hdVe99h//U9vchy6DHtRDEVAEFIEoIwBdCM/OwKGjXXUm99eLf1ih09bzWFCUAAWsHHT8YPxLlq3wFXgWfIRzixapFyggvlGNxnK1cvXawHK1bsMmI1dRxUzLnTkEWP42l271lb/2HbuoFyhzsGtKioAiUMAIYJoXlkLUqX+fr97EQDiIkg6E52dlKwEKWC8QYMz5HDtxGvUfPJy69+hF73/Ujf737y5zGsAV9RqYbQ/xMVR8DwgLA9FQdNQ0IMgRjMZyNX7y9zRwyEj6vGcf+rBrd6pVr4EjV7Wvv4e6dPuCevbuT0NHjqXJ02erMRpBWamJIkP+sGC19audHHmTgzjy+uDBMh3QqYlK0DQVAUWgYBCAPYcp66PGTkqoM2MD4QvVDszT2lUCFKBiIPA8TQS7fmCrx7lz5xL2ef///vsKpxFcdu3dZi98fE8D+6Dzh/OUAAUAOaJR2ADFbjHr16+n4uJimjZtGl15wz8cubrqxobmg5z4Rgs+yIhv1EABq1xFVGgyVGzIDwZo9u370ZE1SXjs65ZtOqjcZQh7TUYRUAQKEwH02Zgy7PYZFFtn4vetDZuqFyhPq1oJUMCK4ZF6bPOHbQ+x5/n8+fPp13+40jEeQIBAfmCkYltA/XZEQHAjHA1yBTIDsrxr1y7zLQV4EOvcdK8jV3Xr32fkatOmTbR3717jelfPYoSFJkNFBwHC9IzX3u7iyJrdgYOIX31zI6p3+4NU/57HzBb/GAxS8p2hStBkFAFFoGAQ4EEjbH6AQe5FixbR1KlT6bGnW8fp0A+6fGoGLWEnYtActiD0ph75hYASoID1wYIPQcYI/M6dO40n6LwL6jqCf0nduwz5wV742CwBa4bUUAgIcESjQT6gGDG9EkoVBAcfEwPpYWP02lvvN1+jh1xh3jEIE4iTHopAOgjwSCY+gIoP9uEr3p0//syRO8hfkydbGY8kvI/42CHkT8l3Oqjrs4qAIlCoCLAdiOUNbAPOmzePmrdqF6c3u3brQQsWLDADmvjIshKg/KxxJUBJ1Asbqxg1BQnCiPzvLrzOEfyL6txhRvFBftRISALYiEeFXElP0ObNm+maWxo7cnXdbQ8QjFSWK8TFM3ooAqkiwDKHjhkdNL7ijg77sx69HLkDAXrs6RepqKiIVq1aZWRQCVCqiOtzioAiUOgIsA2IdZOYDYRBIejG1q90iNObn33Zx+hUzOrA7A7dCCE/a14JUAr1ghF7jABgowP5BeALr7zNuDv1Y5UpgKqPGE8QDEzIlSRAmH4EI1XlSoUkkwigM0fHjA4aBBtr0AYMHhbXkTd9tq3xSGK0E4M+kEGdApfJWtC0FAFFoJAQwAAkBrhBgqATMRXu1Q4fxOnNnn0GGlsQszqgY1Vn5mcNKwFKoV4gzDBUIfh/vPQmR/D/fHl9hwDpCH0KwEb8EShWeHkgV24ECIpU5SriQpLB4vNoJuQKAzqYqz7yu/GOPoMH6ImWr5ipHiBJSn4yCL4mpQgoAgWLAOtOkCAMTr717sdxerNn74GGHCn5ye8qVgKUQv0wAYL78y+1bnEE//xLblQClAKe+kgMASVAKgnZRgAdOeQOHXVZWRlNnDzd0WcgQE2ffZmw9kzJT7ZrRt+nCCgC+YwAdCf0Jja86tS5e5ze7NGrv/Gs6zrwfK5BIiVAKdSPJECY9gZDAefvL6qnBCgFPPWRGAJKgFQScoUApnSUl5fT1OkzHX0GnfZo8zaGGOm267mqGX2vIqAI5CMCIEDQiyBAH3T9Ik5vftazr9nUSNeC52PNVeVJCVAVFoGvJAG6uO6djuD/9q/XKAEKjKJGtBFQAmQjor+zhQAToO9nznH0GQjQw0+0NgRIRzKzVRP6HkVAESgEBCQB6tztyzi9+ckXvZUAFUAlKgFKoZLYUMUUOHz7hz1A//dPdQwBwrxQNA49FIFkEGC50jVAyaCmcTOBABOg2XPnO/oMeu2Bps8rAcoEwJqGIqAIhAoBSYA+/uzrOL3Z5dOvlAAVQG0rAUqhkthQ3bFjB9W6/h+O4P/H+bXNbkpKgFIAVR8xazF0EwQVhFwgwARoflGJo89AgBo1aaEEKBcVou9UBBSBvEaACRDWSHb/sk+c3ny/6xdmDZBOgcvrKtQ1QKlUjyRAdW5u5Aj+v//+cocAIY4eikAyCLBcqQcoGdQ0biYQQEeNnd5KFi529BkI0D0PNTffu9A1QJlAWdNQBBSBsCDABAjfA/ri635xevOdD7ubNZVKgPK7ttUDlEL9sKGKb2Nce+v9cYKPD1/h44JKgFIANuKPsFwpAYq4IOSg+EyAlvywLE6f3dX4SbMLnBKgHFSKvlIRUATyFgFJgL7qMyhOb3Z47xMlQHlbc1UZUwJUhUXgKzZUQYCuv+OhOMHH/5QABYZSIwoEWK6UAAlQ9DIrCIAA4aN9y1esjNNnt93bVAlQVmpAX6IIKAKFhIAkQH36D4nTm6+/3UUJUAFUphKgFCqJDVWQnZvubhIn+Fu3bVMClAKm+gjpGiAVgpwhwARo9eo1cfrs5gaPmm1e1QOUs6rRFysCikAeIsAECN9P6z9oeJzebPfWh+ZDqDoFLg8rTmRJCZAAI+ilJEC3Nnw8TvA3by5VAhQUSI0XhwDLlXqA4mDRH1lAgAnQ+g0b4vTZjXc9ogQoC/jrKxQBRaCwEJAEaNDQUXF6s+3r7ykBKoDqVAKUQiXBUMU0N3iA7mjULE7w12/YqAQoBUz1EfUAqQzkDgEmQJs3b6b//bvLHJ123e0PKgHKXbXomxUBRSBPEZAEaMjw7xydic1jWrd7RwlQntabzJYSIIlGwGsmQNjw4O4HnooT/FWr1yoBCoijRotHQD1A8Xjor+whwASotLSUfvXftRyddnX9xkqAslcN+iZFQBEoEAQkARoxepyjM0GAnmv7lhKgAqhHJUApVBIMVXzrZ/fu3dTw4WfiBH/Z8pVKgFLAVB9RD5DKQO4QkATo13+8ytFptW9sqAQod9Wib1YEFIE8RUASoFFjJjg6EwTo2RffUAKUp/Ums6UESKIR8BqCzwSo8WMt4wR/8ZJlhI9ZgiTpoQgkg4B6gJJBS+NmEgFJgP7rz1c7Ou2KevcoAcok0JqWIqAIhAIBSYDGTpji6EwQoKdbtVcCVAC1rAQohUqSBOjBZq3iBL+4ZLESoBQw1UfUA6QykDsEJAH6/y+81tFpl1xzlxKg3FWLvlkRUATyFAFJgCZN+d7RmSBATZ992RCg06dPE+LpkZ8IKAFKoV4kAWry1Itxgj+vqEQJUAqY6iNKgFQGcoeAJEC/v6ieo9P+dtXt9OOPP9LJkye1I89d9eibFQFFIM8QgB0IgnPo0CGa9v0sR2eCADV5+iXC9thKgPKs0qzsKAGyAAnyUxKgx59pEyf4s+YUKQEKAqLGqYaAToGrBon+I0sInD17lo4cOUJbtmyh8y++wdFpf6l1ixKgLNWBvkYRUAQKBwFJgGbMmuvoTBCgh5q9oASoAKpSCVAKlQTBr6yspD179tATLV+JE/xpM2YrAUoBU31EPUAqA7lDgAnQ1q1b6c+X3+zotD9ddjPt27fPeIB0XWPu6kffrAgoAvmFgCRAc+YVOToTBKjxY8/RwYMHST8gnV91ZucmIwQIgoAOFB0kTlyH/cQmCPv376cu3XvRHfc96Zy8BghTSsKOgV3nkAM9UkdAPUCpY6dPpocA2jI8QCBA8PqgE8f5PxffoAQoPWj1aUUgkgjAHmCbMIx2IWw8DISD6MydX+zoTOhN7A6MqcOwE6NiC9r2rqz7fG0AaREgJj6oYDV+87WKs5cvbgCZkAWWLU4zCiFGiyoqKgjfYql7S2NHoV532wPG24iP70ZJmbICzYQ8Za8VFOab0L6YAF1U5w5H9n534XWGAKGjR33ooQgoAopAIgSgK8JuF6Jf4jVAJQsXOzoTBKjBg0+rB+ickEAW0L/kYz+eFgFCwbA4FqceigAQgBGfrqGE56FY0k2n0GoE5cUW6tu3b6drb73fUaj17ngosuswgEm2lCcUdFRPGCsg31gDdEnduxzZ+81frqG9e/eakUyuhzBjVGg6Q/OrCOQbAtAP6L8xaBLmg8uJTRAWL1nq6EwQoDvve0IJkKh8yAP6j3w7UiZAXPnoNMMu6PlWafmcH/ZSQD5SOViuQKSidsDY5+9LLVu2jObOnUvFxcW0fv36SCtT9nrVlDxA5phoIeRrKOyonBjEOnz4sPE+XnZdA6cz/88/1TEffAYx504sbJhwnUMOULZUdVdNyaemqwgUEgJoQ9AXR48eLaRsJ51XtlVAgJYtX+HoTBCgWxs2jXSfbYMJjpCP66HSIkDoNDHP0Rb0Dz75ivSMDgZS2GFEQS5SNSKiojwlZnwNQwxKAttnYiR+9erVtG7dOtq1a5eZnhT2KQWMgx2CFKYjU3Z68jcbvVEk3BIHtDueAlf7hoZOZ/4f/1M7UlPg0Aaj2s6kPOi1IpAKApIUoB+TR/hswp70bpce9NZ7n1DrV95ydCYIENZRvtnpY3q3yxf0/sc9I2kPy7oHR0A/Dv2aT0fKBAgFQYFgnMELJA8IgJ7RwABKTR5QemD7qRIgjDKXl5cbEiDTDZ/ydCfI73X9kt7p/LlRqm+804Xe6NSVOr7/Kb37UbQUqax7GObwLNaE8kSakNeoe7ElAbr65vsc/f3vv7s8UgQIcgeyDTz0UAQUgeQQQL+PwSRsEIXBcXnc0fgpR6+Exj78zSX0L+f9nf75//41Vjb8/s3FsfO8i+hff3NJ+MocwL5HXcsDHAFewZrow+V7kr1OiwDBKNm5c6cSoAACEZoGb5XVJkDYESVVAsTKE2lAgcojlMrTwrJKRi6JKU4oTz4944aPaNvKE15FjCBlWnlC3njtC94R5QMGPzA2689ue4Dq3NyIrq5/H8n1Z5nGP1/xBuHOx+ka+YqX5ksRYAS4Dwf5wfb58ghnH84E6EL6p1//mX75f/5Iv/w/fzKECMRICVBMAtC/hpIA7dixQwlQhIzTKiM9ZnjXBAE6cOBANEaPIiw3thzJ3zYBgkewpggQRvshb26EG/mIynn7fU+aees33vUwXVP/Xqpzw9109Y0N6LrbGtPN9zSh2xo1Cy0Wtg5jL3ZUCJ80UvVaEUgHAUmAsHmKPEJLgODx+a+LDOn5p//8C/3zf/6V/vm//ma8QP/6W/UAQQbQhysBUoMvdO5Q23jIhAcIxmg0Ro/C572RRCbVa5sAYZFpTREgeCvRUduddap5L+jnzHQOdOQXUlxHHvKRTFuHQf/U1JRLaRDqtSIQNgSYAKH/tnVqOAnQpcbLY6a9nfd3Q4T+5byLYuQnotPf0AfafbgSICU/oSM/EHTbeFACpKQmXRJgK8+aJkB79uwx31mSxki6ZSjM58/NX0dHjk7cnOGfxmHrMEzfUQIkW4NeKwLBEIgkATK2rZiubohPND0/3O/ZfbgSICVASoAS6FBWnuoBijaJspVnNgjQ7t2746STFXn0wuh15EqA4kRffygCKSPAfXikPEBq21azbe0+XAmQCkk1IQmDcWUbD+oBijZ5yYRM28pTCZDKVCbkyisNW4epByhl+1cfjDgCSoBUV0PP2n24EiAlQEqAEnQOrDzVAxRtJWorTyVA0ZYHL+KSqf8rAUqgmPW2IhAQAe7D1QMUbZ1t9+FKgJQAKQFKoERZeSoBUuUpRUUJULTlIVNExysdJUCytem1IpA6AtyHKwGKts5WAqSEJ5SExzYibOMhyBS4M0d2mA+lQlnKg5WnEiBVnlIulABFWx5snZPp37YO0ylwsvXptSIQHAHuw5UARVtnKwHKQwL0XtcvqbLypGtrPnGikhYvXUkNH2npSly+7PMtnT37k3n2UHkFvdj+Pdd46JzXrNtk4s2csyAuzvV3PUJbtu0096AocP+8C+rGxeHO/f2Pe9LJU6eo4shRevL511zjcNxchrbxkIgAGQV5cB2VzWxFJ3YVxxEhVp5KgFR5ykaqBCja8lDT+s3WYUqAZOvTa0UgOALchxcaAUrHNly6fI0BCGVftXYD1b6xkau9BjsO9hzsOth3Uq+FzTZUApSHBKjzp73o1OnTxug+fvwEHT123Jxnzpx1WjgEtP3bXeOEE4I6cOgYJw4ulq5Y40le1m8sNXFnz18Yl87dDzSnPXt/dNIBGev4wWdxcbhRcF6Rx+at33CNw3FzGdrGgx8BgoI4c+YMnTi0nfb0/aM5949vREfXjyR8dBAnvsBeaATom0EjjUw5FXvuAuU9cvQYTZs5n25q8Fi1OsyU0ssEOc+lDNnvtpVnsgQIuB+a05Yqd86LI9hu9YPvAGEbbN0FLroky9ZhSoDslqK/FYFgCED3og8vNALE9hbyn6xtyPYeEMIgea/+w6r19ejjYMfBnoMNivfJfi9stqHdh+saoDwgRCzkNqmAF6bPwBF08uQp08o3lW6jy667J05AmQBBeCHkbkLMAs0NwosAwdDHKAAOeIvsdyEdr7zyO/IltI0HPwKEcsPgPHRwr0OAmAjtG1qHDv/wKR0r20EwQKBA5YEGlS9ltvPBsgG5OHaOVCNkjyHKAeL7eIt2cWXIlNLj9zNeqZBzu0y5/G0rz2QIEGQMJLts/ltGxn787k46sqKXIdfo3HDygWslQNElPizjtg5TAsQtRENFIDkEoFMLmQClYhuyvcezi9DXo29n/cJhEAIUFtvQ7sOVAOUxAWIBnTpjnmntMF5btu0YJ8BsZJZu3UHbdsS+GbJ5y3ZXdyc3CC8CBPI0b8ESYyCfPnOGvug1KO5dyE+hEqD9+3bRsUM76VT5Fjp1YDWd3LuIKnfMpuOlE+nIumFUtrQX7Z7XuRoBYiKEcPfUVrRn3Yw4Y7UQCJCt+P5S+1YaP2WmQ4SKFy6Nq2cmQOkqPZbNdMg5t4F8CG3lGZQAofMF+cFHLPeXdKsmY2Wz29DJg+sduSpUAoRpFm4HCPe+Hw9Qj96DXb3T6nF0J3tKgNykSf+nCBCdPbojKRjCRoC4P/SzDdneg00HEgQMvhs/La6vRzpBCFBYbEO7D1cCVAAESBqStouS70HY4eKEsYHTzd3JDcKPAH3+9UCHSIFQ2VOkCpUAbRzxIG39+g+0+5vYFDdJbJK9PjCpCR3fPMEolEIkQFB68C6uWL3edCIwTu95+FlHMTIBSlfpsWymQ85Z0edDaCvPRAQIHQ5OQyRPnqSysjLaMb97NQLE8gev0LFNY+ns2bN04sSJgpsCx/oFXmSMWOLEGkZggAPh/JIfqnmWWd7Yokl1Ci7LG6dT6B5HJUBckxoqAlUIQI+cPLiWDs1sSSe2xw9IVsWKv8IzYfIAcX/IOg99tW0bsj4eNGysWRoBRLBOvFnLV52+HukEJUBhsA3tPlwJUAEQoNHjp5rW7OcB2r3nR3qgWWuC9weHPeoPQecG4UWAoCT6DR5lPD/wAMFw+3bk+LjGUqgEaNOYpzNGgPaOuo3KV/SlM6dOmA9rsTLKt5CVo5ssIK+QAxy2e50NUijVdJQevx9ylyo5zydMWXminaBtYFpl+cE9VHlwE1XuWUjHN4+nI6u+ocOLOlPZnJfp4NQnaP/YBrRv2DW0p98FnsSHCZAMDxa/R7u2lxbUGiAv/YINXDZu3mpkDXqlW49+cTqF5U09jvGeICVARmT0jyLgIADda7zpZVsdfbpv1C1UsbynM53YiSwuwkqA/GxD1sfY1Kpdh4/MGiJAYtt/QQhQWGxD7sNZNJQA5TkBeuiJF830EVSY27ocNjLZyGWCAoG13Z3cIOwGwAYI3oH0sPaHd4zbf7CMHnmqjWOwcPq20ZxPhiryYhsPG8a/lDYB2tH/Mloxph1tL11DaDgYUbr9vicdbPINA1s27PwtXLLC6AEvD1C6So/fnw45t/Ocy99QnjDSK5Z/SfDW7Pn2StrT909ORywJTLrXeyc1pe2rvqedO2O7M7LCzmX5E73bS7/gOew0VH64whRjbvHiuDbD+kc9jkqAWM41VATcEGDv+MEfd7rq3bLZbc00d+hp9F984DpsHqBEtqGtj3nAE5spgBCxPg9CgIBjGGxDJUB5QHhY8DhkUoHNDiZMnUVDR02kEWMm07KVa50NEEBE7PU/eJ6NTCZAmNqEqR84bHen3SD4/WyAsJDj/9gFjueNTp4+x2ksnNdCI0DrJrxCawbfS2uHPkyrhjWjFSNa0NKRL9GS0a9RyeiOVDT6A5o1sgvt7Pc3V8W6ZMizNHFUX5o6dSqVlpaa6UwnT54sWALUvmNXs/Ul6rymCLEtmyw7yZBzltF8CO+470k6ffo0la0c7Coj6ZIePL9j8DW0eX4v2rZtG23evJl27Ih9m8o0aCKnHeYDHnYevPQL4kkdg3jyWb6nHkclQCznGkYPAfQLbifIDM6zZ89QZflOOrSthHb8MMxXB+8f+w86tnW64xUqdAKUim1o62NMfYNNiENOD06GAEFvF7JtqAQojwmQm8rDfHqQInstDhsQtpGJ/3u5O+0GwWmwAYL3Iz38H0SqZPFyk6XDFUeoVbtO5v9sxBYaAVq2bBnNmjWLvv/+e5o+fbo5cc2/p02bRhMnTqSt/a6IU6zLBj1C08f0oSlTpphz/vz5Zm3GkSNHjDF8R+P89wBh5B3eQBBruMw3bNpCvMU6pkw2evQ5V4NUykMqSs+WzVTIOctoPoTw9mEjg73LhsbJiBfx2d33Ato+oDaVDrqJNn57F6399l5aOfhBz2eXjnyR5s2bR6tWrTIke9OmTaEhQJiei0EcHF4eIBgp6UzBZXkLi8fR9mLrLnBGfPRPSBD46ewpOnN0N1X+uIKOb5tBR9cPp8PLelD5gk50cOYLhLW2P46+jfZ+W9tTZ3rpXv7/3gEXUcXyrwwRKmQPkFuVJ7IN3ew92AGGDIo1Q8kSoEK2DZUA5TEBkiwfi4V5+2sv7w+MQu702QPEhqKbu9OtQSC+GwHC/0F6QH5wgAxB8AuVAK1du5aWLFlCS5cuNefy5ctJnj/88AMVFRXRlkHXG2W7buBt9P3oHgRiBKN08eLFtGjRIlq9ejWB/ECZYlQKDYoxz7eQZcNNeYLAYnEkdoSz8+0mD6koPX6/lM1kybmdt1z+BgFC3e9YOoq2jLyfVg55hJYMaU6LR71KC757j+aP+ZjmjPmCZoztTdPHD6apE0bR5MkTHfIMgj127NhqnfnqQQ1p+riBNHv2bCOb8Pzs2rWLtmzZEoopcJCdUeOmxowQ0fFyXdryluoUXFveWFehww86HZjzlA+hEiA3zaX/CwMCxgg/XPXdPSYsmQ73zWhDRw9sMYOVmLFRqN8BSsU2dLP3pK7l3YKTJUDQjYVqGyoBymMCZHtVXmz/nkNAsMWs23d57E6fO243d6dbg0B82SiQHqeBENPfoKx4ZyY2Kuy8ymfy4do2HjClaPv27WZROT4wCUUoT6y1WLduHW0ZehcVj+9uprrNnTvXGKQbNmwwz8Iw3bp1q1mECfKDoxAIkPQALV+1ztf7g7rzkodklZ6XbCZDzvNBljgPTIAgSyDMEyZMMCc8ijNmzDAeRpCYOXPmEGQH3sLi4mJasGABlZSUmN/wOG4bUMuQoB39L6HZo7oYko14INeQ0QMHDpgplpDJQlwDhA0P4G3EOeX7uWYNI3QIdqfEtq0gRIypl7ypx7H6Okb1AIXB9NcyQBdgI4OjFYeqDQZligBtHnInrZn/rZlKjN068V01nIVKgGx7K4ht6GXv2RsSpUKAoLcL0TZUAlRABAhChl3YYGxjBODtDz+PMxxw38vIxD12dzJ58WoQXgYv0sAGCDx1BfNGP+re23xs1W6QiJtPp02A9u7dazYuwBQmVoYYEcKJ3xjZ3717tzFC4emBpwjEByPx2L748OHDZuobyBMUOB+FQICkBwZ11PXzPkaeUA7sEGPXm588JKP0vGQzGXJu5y2Xv0GAsD01DFFMT4PHEOQGnkSQF3gZQaIhN7iPtWIgzCBMCPF/EKfSb2+lJcNfMCQb3kVMzwS5RucMOYQ8Hjt2zMgjZFIeuSx/onezfpHVe2UVAAAgAElEQVT55Wtsqf96p0+qyRrSdJM39TgqAWLZ0TBcCKDfQb+LgZ5dw272JEE7+19MWwddSxsH305rh9xPq4c9TsuHPU0/DG9FC4a2obUD73B9dsGI181gFPpw6F0QIO7nw0KAoDcT2Yasj+01vrVvbOTsFgy9/Mpbnc1OsG5babvpZu4HCtE2VAKUZ4Y6hMnPq4L1GfsPxObO2x+sxLNeRibuSeHFN1+w9gOH3SBkPKTHAs6hbGjjJs8oSAIEZQvjFWQSClie+B+mtIHkwOCEwYoQxAeECfdglMLwhQKVRyESINQre2HgHcIOXVzXCP3kIRml5yebQcm5zFeur3kTBJATyMbGjRsNyYFMsNcGuwNCjkBkcB49etSQGVzv37/fEJ0fFs0zniFMyQQpgqzhGcgYRkax0xGuQbYLkQCxBwjrzXbt2Wfampf3x0/eou5xtAdx1AMkNa9eFyoC6G9Zv22a2oGWjX+Xisd+THPG9aSZ4/vR9AlDafKE78ya3EmTJhFOrMHFBkQ4cQ3v+6pB98URoB8GN6VJE74z5Ae6FYNQaDPQ19i8BiQoTAQokW3oRYCgc9nmRF3ApsOAdrIECOkUmm2oBKjACBCEjL/262as+hmZeJbdnRBuCDmOZAkQNmDASAGOiiNHzRSqQvMA4ZstULogPvaB/0ERQEHCYIWxysQH/8cJEgQDNiwEqM3rHxC+K4Wyjxw7JTABSkbp+cmmJFl+5Bzvy5cTyhOyAJICeYAssFzhNzpZJjAgMThZfvB/yBQ6ZKztgccH09uYZONZxMWBOmEDoRAJkNQv8OTMnrfQlAlebBj1dn1KWbAHYKLscVQCZGtq/R0GBKDnoN+gC+E5x/RfTB3GJkUIefowPOy4h+nDmJUBUgMPOtbjwnO+engzQ4A2D6hH00d/FTeVGLtnsm6F7sU7WWdjNog88nkQk4mKl73lZxv6ESDoZd4tmMlPKgSo0GxDJUB5ZFCxIZBIyDH1DcaDm7HqZ2Qifenu5EYvDRTE8TNAOI9f9BpE+IghH14NkuPnOrSNBzZUgaHbgf9DScJoZSMW/+MzbAQIChDEA8eOXXvidhlMJA9BlV4i2QxCznMtR/L9rDxZJiBTIMuQGf4fh7aM4f8gORiNhLeHP8AGueLn+RnEDQsBAn5SniBzkD2Jq7xvE6AoexxtHaYeIG4hGhYyAqwLoQcxCARPOqYOu00fxppIxMFUdPaI439muvHYF2nx2E5mR1cQIkxFxrRjtBPoZalb8c4wEiA/29CPAEH/yg2JIE+pECCkU0i2Iffh3H64H+bBR/5/rsNfpJoBFAQjrRgBqKiI7XnOaclON5+uExEgbH6wqXSbKYZtrCYyMlFOTp9xSIUAyZ2ZkE7YCBDKBCUpT8aLlWeYPECQi28GjTSkD8S2W49+jlHqZ5Byuwmi9BLJZhByzu/Lh9BWnphbDo9hUOUJOQK5RkcMMsTEB/+XB36HiQCh7uBlRLngdYT3UdZnInkLOs3CT97kOwrF46gESLYKvQ4LAtAD0JmYcQGiAj0Kbw1COYUYuhW2HE5MX4dORAi7Dp5xkCB4hLCGEiTKnkos9Squw0iA/GzDRAQIOhi2IB+pEqBCsg3tPlwJUB54hJig+JEKNlaxocFr73zsGBDc6ePbF7c1esL5vzQwMOKKxe5QAjiHjZ4YF4+NA9zDdzjks/IaO4/wh7T27ttf7fsxMm6ur23jIZEHiJWAW8jKM2wESM4hliPzLA/AAvLlVpdBlB7Lpr0Jg0yPZZ9xt8m5jJvra1t5JkuAUEZugxxyuWWIe2EjQFhnhim8ODBtQ9ZlInmLqsfR1mHqAZKtRK8LGQHoOJAgDAhhMIhP/OYZGBggQhycHB//A5GB9wgeId7Zlae7IR2OL/HB82EkQNCjXrYhEyC3jY5Y/8IGwKA6jiNHjxFsPL6HkHUz8AuDbWj34UqA8oAASYHT68ys+bCNByVAzeMUG8sZzyEuO3TY7PiXSaXHBCgdcs75zIfQVp6pECDZKXtdo7MJGwGSUy6xqQs6Xq5T7mSBhxfhjqLH0dZhSoC8Woz+vxARgJ5LdNrlQnyQIHiPMJ1Yrtf18qgjDTxXiASIdaSGmbEL7T5cCZASIMcQCVMjs42HKBKgMNVnPpTFVp5KgOI7JR5x9PLiYZolpltihBajllynQQhQFD2Otg5TAmSbw/o7igiAzECHSG8REykvPJQAxetq1r1RC+0+XAmQEiDHEAlTY7CNByVAqgDTlW9beSoBipepRAQIU9l4usXS5WscvcMECEZKOtMswuZxtHWYEiAv81b/H0UEEpEeiYkSoHhdnW5fWKjP2324EiAlQI4hUqhC7ZZv23hQAqQK0E1OkvmfrTyVAKlMJSM/yca1dZgSIGnS6rUiEBwBJUCqq6F/7T5cCZASICVACfQoK89C2wQhWYNL4/t3ErbyVALkj5fKU3r4KAFKoJj1tiIQEAHuwwvtQ6iqQ9PToTZ+dh+uBEgJkBKgBEqUlacSoMwqI1s55ftvW3kqAYq2PNS0vCoBSqCY9bYiEBAB7sOVAEVbZ9t9uBIgJUBKgBIoUVaeSoBUeUpRUQIUbXlQAiRbg14rAvmLAPfhSoCirbOVACnhCSXhsY0Re/RU1wBFW/HZ8pHKb1t5KgFSmUpFjoI+Y+swXQOUvwa25iy/EVACpLoaetfuw9UDpIQolITINh6UAKkCDGp4esWzlacSIJUpL1nJxP9tHaYEKL+NbM1d/iKgBEh1NXSy3YcrAVICpAQogd5m5alT4KKtRG3lqQQo2vKQCZLjl4YSoASKWW8rAgER4D5cp8BFW2fbfXikCFDAtqLRQoiAeoCirfj8DM2g92zlqQRIZSqo7KQSTwlQCDsiLVJOEPAjQDnJUI5eChzwAdnmrd+g5q3fNB+Uxf+ieigBimrNR6zcNUWAIgajFlcgoARICVAqxCboM0qARGPTS0UgDQSYAGEaKbxAUTyY/GzYWOrM+tm4eSudPXuWokqC0IcfO3bMEMF8kolfpJqZn376iU6cOEG7du2iHTt2pJqMPhciBE6ePGmUXmVlZUoNnZXngQMHCApUD0UACGD06OjRoxlXnpA3yOqePXto9+7dcWDDC6VnNDBQAhQn+vpDEUgZAe7DMY19y5YtKadTyA+C6Bw/fpyatmjnEKAnnmtPsI9gN0fxABkOFQFi4wEF27x5M+3cuZMg9DBco3Si/GvWrqe58xcY4x+/o1R+lBX1DiIMhYffaOiQj2QPPHP69GmCFymqypMxg6IcPHysOVPBktMJQ7h3794aUZ6sw9wIUBhw0zKkhgB0GAyYqBorqaGmTykCZPp99OFlZWW0YcMGMzjOA5pRsItg/8EWnj13vkN+2BM9v3ihGWiLko0IWw78APYhHCb5plNT9gCxsYrRWRRw4cKFNGPGDJo6dSpNmTKFJk+eHOpz0qRJNHHiRBoxYgT9r1+db84mzZ6l8ePHE+6FvfxcPtT1zJkzadGiRVRaWmpG66EAUzXaMW8WMgXluWLFCoqS8uQOAgoSCuMvtW41J7wTUVKawAGkmpUnOpSaMEglAQIJ0qMKgW9HjK/6EbErJUARq3AtbsYQgE5FHw6PPWYGLVu2jGbPnk3Tpk2LhF0Im/C7776jWxo8Uo0A3faPR2ncuHGRsA9hF37//fdUVFRE69atM3ZcPnrA0iJAcPVhCgkMlW3btpmCrl69mlatWhX6E8b5kiVL6NGnWtE//ecFzgkisHz5clq5cmXoMUA9o77Xr19P27dvN3IAeUiH5UOm4CqF0Q+M582bZxoSFCjIddhPKA4o0TbtOzoKtO1rbxuliXthLz+XD4MpxcXFNao80VnztM2NGzdmzAgo5ISACdrg3+rcYdoxfkfpgP7K19HKKNWDlrVwEUD/f+rUKTOQicErDGauWbMm9PYQbL6lS5fSqO/GOn03e384HDt+YmTsw7Vr1zqzgjCAiX4l3/qTlAkQmid3ljAiwPgxcg/XJwhRmE94JTAtZ+Wq1fQv510Ud7Z5rZNZU8Aj2GHGAWVDffMaDchBukIO5QkPUkVFhcER3kWMIECBRuGEEl28eDH98dKbHCWK6x9++MF0IFHAAIoTnebWrVuNJ6imlCf0FzpqyDHeCcINMh8FjN3KiMEMyN/7XT43svdWp67mt1vcsP0POga6BnIA3Z7uQE7hmq+ac0UgfQS4H4fuPnz4sLETwm4bwnOMgeAHm7Zy+m4mPhw+/ERrM7gLGzLstiE2PoAdB10Kr2C+kR9IeVoECAmgUBB2GL4wXHHCqAjziQpFxbZ46Q36199cUu1ExWO+IwhBmHFA2bjOUf+Qg0wIOStPYAicQbCAKRRomE8oREzF+uSL3tUU6KdffmOmwTHpDDMOqGvU+ZEjR2pUeUJWeboGRipBgObPn0+zZs0y03nhwo/SCS8rvI/nX3yDkb/zLrjaeB7x/+nTp4cWC3gb586dazz6mZjGm775qCkoAoWPANuG0LFsJ4TVHmInwOIfllXru5n8cLhi1ZrQ24dc36j7TNmFNdEi0iZAnCkIexRONs5Lt2zzFPQ33/3E2fEjCphwGVkWMhEiTWCNBiQVKDessIXoGED4QGyw9oeVJYf4H8ggFG3Yyu5WHq73mlaeSB/YA1uM/GPzjU2bNhGmxMELFZUTni9Mae3U+dM42Xup/TvONNewYoG6BvHB1DeQb/b+QAfpoQgoAqkjwLZB2EP0V5gF9VjztnH6k/tvGTZr0S5Ox4Qdm9Slp+afzBgBqvms5scbYDDBUG3+wuu+gl5T03byA4Xs5SLsyoHLB7mC4dXzm8GectWr/1BjrDMp4GfDHGZD0oAnSBjaNbxOmLKBE16oqJzwLMIL9ufL68fJ36/+u5aZsoH77JkLGyaoa9Q71h6CDKc7jTcbMqvvUAQUgfxAAP0v9AY8O5Lo+F2v37jZDPDmRwmimwslQEnUPQQdTH/T5i0JBf2t97qZRoFn9FAE/BBguYIR5ub9YUX619q31chuaH55i8o91AGIEIxf9j5FJUTnDeO/R++Brnrt1Q4fGbkDSQwzJkx8VGdHpdVrORWB9BFAv4HBs3GTptOAb0dQ9x696IMu8Z702tffQ126fUFf9RlAw0ePpynfzzG6VHVN+vink4ISoCTQY0Fv1rLqA1dsnLqF6gVKAtwIR4USxNQ2P+8Pyxe8QDBEVXHWjMAA16idIDWYAuhHvg8dKjcdNnRgmPGpGanSVBUBRSCsCGDgBAQIn6rAJirY+hm7mXKfjfCqGxuaz4VgMyNMs4YXXfvx3EuEEqCAdYBOH4bC2vUb6aa7m9C1t95vhBrMXgq6vO7w/qcq5AHxjWo0lqtE3h+WK3iBoGxhiOqhCKSLAOQvCPl+8dV31aOdLtj6vCKgCIQOAfTF0KEgNZhGDBK0YMGCOLuwbv37zDeRsNMkiBI87rAnoX/1yB0CSoACYg9BxVQRLFLHYmCwfHzwtJYgQP/+u8vNLkr4dg22VcW8eTyjQh4Q5AhGY7n6ut8QuuHOhwmK8sob/kG16jVwFChIdp2b7jWkG+S7z8DhRnlGEC4tcoYRwOhlUPK9a/deXR+TYfw1OUVAEShsBNCHg8zw7sDYyRWfFOBBS4QYMIfnBzYhyA+8PzqImft6VwIUsA4grCAzWAgMQcYHr/CFY0mAfvU/tcwWuvhAKpg+yBJGBpQABQQ5gtEgV+w+x7dV5syZQxMmTCDpWQQhwkdQMaqEHcogg0qsIygsGS4y9BL001d9v43rrGXHLa+fefFNlbsM14EmpwgoAoWPAHQpBpPYE4R+WurOerc/aD5xgWURIEvo99UuzH29KwEKWAcQVggu2Ds+cgoStHz5cqp9Q0NH0H/9hysNMcKWqvgoFrZF1HmeAQGOaDQoQowc4cNoUJogz/gODUgPK1B4f+BVxHdqtm3bZnYoU7mKqMBksNiQPeizmxs8Stfd9oDxMkriDfmDHF59cyNCB17/nsdoc+k2HbnMYB1oUoqAIhAeBNijDhuQ+2+E0J+wCXVAPL/qWglQEvUBgwGGJ1g8XJnw8kgCdN4FdY0RC2MWcdTNmQS4EY0KYg05wTQkKEgoTngXQXpYgWJaHL7Rgq9MQ+54c42IQqbFzgACPGIJuYNcYXEuPvr6YdfujtxB/h5v/qJZvIvBHnwnB4QJnbyOXmagEjQJRUARCBUC0I0Y+MYAOfffTIAwcK4zN/KrupUAJVkf7AmCkO/YsYNqiZH63114nTEmcE8XuCUJbESjQ54ksQZ5xhozkB5WoJg/DM8PFlnCW6SyFVFhyWCxmQCB0GDOOtYslpSUULfPv3LkDvL39PPtaPHixUYmsXiXybcSoAxWhialCCgCoUBACVBhVaMSoBTqC0IOQwA7ftS6vmqq0n///XrzP9yDUfv/2jsPLymqfd+/f+Gu99Z7b923XGfde9c513XuvZ5zj3gMiIoHBMk5iTkgkgQkSFCyOGSUJAiScchZkCxhggyDEmSGZIABYRhbhAEJv7e+1eya39RM93RX9/R0VX1rrVpV1ZX2/uxf7f37/vauak4kECsBOJSmJwjRoycatrcdUTN+GO8KMfoeK1EeVxUB1FEYkoE/Ai0qKpKCggLrfyqM8MayZ7/3rV5J9E6itwgRTNZtVZHlfhIggSASoADyVqlTALkoLzgARgA99FTZ17r+9Ld61jARCiAXUHmK5Vii9/DMmTMVBBCi7xw/TCNJJgHT+whRA7tDT9C8hcts4Q0B1HfQKOv9NNRpZkgve3+SWQq8FgmQgF8IUAB5qyQpgFyUlxZADz7R3HYY/vz3BhRALnjylDAB2BWGJGG4m7MHiAKIVlJdBGB3ENcYfol/Mtc9QP2GjCn30Q2Kn+oqBV6XBEjA6wQogLxVghRALspLC6C/PdHMdhj+85FnKYBc8OQpYQIUQLSEmiKA3h18YGPh0hV2fQYh1H/oBxIKhayhlzWVNt6XBEiABLxAgALIC6VUlkYKoDIWMa9pAfTXx5vYDsNfajehAIqZIg90EqAAchLhdqoIRBNAePcHDTsnEiABEiCByAQogCKzScc9FEAuSsUIIHwF7oHHGtsCCMPhzp07Z/2xJY7hRALxEKAAiocWj00mAQqgZNLktUiABIJIgALIW6VOAeSivIyjiv/PwLA3M2b+73Vb2gKIY+VdgA34Kcau+A5QwA2hBrJPAVQD0HlLEiABXxGgAPJWcVIAuSgv46hCAP3Hww1tAfRovTbWl5TwXy0UQC7ABvwUY1cUQAE3hBrIPgVQDUDnLUmABHxFgALIW8VJAeSivIyjCgGEL7+ZHqDHG7SjAHLBk6eECRi7ogCiRaSaAAVQqonzfiRAAn4jQAHkrRKlAHJRXsZRhQC6v9YztgB6qlFHCiAXPHlKmICxKwogWkSqCVAApZo470cCJOA3AhRA3ipRCiAX5WUcVQgg/Pmp6QH6R9POFEAuePKUMAFjVxRAtIhUE6AASjVx3o8ESMBvBCiAvFWiFEAuyguOKv45HY7qH//2D1sAPdPiRSkqKhK+A+QCKk8RCiAaQU0RoACqKfK8LwmQgF8IUAB5qyQpgFyUlxZA//bXp20B9Gzrl+XChQsUQC6Y8hShAKIR1BgBCqAaQ88bkwAJ+IQABZC3CpICyEV5GQF05swZ+ZcHnrIFUNN2r1kC6MaNG/wKnAuuQT+FPUBBt4Cayz8FUM2x551JgAT8QYACyFvlSAHkorzgqOLf0U+fPi3/8l9lAqhFpy5y8eJFoQByAZWnsAeINlBjBCiAagw9b0wCJOATAhRA3ipICiAX5aUF0B/+8wm7B6h1564UQC548pQwAfYA0RJqigAFUE2R531JgAT8QoACyFslSQHkory0ALrvz3VsAdTuxe4UQC548pQwAQogWkJNEaAAqinyvC8JkIBfCFAAeaskKYBclJcWQP/v/tq2AOr0Si8KIBc8eUqYAAUQLaGmCFAA1RR53pcESMAvBCiAvFWSFEAuyksLoH/+98dsAfT8G30ogFzw5ClhAhRAtISaIkABVFPkeV8SIAG/EKAA8lZJUgC5KC8jgE6dOiX/90+P2gLopTffkZ9//pkfQXDBlKfwM9i0gZojQAFUc+x5ZxIgAX8QoADyVjlSALkoLxg5vgJ38uRJ+T9/fMQWQK91HyCXLl2Smzdv8jPYLrgG/RT2AAXdAmou/xRANceedyYBEvAHAQogb5UjBZCL8oKRh0IhSwD973972BZAXXoNogBywZOnhAlQANESaooABVBNked9ScC/BO7evWsFg9G2BWFGPQrfEKOD/ukPD9pz3cadrP+ILC0tFfiPQWDhzKOxhXSydgogF6VhBFBhYaH8r395yDbyt/oMpQBywZOnhAmgwvjtt9/k7Nmz8kTD9rZdmcqT/y9FS6kuAhRA1UWW1yWB4BGAs4v2rDgrR755vbsc7T0gEPORt/vL4R59JbdLDxn1x7/Z87RaT0neW2/Lt736BYJDZeV9rM9A+eaNHvLb6TNpM0KKAshF3aQF0P9UAqhnv2GWAIIzgQqAU2wETGTAGTEI2vatW7fsP9h1CqCioiK5fv16IKJHxh5is57EjzL3C5q96fxCXOP9xQVLltvCGxHMfkPGyC+//CKo0/Txfl03toAlJxIgAXcEUD/gVYBvuveRE736ycW5C+XivIVcBpxD1pMN5VJ2rtWWuLOs5J6VkAAyjYVfG8NI+YIzAKegoKCgnLPQe+AISwDBmQhyN2c8JgrGpZcvS/7Lb8rXzdoFes5t2lZymrSRA8+2lBV/fUwW3f+gNWM9q1Erwf4gMMp+urGcnjKt2itJ1F94Ts+vXCNZTzYIBNtI9gPbgo1te7KhbXewv42PPm3ZZFBsD3yy6zeVy/sOMIgVT0XOY0ngHgHUqyaYl997gBTNnCNy5nuRM2e5DDiHQ607yfm9+62AWjoEmVwLIOM8FIz6UPbVqiMHHq8vB+rUD8RyX536suPRp2X13+vIvL8+LNMfqCWzH3hIMmvVlq2PPCV7H68n+wPEw5T7103bSui7EzE7rqaiDJ0/L3v/+zG5umO3XN2+m8uAc7jwyVzL6a7unlSIbwQrjvQfLKeHjKTdBdzuTP2T16ydXNj9lSWO6dWSAAnERwDtOnp/8EGo3O59pGjGbJFTZ0ROn723PMPtgPLIa9VRzmzfKXgXytMCCJFTZOJgm85yacFSuXkwPzBzaW6e/LI/W37Yul2+XbVGcpdmSv7yVXJ681Yp/mq/XM85KDe+PhQYHqbss//RWC59860V/YnFuK3en9JSuXDqlOx76Il7ESJEiTgHmUFow2b5um1nuXbtWsxiOr4mWqzK145S9h8i56ZMp93xubNs4HDHl+SHL7fza57xPlQ8ngQkXLcisHThwgXJfuttKZr+icjJU5zJQPJadpSTW7dZw/lj8RGr+4Fy1QOEhCM6W1JSIrntX5CSNRvC6h4NKFQ+uvmspT+3754+K7+fKJRfjxyVi7kH5aesbCnKzpWSw9/Ije8K5I5R9wHhYco7t0Fz+TEn1xLGEDdVTRDR1kv/R4+GBRC4GXZmndtlkTPDxETSfLodWr9Zslt1sr6mA5FSHRPqMEQpL1++LLm9B8i5ydPCtkd7C5y92ZHpe88TBNCpTV/EXI9Vh33ymiTgVQKoWyGA8N7qgTd7StG0T0QKTpbNhafK1vE7twPDI69lBync8qX3BRB6fy5evCjZbTtLyer195yH02VOhNWY+Hf7TuEp+f1EgVw/elyuHTlmLW8e/05uF56Su8YxPeXf/FtCxZG/nGeaScHuPXL16tWYho9AAOHYk4cPy95adRghYoTIsoHQ2g2S06qTFWCpTgGERtqqw3r1k3OTPqb90f4sGzjc4UUpWL/RaqRjCeR41VFlukmgOgiUE0BdekrRxzNFThRwJgPxjQDCF6nOnz8vWW2ek5JVa8Mq/uTpe8tTgdi+W3hK7pwolDsFheElxA+iGXAkrGWweEAAHd+x0/qSGcRNVZMRQAWHDoUFUKGKElmRIW6Xj5wFg0dozQbJbtUxJQIIwzSyer4j5yZ+JEL7c0Qig2Fv9jN2r/zz278gJ9ZtoACqqgLnfhKohEB5AdRDij6aIXL8hMh3J8JLrHM7kDzyWvikBwgC6KeffpIDrTtJyYrV9xrOQi6trt5gcoAAOrptu3sBxAgJo2QnCiS0ep3ktOooV65csd4nq6SNTfgn00iXCaCpZM/nz7KB/HYUQAk/YLxAYAmYutUaAtelhxRNnS5y7DhnMpC8Fu39MQTOFkAYqgIB9N29Lk5rWcjtAPLIqd9Ujny5zaUAerwsOmRFir7jdrlIWXB4hFattXqAUiqAxk+hvQXU3sKR6bLniwIosL47M54EAhUE0JRpIkeOls1Hj5Wt43duB4aHDwVQRynJXClyHA0I5yAzSFgAMULCKNmx4xJauSb1PUDjJ5M9nz/LBvLbPc8hcElwhHmJYBKoIIAmfyzy7RGRb74NL7HO7UDyyGvuux6gjlLy+QqRo+ji/C68xDq3A8cjYQHESFBgIkFWRDBCeYeWr0p9D9C4SYxERigPO3obkP35bSmAgum6M9fJIFBBAE36SOTwN2GHH0szQxCZde4PBJ+85u38NgSuo5QsW17ecdPdnVwPDJuc+k3cD4F78PHy0SETJeIycFxCmStT3wOUMTFwnK0oLJ+vCuWe37Yze4CS4QnzGoEkUEEATZgqkn9Y5FC+L5fXDmTL8qmzpEPrV+Rf/+MJ+ac/PGgtR/d7X27kHvRtvt2Upz8F0NLMsHqFojfdnGad2+EGNgA8EhZAJhoEVmadkaFARIZ0eUMAZbdM5UcQ+sq5DyeU2Rztr4xFAJ8/CqBA+u3MdJIIVCqA8g6JmBlCyKxj6eHtkxu/kFbNX7BET8bAkXJtf5aM7PuetT3vg8nhfHo4f1Y5JTH9ec382AO0+GE2tlcAACAASURBVPPyDaZ2XrkeGDY59RLsAUKUyOfztayccLSojTNaNCwcLfJ5/mMp39Cy5ZLTskMKvwLXV86NneB72wP7H7/cIeMGjZL6z7SzGmlEKx96rLFsmbMwEPmPxf7y2zzHHqAkOcO8TPAIVBBA4yaLfH3Qd/OlHbvk+XavWfVorUeelSMr11p5zJw43fpt34JlvstzouWY16ytD4fALVpavpsPihGz6fb0+PbdvENyYv1mGdxzkDxZt6XtONR9urUcXrkunFcf5dcqOxf5SVgA4Z4+iQxVFjk5uXFL+WjRgezy0SKf5z/WSF9oyeeSnWoB9ME4T0ciK7M3zfvuwTxZN32uPFDrGbn/r0/LrnmL5dy2ndKsyXNS65FGcmTVel/nvyo+ej8FUPCcduY4eQQqCqBJIrlfh8UAlmaGKDLrHty/7uPZti/YvtVLUrLrK0/nxxI31Vwe/hRAC5aEHdeDeb5b3so9KJ+MHC/33V9b6jzZQo6uXi9HVq2Th2s3kaaNO0nRtp0iPsy35RDEma+ceo0TeAeodpgj7unD+dLO3eWjRavWWfnMnDQjHC1auMyX+XZTlhBAKe8BGjPO1/yzFmda4gc9Pm++0F2u7t0vpQeypffrfaRF084C+3RTVn48hwIoec4wrxQ8AhUEEN6vzMn11Vy6d7/0fu1tWwCNeHuw3M7K9lUeq6PM/CeAWnaQEggg3cUJB9Yn2ysnz7TEDxyHEb2HyO2cXCnetUdaN39Ber7aS67vOxB2HHySX6vcXJZfwgLIRB88HhmqLJJSMVq0pyxa5MP8JhLZCy1elvoeoNEZvi2PGweyZWDXfnZjPXHQqLK8ejDyWtnzlYi9Oa+X37oTh8AFz29njpNEoIIAwvuVB7JEIBDM7PHtH9Zvkgb129p16ryR48vyhjx6PH/Vlf68pm18NgQOAuizRSI5X4dnNKhmHUsPb5fs3GNH7SGAVk2e6en8WOVSjeWR848Ee4B8FiUyERRGi+KL/oUWLkl9D9CoD30bvStcs0HqPNncbqx3zlng27yaZy6RJQVQkjxhXiaQBMoLoO5SNHa8yP4DvplLd+2R6UPH2PUpfEMz1368qRSsWOObvCa73HwqgBaKZOeEG1UfLbMWLLN7f6yX3Jav8WU+LWchCeWWsABCGhA98dnyhw1flI8WjZrgy3wmq9xCCxenvgdo5Fjf2Z0pj03T5tgNtFWPZa6m/UWpZyiAAum3M9NJIlBBAOH9yr37RPbt9/wy88Mpdl1qRI9etm7aWYq3fOn5fFZXeeU18WMP0Nz54S4/dPv5aNYqv+mzHaRo0xZf5S/ZZZWQAPpb7XDUBPaDaJFPlqW7v4oeLcKXY3yU32SUW+izhantAerRV86N/MCX5XB+4xfyasc3Km20e77UQ67v2uPLfCdih/mtOnIIXJKcYV4meAQqCKAxGSJf7fXNfHv3HhnRY6Bdp77Y9lUJQfT4KI/VlRf/CiDdxWmcV/Obx7bv7j9gdWM2alD2qVit8sf2HSp3Ec3waP5s5yDJ6c95upH7jyBAAIGp4WrWPbydmTHVriS1/Zh1O1qEcvBBfu08JJgfCKCUfgUOAmj4mLIySDD9Vr1Qw+VZunO39H6lZ1T7y3jnPV89b8myPwqg4DntzHHyCFQQQKM/FNm9R2TPV75YlmzeIu2bP2/XrSO6D5DbO3f5Jn/VWU55TVr78B2gOfPCXX7o5vT4fHnzVmnTrLNt3MZZ1Ut0g3o9n9WR/oQFECIosB8dSfH49u09X0WPFnk8f9VRXqG581PfAzRsdNjufFYeRes3SdOG7e36bPrgUb56vqrD/vJbdmAPUJL84d9+uyYlv4Tk1q3bcV3xzp071nk4F+vxTLgXzsO9OaWeQAUBNGqsyK7dvpkLli6X2rWb2HXq0tETfJO36i4nfwqg2XPLGlWofO3AenT7yOJMqfVwQ9vIN02ZGc6XR/Njl0k1pj8pAgjpMzPsyKybdHtsu9JoERoDj+YnFeUR+vQzyW6Rwj9CRQ/QsHvCwGP2VVV55M5dKPf9+2N2PbZz+pyyZ8oHz1dV+XeznwIocad5554sqd2gTHj/61/qythJs+TatetRLw7xMm/xSnngsTIHE+ufr9oot29HF0K49ujx0+W+P9ex7f3pZs9Lbt43Ue/JncklUE4AvdFdijC8eMdOEfSS+GC58+NPbPtC3Zo7Z74v8pWK8slr3MqHPUCz5oS7NtHN6ZN5VUbZy24QQkcWLfNN3qqrjBITQI+Vj6LAjnTUyKPbEaNFHs2PXSbVmP7QnLmS06K9XLlyRW7dupXc1vne1UwjfeHCBcmCAHp/pC/szVk+iE6a3mtELWGP1jHVWH5ev35+i/bsAUrgqbtcXCJtXuhh252xPywhjKJNx0+ckof/0abCuXUadpTCU2ejnSqbtu6ucB7u+WLX/hL69WrUc7kzeQRM3VpUVCQHIIBGjBHZvsM388Q+Q2w7a9qgnRStXOObvFV3OflPALVoLyUz54TVPRxWqHyPL+/u3CVje71rG3mLRh3l0roNns9XdZdLTt1nE3gH6LGw7Wj7wbrHtyuNFpnnwwf5q47ygQDKTrUAem+E7+zv9vYdMqJbf7sea9+ss5Rs2FxWj9H+Kq1fKIASc4YhVCBYtPAx69PnLI568UgiBudn5eZHPTdjyuxK7xmLeIp6Ye6Mi0AFATR8tMiX28rmbdvL1vG7h7avb9wsPZ/vattZz85d5fqmL6LmJ2vmp/bx5jmobNn7hbekdONmT/GwyjWO8str5LceIAigGZ+EuwDRvemD+eqmL+RN9eWkd17uITfwlQ8f5K0685CwAEKUCIx1tMjj2xP7DrUrPytatGqtr/JXHeUVmjUn9T1AQ4aHy8Xj9qbLo2TdRoHoMY3tiLf6yW00Vj56vnR+rXwlofwogOLydyscjCFnehiasT8sIVKiTZmrN9n2qs/DOgVQNHLps6+CAELv+patvpiLPl8hTZ8p+wPUiW8PqjJfmSMypHWjjlK8ao11bNa0sFDH7+CC37E/o+fAKq/ldY55jVr6bAgcBND0WWHViobVNLAeXp5amil1nyj748B5Q8fEnC+o/TqPN5VCDJnzCY9Y85GQAPrvx8p4wXa0/Xh0+/rmLeWjRc+/FY4WVZGfDNX7aJwA2JVdDlWcb7HzML/QzNmp7wEaMsx3fAsWLZPatRvbDuXSkePCz1UU+yles15aNy4fva/zeJPy9VmU823b87D95bdoxyFwCfjTBSfPlHv/x9RhWFbVA7Rq/RbbXvV5WKcASqBQUnhqBQGE3vXNX/hizp0+u9w7lZs+nFRlvjK6D5DM4R/ax2G7Tu0mUvjZIuu30nXrpXfnrpL18Sz7GL/wcubDnwLo4xnhLsCtX4qYGV2bZh1LD23vUy+5oeLFdlXpL1y4VOAo4HjLYVi41LP5t8rNRXnl1G3ofggcBJC2Fx+sF2WurBgtipKv4lVrrUiQZT8Lllg8zG+WA4AX2KOc75d9oRmfpL4HaPD7vmO7c/J025m0XtaFiI5iP8bWwtHKteXsT9tktGv4YV9+cwqgRPzlogs/S9P2lf/3FIa4RZsOfXNM7n+ogW23qPesNjWGd4Ai9R7xHaBoxJO/r7wA6iZFQ4eLbNzki3npe2Ns26z9aCMp+HR+XPkqXbNWend+U1o/216KM1fEda4fGOY928KHPUAfTRf5Yku4+84Hy3mDR9pGXveJZnIKYiZCvopXrrac1t7Pd5XSDRslo8c9dT9/sW94WN2uEfKvueQ8laAAwj30fTy+nTt9jiNaNDlq/sK21EGsrnGT9y+2iOkyh235iY+Vl0rKOzR9Vop7gPrIuXffC7OtJD22/ZsywXAOs57Gx09UPYkNnm4pPyz+PKr9GPuDven8wR4tAT7tk6jn2+d4hE+k9FIAJeYUwwFeuW6L4MtvRsBg+e7w8VV+Ba609IaMmTCj3BA6DKebOnN+lZ/S/vlSsbz01oBy93yobivr4wiJ5Yhnx0OgggAaPExk3XqR9Rs8vby9dp2MeP1t277aN+4oJcsy48pX4Zx5UuexxpLx1jtxnecHfih/fwqgKR+LbNrsi/nG2vXyzgtv2Ub+YssXJLRydcx5y+jeX+rUbiyF8xbGfI5f2CUsgHxiQ6Y8y0WLHmskBXMXuLKJrI9mWvaYOWysq/NNeryyDH08Q3Kap/IrcH3k3MChvmJ7fc066dmpi12PYR2/xWsDpWvX3YtYdpDi5SvjPj/e+6XD8fnN2nIIXDweb4RjIWaKr/xizfgKGxzjWCf8h485N57/8zH/H2TOvfn777HekscliUAFATToPZG16zw/lyzNlPaNOth16qCXusnNVavjylfWpGnhtnzo6LjO8wM/5CGvoe96gNpJyaSpIhs2hrvzsDQzuj3Nukf2X1qWKS0atLONfGzXd+QuIhcxpj+jWz9L4RearlGP5T+R8kpMAD3qC/sx/G6vWy8j3uht25EVLYID6cIeYFOtG97rMndxvn1PjzyPoanTJLt5uxR+BruPnBsw2Ff2V7R4mTStX/Y54YnoPYyz/DFEA3Zn2d7ny+M+P977pcvx+U0pgJLkC/MyASRQQQC9O1Rk9RqRNWvDS6x7cLtg1qdS+5Fn7TY9c8iouPOTOSQ8uihr4kee5+GmPPMaNvfbELh2UjJxSrhrE92cHp+PzJgjtWo9Yxv5qmEfxJUndG2iixNdnV5nEW/6c55skMA7QI/6ihe6xts3KnuZfNBL3eUmKv04n4+sydMCZ0+hKR9JTqoFUP9BcZdNvGWZyuNzp86Q+/5U9geom0aPjzl/sDk9dMkarhGn3aYyr8m+V37TNuwBCqDjziwnh0A5AfR6NykaOERk5SrPzzs/KPtPtYf/3kCOT5sVV55Kl2VK7w6vS+sGbaV4weK4zvUDP+TBfwKoeTspmTA5rGZ1NyfUvge3V73/gd343/9fT8ohfJkD+YgxPxldIYAaSeEncz2Zf7vMYsyvPj4hAfTXR8OMzX2xNLPh76FtK1r0qCNaFGf6YUOtGrSTwtnzyljgGh7kYZdlDOkPTf4o9T1A/d4tY+wDvkvVe4x2fRSn/aHMsiZ+bNWHrRu0k+LFy8KMfMAnmj3mN6EASo4rzKsEkUAFAYTe9eUryuYVK8vW8bsHtu9mrpBxXcpGdAx8rovcWPp5OB8xpr94/kJp/Uxb6d3hNSldsqyMQYzn2ww9fHxeg2Y+6wGCABo3qXxXoOni9Njy7qo1MlYZeZuG7eTyoiVx5S3jzb5hATTr07jOs7qEPcbLmeaEBRDG0/pkrhAtmv5JXHkrnDlHWj3TVrD0C5NY8xGaMDn1PUB9B/qG8+0Vq2TEqz3tQE63tq/INbys6/LZyrwnprB0ew0vnZffpDV7gILouTPPSSFQQQAhuPR5pqfnC7PnSot6ra069eGHnpEjk6fFnZ+sDyda52cOHBb3uV7nZ9LvTwGUMSGs4tHNCXXq0eXVJcvkzVYv2o7DmNd7yW18qjCO/GR06S11Hm0khXB44zjPy9xMPnOefMb9EDj0AMF2fDDfXb5SxnXpY9vRwM5vyg04oDHmDbbT6pk2YRuK8ZxYr+2F40LjJ0l2sxS+A9S9j5zrMyDm8kl3hpfnL5Q2Dcr+rG/d+2MSylvWuCmWLaNuS/e8JyN9FEBJ8YN5kYASqCCA+g4QWbJUZOmy8BLrHtq+sWCRzOjWz6oDH/hLXfny/TFyd3H8+cl4rZd1jcz+73kq/1ZZJam88p5p6sMeoLHjwt15mcs9vSycNkvq3Bu2hOFvBydMjTs/GW9AAD0rhfhzWI/ziDf9OU8kKIDAywfzhTnzykeLpkyLOV/Fny2QV5t3FtiiZpGVMUl6t39VSlHx+oBRtDyEMiZITqoFUO9+vuG6a/Q4ue9Pj1qNbddWL8qvCxbFlDfYHGwPNqjLB7aHd4KsyKXPbQ/5zm/cij1AAXXeme3ECVQQQH36iyxa7Mn5yPjJUutv9eShB+vJ4E5vyPcfzXCVj4xXe1h1qHm3MrPfEFfX8SpHk25/CqAPMkSWfR7u1vPwcmHfIbaRjny5m9xCxCKO/BTPmy+t7315KWvsRM/zsLot48h/YgLokTBr3A+zYe+x7RuLlsiM7v0tO3rgr3Xly+EfyF2MFY4xPxmvh6NEpqLUS+yz2HiYTyzpD304PvU9QL3eKWPrYb4XZ8+VV5t1tuyvfYN28sOM2TE/T4UfzZA6jzwrlp3ds9fiuZ9ZdRrqNazHUn5efn6RdgqgxJ1gXiG4BMoLoLek6O1+IvMXiCxYWDZzu4wFuASER179Jn7rAWorJaPHlu/aNF2cHlqenDpdGj3Z3HIcerV5Wa7gK24xph9dmtpR1eu9270ipfMXxnytWO+ZjsflPFE/gSFwj3ie0ZGJU6XWg/XkoVr1ZfBzXeT7aTPjylPhlGlS55GGEW0p47WecV0vHW0kljSFPsiQnGZtU/cZbAyB69nX82zn9Rpo2U6Tp1rIpz0HyDX05MRYh5njYGO6/sJ6UOzOMMhv1JI9QMH135nzBAlUEEC9+op8Np8zGYj/BFCztlIycky4O2/xkvLdeh7Yvr1gkYx4oas1ZKRN/day7t1hcgOOA7psPZB+07WYDunNqZOAAPrLI2HmCxeVtyFuB45HaMxYyU61AOrRh/aHOo/Pm1AAJegB8/RAE6gggHr2Efl0rsjceVwGnENePb/1AEEADR8V7sJDNx7nwDLIqVPPfQ8QBBBthwzmL5DQyDGp7wHq9jbZ8/mzbCD/2RbsAQq0C8/MJ0KgggDq/rbInE9FZs/hMuAc8uo19tkQOAig90eIzPuMc8AZ5DyeoAACP0SJuAw0h9CI0anvAXqrJ+2Oz5313OU3bE4BlIgHzHMDTaCcAHrhNTnepJWcf7NH2dxVreN3bpex8TmPAw/W9qEAem94uGsT3XucA8sgYQFE2wms7eh6IzRsZIp7gHrLua49yJ7Pn2UDFECB9t+Z+QQJGAF08eJFObh6jezo+KJsfe4l+aLjC7K5w/OcA8YA5b6l04uy48XXZN/LXeTHo8ektLRUYCc1Pf0PNwlAwq9fvy4//fSTHEAP0JD3RT6ZHe7i5DKwHBITQA8Hlps1NIDPjV3+ofeGS3bTNin8CEJvOdelu31/lsecQNfn+Q2asQfIjWPAc0hAxHJsf//9d6v+PnnypOTk5MiuXbtk+/btnAPKYOfOnZKVlSUnTpyQy5cvy82bN70vgM6dOyf7nm4kh55+Vo42b8s5wAyOtWgnu//4Fzm2fYf8+uuvcvv27SobAxxz9epVOXn4sOy5/78F16AdBfs5gg0crtdY9j/VUEpKSuTWrVtV2pGbA3SU8sAb3eXrvz9B2wtw/WXqHdjf3j8/KAXrN1pBvjt37rgxL55DAoEmgLYdQfJLly7J999/L4WFhZbzCweYc/AYFBQUyNmzZ+Xnn3+Wa9euxeQfpuIBct0DhC6soqIiyVm0RLaM+kDWDnlfVg8aKqveHcI5QAxQ5msGvycbh42U3eMnyenTp2MWQHAufvvtN/nhhx9k3+SPLDva8P4IWTd0GOcAMkDZfzFyjOyZMFmOrFgloVCo2ipKCCBEodBAf3vggOz4cLxsGj5K1r83nLYXQNtDucP+tn2QIdkfz7DaNrRxFECpcEN4D78RQP2K4BWeIQQ5UZcjoMU5mAx++eUXyy+EKIZdwD7SYXItgIzz8N1338n+/ftl27ZtsnXrVtmyZQvngDH48ssvZffu3ZKfn28Ni4TCj8VxwEOAChJjhWFHubm5sm/fPtm7dy/ngDLIzs6Wo0ePCnqXY7UjNxWpaaDRMCMydfjwYTlw4ADtLqB2hzoH7djBgwfl1KlT1vCddBmm4ca+eQ4J1DQB1LHwA9AbBKcXw+I4B5MByh8z7CFdxA+eD9cCyAxfunDhgtVgHDt2TI4cOcI5gAzgsKJb+8cff7QiPLE6DsYJRS8QukZxPpzRM2fOcA4oA/QGQhBjGCUay+qsLFEZ37hxw7LZ8+fPW0M1aHvBffZQ90B4FxcXW8N30MZxIgESSIwA6nDOZFCdbblbC3UlgHAzZAYOCqK0iKJeuXLFajjQeHAOHgN0baOrGw5lLL0/xmBhR4gMoCcIQgjXgPPLOXgMUPawAdgC6pZ47MjYUzxL2B6cXAh2dM3T/oJnc7qegf2hPUMdBrtIxwY7HvvmsSRAAiRAApEJJCSA4KCgoTDdm3AkOAePgenWhi3AJuJ1HHC8sSVcg3OwGRgbiteOIldzkffgHrS/YNubs74x9hfZariHBEiABEjA6wRcCyCTceNAcMkuTtgAJxIgARIgARIgARIgARJIZwIJC6B0zhzTRgIkQAIkQAIkQAIkQAIkQAKaAAWQpsF1EiABEiABEiABEiABEiABXxOgAPJ18TJzJEACJEACJEACJEACJEACmgAFkKbBdRIgARIgARIgARIgARIgAV8ToADydfEycyRAAiRAAiRAAiRAAiRAApoABZCmwXUSIAESIAESIAESIAESIAFfE6AA8nXxMnMkQAIkQAIkQAIkQAIkQAKaAAWQpsF1EiABEiABEiABEiABEiABXxOgAPJ18TJzJEACJEACJEACJEACJEACmgAFkKbBdRIgARIgARIgARIgARIgAV8ToADydfEycyRAAiRAAiRAAiRAAiRAApoABZCmwXUSIAESIAESIAESIAESIAFfE6AA8nXxMnMkQAIkQAIkQAIkQAIkQAKaAAWQpsF1EiABEiABEiABEiABEiABXxOgAPJ18TJzJEACJEACJEACJEACJEACmgAFkKbBdRIgARIgARIgARIgARIgAV8ToADydfEycyRAAiRAAiRAAiRAAiRAApoABZCmwXUSIAESIAESIAESIAESIAFfE6AA8nXxMnMkQAIkQAIkQAIkQAIkQAKaAAWQpsF1EiABEiABEiABEiABEiABXxOgAPJ18TJzJEACJEACJEACJEACJEACmkDaCKCSX0LSre8w+ac/PCgtOnWVU2d+0OnkOgmQAAmQAAmQAAmQAAlYBL7YtkceeKyJ/Otf6sqMT5fI7dt3SIYEYiZQLQIoN+8bue/PdSwxA0FT1Txx2lwpPHVW6jTsaB+blZsfcyZq8sBr167L8jWbpc0LPew812/xkoybOkd+PFdUk0lzde8bN29a+XmobiurLNq/3EsgTjmRAAmUJ8CgTXke3CIBEiABNwRQl8LXqMpXNPuNX5IxZbZ9Tu93R0tp6Q03t0/pOXfv3pUThadl8MiJYvwsiLiXuw2Ubbv2CXwwL0xg3f+9D23+pmyiLdPNr08bAaSdidbPd5cffjqf9jbw86Viee61PhEN4MWu/SX061U7Hz+dv2D1bsFA3uw9VH69+pu9r6ZX7ty5I3v25wrYawPGdvGVX2o6ea7vzwiRa3SBOpFBG+8FbRB8Wrlui+U4wIFAvYVI8EtvDbDqMtRpnEiABNKfgFsBpNv3TxcuF4iLdJ6QPgTMo3UQbNq6u1wWkC8cj7pt64695fbV5MaZ73+Ses1fLOcvat/Rud7iua5y4eKlmkxyhXtXiwBy3gUOtHas000FOtMbyza6WsdMmGEX/thJswQNMgz84s+XZda8ZfJWn/fLCSDdy5VuwuKTzz6386INN13SeePGTXlnyNhK06jTa9brNulsDaNkhCgWa66eY/AswOanzJhvCX9T6T/ZqJNMmj5PLhVfqZ4bu7iqGwHEoI0L0Ek6BeyjBZ9QD0ye/pncunU7SXfkZUiABFJJQLfdWPfDdPzEKXn4H20sPwbiIS//qOUz3vz9dzmQc0he6zFInAJIc8hcvSltMThFrDMf6ZhwCiCXpVJ04Wdp2v4Ny5AbtHo5ph6rdBZAeMgatX3Nik4gymCEhNcFECNELg08CactWbE+aqQLDQCGAqTjxKBNuFTStc5C+bzS7V1ZuGyNFXCC2EYACoEoU3fB0YDDwYkESMB7BLTj7xcBNG/xSrt+mj5ncUyFojmkswDSbUWtp1rKkeMFMeWvJg9KGwGE4WAYFobGCx9BwHAxM23Z/pXlSOHdGnwc4dtjJ6T7O8OtLkF0C/YaOEq+//GcdTiW7w4fX+k+cz2zxPA09Hw0afe6dd94hk/owoYQgiCKNFXVezFoxARBBMBMuLYeH4pxotg2eTTHYanZnD77o3UMrochIYi4YzhI/rfHq+waxrhTM2QEPXTGiUgXAaTzbNZRGZh0emX8r0l7pKWfIkR4KRVjhI8eL7Qi8bAxOKymJwhl57T9SFxS/bsfBZDfgjaV2YTOI+zLC1HIyvLB30gg6AS04x9JACFYC78NbQqGipkJ71/Dr8Pvazdus0YbYNSBeecGPuaur7ItnweBE5xbu0F72/80+8z1zBLvvaxY+4X9zjeuj/e/123aHtO7OzpPEEPRpmijEu5/qIEc+uaYfXo86YKvjeAR6keMzkC7jPSDCX4DI7Td4BLPtC/7oO2PYbhbOo3wiJSPtBFA2uHAxxAgAsxkHF2oSjhMMHgUlJ4btn7FEjNmLLhzn1M84Po4Rx+n16saPuFsaDE8K1KBxyqAEMXEmPbK8oe0IW97sw4aLNbSsMHDO3XmfPsB13mp7LxyF3FsUAA5gKRw088RImBEN/nzb7xjP3fpKrB1fYRnqbJhuwza1GzQprLH0lluO/dkVXYYfyMBEkhzAlosRBJAxv9BHa2P0QHqvoM/qNTXg3hBL0xlH1/Avo1bd5UjhHe+EVDWvpVehw9YlWjQ7TuEBgLYJvBc7mYiEqsAijdduo7sM2hMxFcLYsmPTvPSFRtsNjgXfm+6T54SQMbYho/9yHKkIDjMp7PNPvQG4ferv10TGL75XXc3wnHp2uc9ax/EgTFCOGdDRk2yfscDsGtvdsTywztA+HqduT6WEC4fTv7E6oWBmHFO+qGszPFDl6EZH4qHEi+Z4TrHvjspzTp0Ze11XQAADhNJREFUse7lfJFMVwBIgzkPqv4z1d0KpxP5i2XykwCKFCFyRkFMFMi81IclIj3gaCIkGCIIxnqfkyfK6+tD30qP/iMswYrj43nnRVf6fosQgRUiVeitAxfMlT0HTqY1sa0bCaSzMgGkj2HQpmJQCtwqC76YOiuZQRtjI7oHlUPgDBUuScB7BHRbqMWNzompS1DX6GO0r4V9rTp3k+8KTlltufaLsA/tOd7FgU/3+aqN9ggF/RErvEs4MmOa1WbBN8RIBvgFaM/gW+I6mPF7tMmZLpyDNvDLnfsifsFOc0B+9eQmXbrdwv11fjDCCukx+Vm3ebu+XcR1sBvx4cf2edrfjnhSGuzwnACCwNEqG0McTGHB+ccHCMwEpwWFi/1akerf8SEDFJ6ZtIGiQPU+c4xZXikJWcPvzP3NEvccNnaqYL+e9LWdjh8cZ3w6G9fA+U7xpR90HdXUv+NB1u9U6Gg7romIQiyTnwSQ5qMrSF0JvNr93XK9EqYcsRw1blqlZYx96CbWQheVEXrhjM3p62DdWT6VlYWfI0TIr1MA6UamMh419Zu2D5SdWwFkbIBBm7Kv/+hnEnySEbRBFBX1q+5drKoXv6Zsi/clARKomoB2/HXbrc/UdYk+RvtaCE5pvwhfGMZ726h7MJRs974c+5KXi0usIW3YZz6khJ36egi6ax9Un4PXOBB8jzTBX4DYQWDItA1m+XSz5620aJ8C19EckF89uUmXs21z1pMQPSZN2m/W93Wu43UStOXmPO2jOo9Np23PCSBt5ACpnXXneyDaOPQ+rdhXrd9Srjy0cVRlzDjR2TtgDABLfKUI3ZNm0ulxCiBtQPrBM+fqfOqegUgVAM5zqnLnw2Ou7VzqeznT6Ty2Jrd13nX56jTpY7Tt6HJGWaGS3J+dZ72rAlGNitGUpeklBE9Ulqbycn78YseeA7b4MT2RcMwwbNH07DkrT51WrGsbMfdHGXg9QmTyefaHc+U+nekMQJjjanrptI9EBBCDNg+KbhD1M+kMCsQTtIET4hwBgGcG9Sd69fkFuJp+inh/EnBPQDv+uu3WV9R1iT5Gt6NOH0bX7c6eex2g0/t0oB0jf/Skz6nqfXBznn5X3bTzWCJ4is9kaxGkOSC/enKTrmj5x7UxEgmvmyA9Tnb63npdi0qvfAAB6Q+kANIGpY3PuR5r4QMkGls40B1eedt2nHE9fGTBTLE+lM50OLf1gx6pAjD31Hl1PjzmGOcyaALIGQXSYhTs0S1uKiRd2eleNf27c+gNzjW9e1VVDjjWjxEi2BieEcMBXMEpXb8UoxsJpDURAaSfV3DQz5dTuOs6Qu9j0MZZS4kVha1MAKG8ML4ejgQCVJxIgAS8R0D7Ls461OQmkv+j61GnH6frdi1ycE3djut9+j6oXyLN+hyTxmjLny8XW6NG9Hvf+HgDPuJgJs3B6cO5SVe0/OOe0diZNDmX+n0lr3wAAXmgAIpizPiCFR6IeCY0uPrdIO3ERDMsbZSRHi78DqcbH0owk34AKqskoj085hrOpXbQnJWH89ia3NZ515x1mvQxmo/m7ay0dCUI5k7nVzM1+3QEBEN6nO9b6XToaLhOq173W4QIeUNPmOk9A9f5S1bZwlLnPR3WtX1UZgNIoz7GaUO6vLXd4Tz9fDntVtcRep+2OaQn0hzP8+qnoA244pmbu2hFuY/IOIeppoNtMQ0kQAJVE9B1nrMONWdHqmd1PeqsE6PV27rt13W6vk+kuhe/YwhupI9hmTRXtsR73ugNN9c2fgWO1RyQDj25SVe0/OPa0djpe+t1jKQyaY912Jw+v6bWAymAtEBBF2Kyp0gOTjTDQuMNxxlG5BxaFS19+gFwVhL4tDa+mmcMM9a86vQ7K49oaUn1Pp137SzqdOhjNJ9olYCuBMFOV0a4tq6QzD5dtoZ3pCXSFOvklwiRU/ygktTjqGPlkarjtH1UZgNIhz5GN5bYF8nusE8/X0671Xak92mbi2RX+D2oQRtjF+hBnTZ7kV3n4RO1GKPPiQRIwFsEdJ2n226di0j1rK5HnT5MtHpbt/26TtcOfnW84K/v62xvNAen7+AmXdHyD7b4vLZ5BSCWd3RR5+r/X6sOPrrMk7keSAGkx02isUQBJnPSL5Hpdxx0L4HzodRiBcZ3MP9ITEnSFQA+660/933h4iVBdyQeqKqGXumbaQfNmU59XE2v67xrZ1GnSx+jK9FolUC0ygjX1hVSvAIIPSDmHJ3Oqta9HCE6X3TR/oohbNH5blxVea+J/do+nA2SSY8+RjeW2B/J7rBPP19Ou9UNt97HoI2hXvVS83WWS9Vn8wgSIIF0IKDbWd1267RFqmd1Per0YaLV27rt13WHFgUDh41L+tBa/e6jcwi9Hv6M/OrJTbp0/tG2OXvJ8SU7/I65qg+BIS3Xr5dKzwEj7XNiGeGi81CT64EUQFoYwCHFOx4wfExYHsg5JK/1GCTODyTogsKYR0QXIaZgUBBReOEdv5uuTFw77/BR+zT9bgkECfZhGIqJhOPLb+YLYvhsI/bjmrg27oEx7eglwqcKzaQrABgsvlOP8aPIx4SPP7WNUjtT5txIS+1AOCuPSOfUxO8675Hyp4/RlaiuBHRFh3zoShBMnYJFV8xmH8oEL1/j+FiiJvHyipYmnR5nBVnTESLYNnp7wAWzF8QPykbbR2U24DzGaUOR7A7n6efLabe64db7GLSJ/YnRX1KsbDhq7FfikSRAAjVFQLdruu3W6YlUz+p61OnD6LrdWW/rdlbv0x9cgY/20awFgr/SwITgNd5lRe87vqgWacL/4gx4P8P6qxQENHEeJlxH96Dgr1j0u4u67h86erK1D+kxvqN5DzLWdOn8o23D+0fwgX+/dcvKh/FfEYjHn5tWNen/xMRfGxScPFPVKWmzP5ACCPRRsOZfgY1zppcwJuenqHWpaSdGn2fWYVR4V0f3LmEdatscY5bG0YFBR/uMMo6P9j9AeGD1Oxbm+s4vLel8mHUdYTbnOZfV4dib+7tZ6srPMHReRx+jK1FdCeiKDufrShAMjMgx19YVs9mnxS1442tnyZy8GCGCPaNBMHYUix0mk1ki19L2UZkN4Nr6GKcNRbI7nKfrDqfd6oZb72PQpnxpwilAkOqrA1/bASQ4FNt27bO+6GhsLp3fMyufI26RAAloArqd1W23PiZSPavr0WQIINwT18QoG1O3VLaM9j9ATr+isvMRwNZfDsZ9Maqnsvsa3yPedOl2CyIH/1NYWVqcn8fW3PW67oXyWsApJQIICheflAZkOAr4QyrnpI9p0amr/HT+gn1IpD+zxAHfHD1hCxn8Can+357TZ3+U+i1esu777vDxFf6ZFj0l+CqVOQbpgzGMHj/d+vNR9L5EmqDQ0dUH9a0NCNfCNfVXPPQ18BDMnLvUTjNEGL4UZ9KNe+YcPGx1KRqBBjGG3ib8gRfeCdGTswI4erzQ+pNXnAMxhHeA9LA4fa5epwA6a+NwVlSmojEH6IpZ79Ndx3DOzB/ZokwRJZk1b5m80u3dCh9IMNf1U4QIYh8BAPN1G9gi3gPyyqQbCdQLupxNHvQx1S2AcE8GbQz5sIhEHVdZw43fsA/1uOldLzuTayRAAl4g8OnC5dZzjGd5yYr1lSY5km8I/xF+JOoC+J6mtwYX0b4m/CotONAGw1fEefDl4EPqCXU+/DVc29Q/8NMwLA5tBIJ+0Sb4hVNmzLc+L63bxpe7DbSCN7rnR18Hfi7EEe6JGV8b1iOB4kmXs93CtfHxGPTeIN/I27pN28v1Qum0ONfRLiEvSBde+TA9W87j0nE7JQIoHTPulzQ5BZBf8hVLPnTedbRcn6uP0VEkZyWAKIqZ3AogOFum8kRFUtkMwRzJKXPet7LzvRIh0j1ileVD/5aOL03qRpJBm/QL2kBg488NIXKcAazBIydWGcAyzzqXJEACJBAkAtF8nyBxQF4pgDxe4pEcfI9nK6bkm+gPnGln75+5gDkG0QlElMykHVxnjyMiGCMzplkCpjLnFxEc3BNRD+fHKhDBwTAcRHTMcEQcB+GCHhEdiTJp0Uu/RIgg5jAmWgudSOvofeQUHAJBrrOCU8rMKQmQQDoSoAAqKxUKoDIWnlyjM+HJYkvrRLOCTOvi8XziWGd5vgiZARIgAY8SYPteVnAUQGUsPLkWqYfDk5lhotOCACvItCgG3yaCAsi3RcuMkQAJpDmBaKNf0jzpSU8eBVDSkfKCJOBtAhRA3i6/dE89gzbpXkJMHwmQAAn4nwAFkP/LmDkkgbgIMEIUFy4eTAIkQAIkQAIk4DECFEAeKzAmlwRIgARIgARIgARIgARIwD0BCiD37HgmCZAACZAACZAACZAACZCAxwhQAHmswJhcEiABEiABEiABEiABEiAB9wQogNyz45kkQAIkQAIkQAIkQAIkQAIeI0AB5LECY3JJgARIgARIgARIgARIgATcE6AAcs+OZ5IACZAACZAACZAACZAACXiMAAWQxwqMySUBEiABEiABEiABEiABEnBPgALIPTueSQIkQAIkQAIkQAIkQAIk4DECFEAeKzAmlwRIgARIgARIgARIgARIwD0BCiD37HgmCZAACZAACZAACZAACZCAxwhQAHmswJhcEiABEiABEiABEiABEiAB9wQogNyz45kkQAIkQAIkQAIkQAIkQAIeI0AB5LECY3JJgARIgARIgARIgARIgATcE6AAcs+OZ5IACZAACZAACZAACZAACXiMAAWQxwqMySUBEiABEiABEiABEiABEnBPgALIPTueSQIkQAIkQAIkQAIkQAIk4DEC/x+sEhkPVzlJdQAAAABJRU5ErkJggg=="}}},{"metadata":{"trusted":true},"cell_type":"code","source":"import random\nfrom random import sample\nimport numpy as np\nimport pandas as pd\nfrom sklearn.model_selection import StratifiedKFold\nimport gc\nimport torch.nn as nn\nfrom torch.utils.data import Dataset\nimport torch\nimport matplotlib.pyplot as plt\nfrom scipy.ndimage.filters import uniform_filter1d\nfrom scipy.interpolate import interp1d\nimport numpy as np \nimport pandas as pd  \nimport torch.nn as nn\nimport os\nimport time\nimport torch.optim as optim\nfrom torch.utils.data import DataLoader\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom collections import defaultdict\n\nEPOCH = 50\nBATCH_SIZE = 256\nDEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nMEMORY_LENGTH = 5 # max sequence length\nWIFI_NUM = 50\nFOLDS = 10\ninpath = '../input/indoor-location-navigation/'\nmetapath = inpath + 'metadata'\ntrainpath = inpath + 'train'\ntestpath = inpath + 'test'\n\ndef set_seed(seed=42):\n    random.seed(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed_all(seed)\n    torch.backends.cudnn.deterministic = True\n\nset_seed()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission = pd.read_csv('../input/indoor-location-navigation/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Procedure to produce floor results through KMeans+GBM @ https://www.kaggle.com/oxzplvifi/indoor-kmeans-gbm-floor-prediction"},{"metadata":{"trusted":true},"cell_type":"code","source":"floors = pd.read_csv('/kaggle/input/indoor-xy-floor/result_floor_feb22.csv',index_col=0)\nfloors.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The main difference between my dataset and [@kouki's dataset](https://www.kaggle.com/kokitanisaka/indoorunifiedwifids) is t1_wifi column.\nWe can treat as time series data by this column.\nI haven't deleted bssids by the frequency of bssids. "},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = pd.read_csv('../input/time-series-unified-wifi/train_all.csv',index_col=0)\ntest_data = pd.read_csv('../input/time-series-unified-wifi/test_all.csv',index_col=0)\ntrain_data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Long sequence is not learnable. make sure each paths length is shorter than memory_length. "},{"metadata":{"trusted":true},"cell_type":"code","source":"del_path = []\nnew_trains = []\nfor i, (p, g) in enumerate(train_data.groupby('path_id')):\n    if len(g) > MEMORY_LENGTH:\n        del_path.append(p)\n        for j in range((len(g) // MEMORY_LENGTH) + 1):\n            if j == (len(g) // MEMORY_LENGTH):\n                tmp = g.iloc[j*MEMORY_LENGTH:]\n            else:\n                tmp = g.iloc[j*MEMORY_LENGTH:(j+1)*MEMORY_LENGTH]\n            tmp.loc[:, 'path_id'] = p + '_' + str(j)\n            new_trains.append(tmp)\n        \ntrain_data.drop(train_data[train_data['path_id'].isin(del_path)].index, inplace=True)\ntrain_data = pd.concat([train_data, pd.concat(new_trains)]).reset_index(drop=True)\nprint(f\"previous path len:{i+1}, current path len:{len(train_data.groupby('path_id'))}\")\ndel new_trains\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"preparation"},{"metadata":{"trusted":true},"cell_type":"code","source":"BSSID_FEATS = [f'wb_{i}' for i in range(WIFI_NUM)]\nRSSI_FEATS = [f'wr_{i}' for i in range(WIFI_NUM)]\nX_train = train_data.loc[:, ['t1_wifi', 'building',\n                             'path_id'] + BSSID_FEATS + RSSI_FEATS]\ny_train = train_data.loc[:, ['t1_wifi', 'path_id', 'x', 'y', 'building']]\nX_test = test_data.loc[:, ['t1_wifi', 'building',\n                           'path_id'] + BSSID_FEATS + RSSI_FEATS]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### building weight\ntest_data frequency / train_data frequency "},{"metadata":{"trusted":true},"cell_type":"code","source":"test_building_weight = defaultdict(int)\ntrain_building_weight = defaultdict(int)\nbuilding_weight = dict()\n\nfor building in [x.split('_')[0] for x in sample_submission['site_path_timestamp'].values]:\n    test_building_weight[building] += 1 * 24/len(sample_submission)\n    \nfor building in train_data['building'].values:\n    train_building_weight[building] += 1\ntrain_building_weight = dict((k, len(train_data)/24/v) for k, v in train_building_weight.items())\n\nfor building in list(train_building_weight.keys()):\n    building_weight[building] = train_building_weight[building] * test_building_weight[building]\n    \ny_train['weight'] = [building_weight[x] for x in y_train['building'].values]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### normalization"},{"metadata":{"trusted":true},"cell_type":"code","source":"le = LabelEncoder()\nunique_bssids = np.unique(X_train.loc[:, BSSID_FEATS].values.tolist(\n) + X_test.loc[:, BSSID_FEATS].values.tolist())\nwifi_bssids_size = len(unique_bssids)\nprint(f\"bssids_size: {wifi_bssids_size}\")\nle.fit(unique_bssids)\nle_site = LabelEncoder()\nle_site.fit(list(set(X_train.loc[:, 'building'].values.tolist())))\n\nx_min = np.min(y_train.loc[:, 'x'].values)\ny_min = np.min(y_train.loc[:, 'y'].values)\nnorm_x = np.max(y_train.loc[:, 'x'].values) - x_min\nnorm_y = np.max(y_train.loc[:, 'y'].values) - y_min\n\nX_train.loc[:, RSSI_FEATS] = (X_train.loc[:, RSSI_FEATS] + 99) / (np.max(X_train.loc[:, RSSI_FEATS].values) + 99)\nX_test.loc[:, RSSI_FEATS] = (X_test.loc[:, RSSI_FEATS] +99) / (np.max(X_train.loc[:, RSSI_FEATS].values) + 99)\ny_train.loc[:, 'x'] = (y_train.loc[:, 'x'] - x_min) / norm_x\ny_train.loc[:, 'y'] = (y_train.loc[:, 'y'] - y_min) / norm_y\n\n\nfor i in BSSID_FEATS:\n    X_train.loc[:, i] = le.transform(X_train.loc[:, i])\n    X_test.loc[:, i] = le.transform(X_test.loc[:, i])\nX_train.loc[:, 'building'] = le_site.transform(X_train.loc[:, 'building'])\nX_test.loc[:, 'building'] = le_site.transform(X_test.loc[:, 'building'])\n\ntest_data.set_index(['path_id', 't1_wifi'], inplace=True)\nresult = pd.DataFrame(np.zeros([len(test_data), 2]),\n                     index=test_data.index, columns=['x', 'y'])\ndel train_data, test_data\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"group by path_id and sort by t1_wifi"},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_data(data, col):\n    train = []\n    for _, group in data.groupby('path_id'):\n        group = group.sort_values('t1_wifi')\n        train.append(group[col])\n    return train\n\ntrain_cols = ['building'] + BSSID_FEATS + RSSI_FEATS \ny_train = make_data(y_train, col=['x', 'y', 'weight'])\nX_train = make_data(X_train, col=train_cols)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Dataset\nI used inverse trajectory and random trajectory as augmentation according to https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7374298/"},{"metadata":{"trusted":true},"cell_type":"code","source":"class CustomDataset(Dataset):\n    def __init__(self, x_train, y_train, transform, inverse_ratio=0.5, combine_ratio=0.2):\n        self.transform = transform\n        self.x_train = x_train\n        self.y_train = y_train\n        self.inverse_ratio = inverse_ratio\n        self.combine_ratio = combine_ratio\n\n    def __getitem__(self, index):\n        x = self.x_train[index]\n        y = self.y_train[index]\n        mask = torch.tensor(np.full(MEMORY_LENGTH, True, dtype=bool))\n        if len(x) < MEMORY_LENGTH:\n            x_out = torch.tensor(\n                np.pad(x.to_numpy(), ([(0, MEMORY_LENGTH - len(x)), (0, 0)]), 'edge'))\n            y_out = torch.tensor(np.pad(\n                y.to_numpy(), ([(0, MEMORY_LENGTH - len(x)), (0, 0)]), 'edge'))\n            mask[len(x):] = False\n        else:\n            x_out = torch.tensor(x.to_numpy())\n            y_out = torch.tensor(y.to_numpy())\n\n        if self.transform:\n            # inverse trajectory\n            if np.random.rand() < self.inverse_ratio:\n                tmp = torch.arange(MEMORY_LENGTH-1, -1, -1)\n                x_out = x_out[tmp, :]\n                y_out = y_out[tmp, :]\n                mask = mask[tmp]\n\n            # # combine trajectory\n            if np.random.rand() < self.combine_ratio:\n                p = torch.randperm(MEMORY_LENGTH)\n                x_out = x_out[p]\n                y_out = y_out[p]\n                mask = mask[p]\n\n        return x_out, y_out, mask\n\n    def __len__(self):\n        return len(self.x_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### model\nmany to many"},{"metadata":{"trusted":true},"cell_type":"code","source":"class ManyToMany(nn.Module):\n    def __init__(self, wifi_bssids_size, input_dim, hidden_dim):\n        super(ManyToMany, self).__init__()\n        self.emb_dim = 256\n        self.entitity_dim = 128\n        self.bssi = nn.Embedding(wifi_bssids_size, self.emb_dim)\n        self.building = nn.Embedding(24, 2)\n\n        self.entity1 = nn.Sequential(\n            nn.Linear(self.emb_dim + 3, self.entitity_dim),\n            nn.Tanh()\n        )\n        self.entity2 = nn.Sequential(\n            nn.Dropout(0.2),\n            nn.Linear(self.entitity_dim*WIFI_NUM, input_dim),\n        )\n        self.gru = nn.GRU(input_dim, hidden_dim, num_layers=1,batch_first=True)\n\n        self.main = nn.Sequential(\n            nn.Linear(hidden_dim, 32),\n            nn.ReLU(True),\n            nn.Linear(32, 2),\n        )\n\n    def forward(self, x):\n        bssids = self.bssi(x[:, :, 1:int(WIFI_NUM+1)].long()).float()\n        buildings = self.building(x[:, :, 0].long()).unsqueeze(2).expand(-1, -1, WIFI_NUM, -1).float()\n        rssis = x[:, :, int(WIFI_NUM+1):int(WIFI_NUM*2+1)].unsqueeze(3).float()\n        \n        x = torch.cat((bssids,  rssis, buildings), axis=3)\n        #(batch, memory_length, wifi_num, self.emb_dim+2+1)\n        x = self.entity1(x)\n        #(batch, memory_length, wifi_num, self.entity1_dim)\n        x = x.flatten(start_dim=2)\n        #(batch, memory_length, wifi_num * self.entity1_dim)\n        x = self.entity2(x)\n        #(batch, memory_length, input_dim)\n        output, _ = self.gru(x)\n        output = self.main(output)\n        #(batch, memory_length, 2)\n        return output\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"utils"},{"metadata":{"trusted":true},"cell_type":"code","source":"def train_manytomany(t_loader, v_loader, optimizer, criterion, model,norms):\n    scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.6)\n    model.train()\n    early_stopping = 0\n    loss_pred = 100000\n    for e in range(EPOCH):\n        for (x, y, m) in t_loader:\n            optimizer.zero_grad()\n            x, y = x.to(DEVICE), y.to(DEVICE).float()\n            output = model(x)\n            loss = criterion(output[m], y[m])\n            loss.backward()\n            optimizer.step()\n        scheduler.step()\n        if (e+1) % 1 == 0:\n            losses = []\n            losses2 = []\n            total_point = 0\n            with torch.no_grad():\n                for (x, y, m) in v_loader:\n                    x, y = x.to(DEVICE), y.to(DEVICE).float()\n                    output = model(x)\n                    loss = comp_metric_xy(output[m], y[m], norms)\n                    loss2 = comp_metric_weighted_xy(output[m], y[m], norms)\n                    losses.append(loss.item())\n                    losses2.append(loss2.item())\n                    total_point += torch.count_nonzero(m)\n                tmp = np.sum(losses)/total_point\n                tmp2 = np.sum(losses2)/total_point\n                print(\n                    f\"epoch: {e}, loss: {tmp:.5f}, weight_loss: {tmp2:.5f} lr:{scheduler.get_last_lr()[0]:.5f}\")\n                if loss_pred < tmp2:\n                    early_stopping += 1\n                else:\n                    early_stopping = 0\n                    loss_pred = tmp2\n                    model_pred = model\n                if early_stopping > 5:\n                    return model_pred, loss_pred    \n    return model_pred, loss_pred\n\ndef comp_metric_xy(output, y, norms):\n    return torch.sqrt(((norms*(output - y[:,:2]))**2).sum())\ndef weighted_mse_loss(output, y):\n    return (y[:, 2].reshape(-1,1)*(output - y[:,:2])**2).sum()\ndef comp_metric_weighted_xy(output, y, norms):\n    return torch.sqrt(((y[:, 2].reshape(-1,1)*(norms * (output - y[:,:2]))**2)).sum())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"train"},{"metadata":{"trusted":true},"cell_type":"code","source":"buildings = []\nfor x in X_train:\n    buildings.append(x.iloc[0, 0])\n\nloss_folds = []\nskf = StratifiedKFold(n_splits=FOLDS)\ncriterion = weighted_mse_loss\nnorms = torch.tensor([[norm_x, norm_y]]).to(DEVICE)\nstart = time.time()\nfor fold, (idt, idv) in enumerate(skf.split(X_train, buildings)):\n    print('\\r', f'{fold}', end='\\t')\n    mtrain, mvalid = [X_train[i] for i in idt], [y_train[i] for i in idt]\n    ltrain, lvalid = [X_train[i] for i in idv], [y_train[i] for i in idv]\n    t_loader = DataLoader(CustomDataset(\n        mtrain, mvalid, transform=True), batch_size=BATCH_SIZE, shuffle=True)\n    v_loader = DataLoader(CustomDataset(\n        ltrain, lvalid, transform=False), batch_size=BATCH_SIZE, shuffle=False)\n\n    model = ManyToMany(wifi_bssids_size, 64, 128).to(DEVICE)\n    optimizer = optim.Adam(model.parameters(), lr=0.005)\n    \n    model, loss = train_manytomany(\n        t_loader, v_loader, optimizer, criterion, model, norms)\n\n    # prediction by sliding window averaging. see https://arxiv.org/pdf/1903.11703.pdf\n    prediction = []\n    with torch.no_grad():\n        for p, x in X_test.groupby('path_id'):\n            window_score = defaultdict(list)\n            x = x.sort_values('t1_wifi')\n            ts = x['t1_wifi'].to_numpy()\n            x = x[train_cols].reset_index(drop=True)\n            \n            for window in range(len(ts)):\n                if MEMORY_LENGTH + window > len(ts):\n                    break\n                \n                x_out = torch.tensor(\n                    x.iloc[window:(window+MEMORY_LENGTH)].to_numpy())\n                x_out = x_out.unsqueeze(0)\n                x_out = model(x_out.to(DEVICE))\n                for i in range(MEMORY_LENGTH):\n                    window_score[ts[window + i]].append(x_out.squeeze().cpu().detach().numpy()[i,:])\n                \n            #  sort by time and get average \n            window_score = sorted(window_score.items(), key=lambda x:x[0])\n            prediction.extend(list(map(lambda x: np.mean(x[1],axis=0), window_score)))\n\n    prediction = np.array(prediction)\n    result['x'] += (prediction[:, 0] * norm_x + x_min) / FOLDS\n    result['y'] += (prediction[:, 1] * norm_y + y_min) / FOLDS\n    loss_folds.append(loss)\n    print(f'fold: {fold}, loss: {loss}')\n    \n    del mtrain, mvalid, ltrain, lvalid, t_loader, v_loader, model, prediction\n    gc.collect()\n    \n\nprint(f\"mean loss: {np.mean(loss_folds)}, time:{(start -time.time())/60:.1f} min.\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Post processing is same as [Indoor GBM+postprocessing XY prediction](https://www.kaggle.com/oxzplvifi/indoor-gbm-postprocessing-xy-prediction) by [@Oscar Villarreal Escamilla](https://www.kaggle.com/oxzplvifi)<br>"},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_submission['building'] = [x.split('_')[0] for x in sample_submission['site_path_timestamp']]\nsample_submission['path_id'] = [x.split('_')[1] for x in sample_submission['site_path_timestamp']]\nsample_submission['timestamp'] = [x.split('_')[2] for x in sample_submission['site_path_timestamp']]\nsamples = pd.DataFrame(sample_submission.groupby(['building','path_id'])['timestamp'].apply(lambda x: list(x)))\nbuildings = np.unique([x[0] for x in samples.index])\nsamples.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"result.reset_index(inplace=True)\nresult.set_index('path_id', inplace=True)\nresult","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.spatial.transform import Rotation as R\nfrom PIL import Image\nfrom mpl_toolkits.mplot3d import Axes3D\nimport plotly.graph_objs as go\nfrom pathlib import Path\nimport scipy.signal as signal\nimport json\nimport seaborn as sns # visualization\nfrom dataclasses import dataclass\n\nimport matplotlib.pyplot as plt  # visualization\nimport numpy as np  # linear algebra\nimport random\nimport pandas as pd\nfrom collections import Counter, defaultdict\n\nplt.rcParams.update({'font.size': 14})\n\ndef split_ts_seq(ts_seq, sep_ts):\n    \"\"\"\n\n    :param ts_seq:\n    :param sep_ts:\n    :return:\n    \"\"\"\n    tss = ts_seq[:, 0].astype(float)\n    unique_sep_ts = np.unique(sep_ts)\n    ts_seqs = []\n    start_index = 0\n    for i in range(0, unique_sep_ts.shape[0]):\n        end_index = np.searchsorted(tss, unique_sep_ts[i], side='right')\n        if start_index == end_index:\n            continue\n        ts_seqs.append(ts_seq[start_index:end_index, :].copy())\n        start_index = end_index\n\n    # tail data\n    if start_index < ts_seq.shape[0]:\n        ts_seqs.append(ts_seq[start_index:, :].copy())\n\n    return ts_seqs\n\n\ndef correct_trajectory(original_xys, end_xy):\n    \"\"\"\n\n    :param original_xys: numpy ndarray, shape(N, 2)\n    :param end_xy: numpy ndarray, shape(1, 2)\n    :return:\n    \"\"\"\n    corrected_xys = np.zeros((0, 2))\n\n    A = original_xys[0, :]\n    B = end_xy\n    Bp = original_xys[-1, :]\n\n    angle_BAX = np.arctan2(B[1] - A[1], B[0] - A[0])\n    angle_BpAX = np.arctan2(Bp[1] - A[1], Bp[0] - A[0])\n    angle_BpAB = angle_BpAX - angle_BAX\n    AB = np.sqrt(np.sum((B - A) ** 2))\n    ABp = np.sqrt(np.sum((Bp - A) ** 2))\n\n    corrected_xys = np.append(corrected_xys, [A], 0)\n    for i in np.arange(1, np.size(original_xys, 0)):\n        angle_CpAX = np.arctan2(original_xys[i, 1] - A[1], original_xys[i, 0] - A[0])\n\n        angle_CAX = angle_CpAX - angle_BpAB\n\n        ACp = np.sqrt(np.sum((original_xys[i, :] - A) ** 2))\n\n        AC = ACp * AB / ABp\n\n        delta_C = np.array([AC * np.cos(angle_CAX), AC * np.sin(angle_CAX)])\n\n        C = delta_C + A\n\n        corrected_xys = np.append(corrected_xys, [C], 0)\n\n    return corrected_xys\n\n\ndef correct_positions(rel_positions, reference_positions):\n    \"\"\"\n\n    :param rel_positions:\n    :param reference_positions:\n    :return:\n    \"\"\"\n    rel_positions_list = split_ts_seq(rel_positions, reference_positions[:, 0])\n    if len(rel_positions_list) != reference_positions.shape[0] - 1:\n        # print(f'Rel positions list size: {len(rel_positions_list)}, ref positions size: {reference_positions.shape[0]}')\n        del rel_positions_list[-1]\n    assert len(rel_positions_list) == reference_positions.shape[0] - 1\n\n    corrected_positions = np.zeros((0, 3))\n    for i, rel_ps in enumerate(rel_positions_list):\n        start_position = reference_positions[i]\n        end_position = reference_positions[i + 1]\n        abs_ps = np.zeros(rel_ps.shape)\n        abs_ps[:, 0] = rel_ps[:, 0]\n        # abs_ps[:, 1:3] = rel_ps[:, 1:3] + start_position[1:3]\n        abs_ps[0, 1:3] = rel_ps[0, 1:3] + start_position[1:3]\n        for j in range(1, rel_ps.shape[0]):\n            abs_ps[j, 1:3] = abs_ps[j-1, 1:3] + rel_ps[j, 1:3]\n        abs_ps = np.insert(abs_ps, 0, start_position, axis=0)\n        corrected_xys = correct_trajectory(abs_ps[:, 1:3], end_position[1:3])\n        corrected_ps = np.column_stack((abs_ps[:, 0], corrected_xys))\n        if i == 0:\n            corrected_positions = np.append(corrected_positions, corrected_ps, axis=0)\n        else:\n            corrected_positions = np.append(corrected_positions, corrected_ps[1:], axis=0)\n\n    corrected_positions = np.array(corrected_positions)\n\n    return corrected_positions\n\n\ndef init_parameters_filter(sample_freq, warmup_data, cut_off_freq=2):\n    order = 4\n    filter_b, filter_a = signal.butter(order, cut_off_freq / (sample_freq / 2), 'low', False)\n    zf = signal.lfilter_zi(filter_b, filter_a)\n    _, zf = signal.lfilter(filter_b, filter_a, warmup_data, zi=zf)\n    _, filter_zf = signal.lfilter(filter_b, filter_a, warmup_data, zi=zf)\n\n    return filter_b, filter_a, filter_zf\n\n\ndef get_rotation_matrix_from_vector(rotation_vector):\n    q1 = rotation_vector[0]\n    q2 = rotation_vector[1]\n    q3 = rotation_vector[2]\n\n    if rotation_vector.size >= 4:\n        q0 = rotation_vector[3]\n    else:\n        q0 = 1 - q1*q1 - q2*q2 - q3*q3\n        if q0 > 0:\n            q0 = np.sqrt(q0)\n        else:\n            q0 = 0\n\n    sq_q1 = 2 * q1 * q1\n    sq_q2 = 2 * q2 * q2\n    sq_q3 = 2 * q3 * q3\n    q1_q2 = 2 * q1 * q2\n    q3_q0 = 2 * q3 * q0\n    q1_q3 = 2 * q1 * q3\n    q2_q0 = 2 * q2 * q0\n    q2_q3 = 2 * q2 * q3\n    q1_q0 = 2 * q1 * q0\n\n    R = np.zeros((9,))\n    if R.size == 9:\n        R[0] = 1 - sq_q2 - sq_q3\n        R[1] = q1_q2 - q3_q0\n        R[2] = q1_q3 + q2_q0\n\n        R[3] = q1_q2 + q3_q0\n        R[4] = 1 - sq_q1 - sq_q3\n        R[5] = q2_q3 - q1_q0\n\n        R[6] = q1_q3 - q2_q0\n        R[7] = q2_q3 + q1_q0\n        R[8] = 1 - sq_q1 - sq_q2\n\n        R = np.reshape(R, (3, 3))\n    elif R.size == 16:\n        R[0] = 1 - sq_q2 - sq_q3\n        R[1] = q1_q2 - q3_q0\n        R[2] = q1_q3 + q2_q0\n        R[3] = 0.0\n\n        R[4] = q1_q2 + q3_q0\n        R[5] = 1 - sq_q1 - sq_q3\n        R[6] = q2_q3 - q1_q0\n        R[7] = 0.0\n\n        R[8] = q1_q3 - q2_q0\n        R[9] = q2_q3 + q1_q0\n        R[10] = 1 - sq_q1 - sq_q2\n        R[11] = 0.0\n\n        R[12] = R[13] = R[14] = 0.0\n        R[15] = 1.0\n\n        R = np.reshape(R, (4, 4))\n\n    return R\n\n\ndef get_orientation(R):\n    flat_R = R.flatten()\n    values = np.zeros((3,))\n    if np.size(flat_R) == 9:\n        values[0] = np.arctan2(flat_R[1], flat_R[4])\n        values[1] = np.arcsin(-flat_R[7])\n        values[2] = np.arctan2(-flat_R[6], flat_R[8])\n    else:\n        values[0] = np.arctan2(flat_R[1], flat_R[5])\n        values[1] = np.arcsin(-flat_R[9])\n        values[2] = np.arctan2(-flat_R[8], flat_R[10])\n\n    return values\n\n\ndef compute_steps(acce_datas):\n    step_timestamps = np.array([])\n    step_indexs = np.array([], dtype=int)\n    step_acce_max_mins = np.zeros((0, 4))\n    sample_freq = 50\n    window_size = 22\n    low_acce_mag = 0.6\n    step_criterion = 1\n    interval_threshold = 250\n\n    acce_max = np.zeros((2,))\n    acce_min = np.zeros((2,))\n    acce_binarys = np.zeros((window_size,), dtype=int)\n    acce_mag_pre = 0\n    state_flag = 0\n\n    warmup_data = np.ones((window_size,)) * 9.81\n    filter_b, filter_a, filter_zf = init_parameters_filter(sample_freq, warmup_data)\n    acce_mag_window = np.zeros((window_size, 1))\n\n    # detect steps according to acceleration magnitudes\n    for i in np.arange(0, np.size(acce_datas, 0)):\n        acce_data = acce_datas[i, :]\n        acce_mag = np.sqrt(np.sum(acce_data[1:] ** 2))\n\n        acce_mag_filt, filter_zf = signal.lfilter(filter_b, filter_a, [acce_mag], zi=filter_zf)\n        acce_mag_filt = acce_mag_filt[0]\n\n        acce_mag_window = np.append(acce_mag_window, [acce_mag_filt])\n        acce_mag_window = np.delete(acce_mag_window, 0)\n        mean_gravity = np.mean(acce_mag_window)\n        acce_std = np.std(acce_mag_window)\n        mag_threshold = np.max([low_acce_mag, 0.4 * acce_std])\n\n        # detect valid peak or valley of acceleration magnitudes\n        acce_mag_filt_detrend = acce_mag_filt - mean_gravity\n        if acce_mag_filt_detrend > np.max([acce_mag_pre, mag_threshold]):\n            # peak\n            acce_binarys = np.append(acce_binarys, [1])\n            acce_binarys = np.delete(acce_binarys, 0)\n        elif acce_mag_filt_detrend < np.min([acce_mag_pre, -mag_threshold]):\n            # valley\n            acce_binarys = np.append(acce_binarys, [-1])\n            acce_binarys = np.delete(acce_binarys, 0)\n        else:\n            # between peak and valley\n            acce_binarys = np.append(acce_binarys, [0])\n            acce_binarys = np.delete(acce_binarys, 0)\n\n        if (acce_binarys[-1] == 0) and (acce_binarys[-2] == 1):\n            if state_flag == 0:\n                acce_max[:] = acce_data[0], acce_mag_filt\n                state_flag = 1\n            elif (state_flag == 1) and ((acce_data[0] - acce_max[0]) <= interval_threshold) and (\n                    acce_mag_filt > acce_max[1]):\n                acce_max[:] = acce_data[0], acce_mag_filt\n            elif (state_flag == 2) and ((acce_data[0] - acce_max[0]) > interval_threshold):\n                acce_max[:] = acce_data[0], acce_mag_filt\n                state_flag = 1\n\n        # choose reasonable step criterion and check if there is a valid step\n        # save step acceleration data: step_acce_max_mins = [timestamp, max, min, variance]\n        step_flag = False\n        if step_criterion == 2:\n            if (acce_binarys[-1] == -1) and ((acce_binarys[-2] == 1) or (acce_binarys[-2] == 0)):\n                step_flag = True\n        elif step_criterion == 3:\n            if (acce_binarys[-1] == -1) and (acce_binarys[-2] == 0) and (np.sum(acce_binarys[:-2]) > 1):\n                step_flag = True\n        else:\n            if (acce_binarys[-1] == 0) and acce_binarys[-2] == -1:\n                if (state_flag == 1) and ((acce_data[0] - acce_min[0]) > interval_threshold):\n                    acce_min[:] = acce_data[0], acce_mag_filt\n                    state_flag = 2\n                    step_flag = True\n                elif (state_flag == 2) and ((acce_data[0] - acce_min[0]) <= interval_threshold) and (\n                        acce_mag_filt < acce_min[1]):\n                    acce_min[:] = acce_data[0], acce_mag_filt\n        if step_flag:\n            step_timestamps = np.append(step_timestamps, acce_data[0])\n            step_indexs = np.append(step_indexs, [i])\n            step_acce_max_mins = np.append(step_acce_max_mins,\n                                           [[acce_data[0], acce_max[1], acce_min[1], acce_std ** 2]], axis=0)\n        acce_mag_pre = acce_mag_filt_detrend\n\n    return step_timestamps, step_indexs, step_acce_max_mins\n\n\ndef compute_stride_length(step_acce_max_mins):\n    K = 0.4\n    K_max = 0.8\n    K_min = 0.4\n    para_a0 = 0.21468084\n    para_a1 = 0.09154517\n    para_a2 = 0.02301998\n\n    stride_lengths = np.zeros((step_acce_max_mins.shape[0], 2))\n    k_real = np.zeros((step_acce_max_mins.shape[0], 2))\n    step_timeperiod = np.zeros((step_acce_max_mins.shape[0] - 1, ))\n    stride_lengths[:, 0] = step_acce_max_mins[:, 0]\n    window_size = 2\n    step_timeperiod_temp = np.zeros((0, ))\n\n    # calculate every step period - step_timeperiod unit: second\n    for i in range(0, step_timeperiod.shape[0]):\n        step_timeperiod_data = (step_acce_max_mins[i + 1, 0] - step_acce_max_mins[i, 0]) / 1000\n        step_timeperiod_temp = np.append(step_timeperiod_temp, [step_timeperiod_data])\n        if step_timeperiod_temp.shape[0] > window_size:\n            step_timeperiod_temp = np.delete(step_timeperiod_temp, [0])\n        step_timeperiod[i] = np.sum(step_timeperiod_temp) / step_timeperiod_temp.shape[0]\n\n    # calculate parameters by step period and acceleration magnitude variance\n    k_real[:, 0] = step_acce_max_mins[:, 0]\n    k_real[0, 1] = K\n    for i in range(0, step_timeperiod.shape[0]):\n        k_real[i + 1, 1] = np.max([(para_a0 + para_a1 / step_timeperiod[i] + para_a2 * step_acce_max_mins[i, 3]), K_min])\n        k_real[i + 1, 1] = np.min([k_real[i + 1, 1], K_max]) * (K / K_min)\n\n    # calculate every stride length by parameters and max and min data of acceleration magnitude\n    stride_lengths[:, 1] = np.max([(step_acce_max_mins[:, 1] - step_acce_max_mins[:, 2]),\n                                   np.ones((step_acce_max_mins.shape[0], ))], axis=0)**(1 / 4) * k_real[:, 1]\n\n    return stride_lengths\n\n\ndef compute_headings(ahrs_datas):\n    headings = np.zeros((np.size(ahrs_datas, 0), 2))\n    for i in np.arange(0, np.size(ahrs_datas, 0)):\n        ahrs_data = ahrs_datas[i, :]\n        rot_mat = get_rotation_matrix_from_vector(ahrs_data[1:])\n        azimuth, pitch, roll = get_orientation(rot_mat)\n        around_z = (-azimuth) % (2 * np.pi)\n        headings[i, :] = ahrs_data[0], around_z\n    return headings\n\n\ndef compute_step_heading(step_timestamps, headings):\n    step_headings = np.zeros((len(step_timestamps), 2))\n    step_timestamps_index = 0\n    for i in range(0, len(headings)):\n        if step_timestamps_index < len(step_timestamps):\n            if headings[i, 0] == step_timestamps[step_timestamps_index]:\n                step_headings[step_timestamps_index, :] = headings[i, :]\n                step_timestamps_index += 1\n        else:\n            break\n    assert step_timestamps_index == len(step_timestamps)\n\n    return step_headings\n\n\ndef compute_rel_positions(stride_lengths, step_headings):\n    rel_positions = np.zeros((stride_lengths.shape[0], 3))\n    for i in range(0, stride_lengths.shape[0]):\n        rel_positions[i, 0] = stride_lengths[i, 0]\n        rel_positions[i, 1] = -stride_lengths[i, 1] * np.sin(step_headings[i, 1])\n        rel_positions[i, 2] = stride_lengths[i, 1] * np.cos(step_headings[i, 1])\n\n    return rel_positions\n\n\ndef compute_step_positions(acce_datas, ahrs_datas, posi_datas):\n    step_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\n    headings = compute_headings(ahrs_datas)\n    stride_lengths = compute_stride_length(step_acce_max_mins)\n    step_headings = compute_step_heading(step_timestamps, headings)\n    rel_positions = compute_rel_positions(stride_lengths, step_headings)\n    step_positions = correct_positions(rel_positions, posi_datas)\n\n    return step_positions\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Extract testing files, buildings and sites:\nos.system(f'grep SiteID {testpath}/* > test_buildings.txt' )\ntest_buildings = pd.read_csv('test_buildings.txt',sep='\\t',header=None,names=['file','building','site'])\ntest_buildings['file'] = test_buildings['file'].apply(lambda x: x[:-2])\ntest_buildings['building'] = test_buildings['building'].apply(lambda x: x[7:])\n\n# How many buildings in the testing set?\nbuildings = np.unique(test_buildings['building'])\nprint('There are',len(buildings),'buildings in the testing set.')\n\ntest_buildings.head()\n# Compile C++ pre-processing code:\ner=os.system(\"g++ /kaggle/input/indoor-cpp/1_preprocess.cpp -std=c++11 -o preprocess\")\nif(er): print(\"Error\")\n\n# Reformat the testing set:\nos.system('mkdir test')\nfor i,(path_filename,building) in enumerate(zip(test_buildings['file'],test_buildings['building'])):\n    er=os.system(f'./preprocess {path_filename} test {building} {0}') #since we do not know the floor, I put 0.\n    if(er): print(\"Error:\",path_filename)\n# Acceleration, magnetic and orientation testing data:\nos.system('mkdir indoor_testing_accel')\nos.system(\"g++ /kaggle/input/indoor-cpp/2_preprocess_accel.cpp -std=c++11 -o preprocess_accel\")\nfor building in buildings:\n    os.system(f'./preprocess_accel {building}')\n# Wifi testing data:\nos.system('mkdir test_wifi')\nos.system(\"g++ /kaggle/input/indoor-cpp/2_preprocess_wifi.cpp -std=c++11 -o preprocess_wifi\")\nfor building in buildings:\n    os.system(f'./preprocess_wifi {building}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.interpolate import interp1d\nfrom scipy.ndimage.filters import uniform_filter1d\n\ncolacce = ['xyz_time','x_acce','y_acce','z_acce']\ncolahrs = ['xyz_time','x_ahrs','y_ahrs','z_ahrs']\n\nfor building in buildings:\n    print(building)\n    paths = samples.loc[building].index\n    # Acceleration info:\n    tfm = pd.read_csv(f'indoor_testing_accel/{building}.txt',index_col=0)\n    for path_id in paths:\n        # Original predicted values:\n        xy = result.loc[building+'_'+path_id]\n        tfmi = tfm.loc[path_id]\n        acce_datas = np.array(tfmi[colacce],dtype=np.float)\n        ahrs_datas = np.array(tfmi[colahrs],dtype=np.float)\n        posi_datas = np.array(xy[['t1_wifi','x','y']],dtype=np.float)\n        # Outlier removal:\n        xyout = uniform_filter1d(posi_datas,size=3,axis=0,mode='reflect')\n        xydiff = np.abs(posi_datas-xyout)\n        xystd = np.std(xydiff,axis=0)*3\n        posi_datas = posi_datas[(xydiff[:,1]<xystd[1])&(xydiff[:,2]<xystd[2])]\n        # Step detection:\n        step_timestamps, step_indexs, step_acce_max_mins = compute_steps(acce_datas)\n        stride_lengths = compute_stride_length(step_acce_max_mins)\n        # Orientation detection:\n        headings = compute_headings(ahrs_datas)\n        step_headings = compute_step_heading(step_timestamps, headings)\n        rel_positions = compute_rel_positions(stride_lengths, step_headings)\n        # Running average:\n        posi_datas = uniform_filter1d(posi_datas,size=3,axis=0,mode='reflect')[0::3,:]\n        # The 1st prediction timepoint should be earlier than the 1st step timepoint.\n        rel_positions = rel_positions[rel_positions[:,0]>posi_datas[0,0],:]\n        # If two consecutive predictions are in-between two step datapoints,\n        # the last one is removed, causing error (in the \"split_ts_seq\" function).\n        posi_index = [np.searchsorted(rel_positions[:,0], x, side='right') for x in posi_datas[:,0]]\n        u, i1, i2 = np.unique(posi_index, return_index=True, return_inverse=True)\n        posi_datas = np.vstack([np.mean(posi_datas[i2==i],axis=0) for i in np.unique(i2)])\n        # Position correction:\n        step_positions = correct_positions(rel_positions, posi_datas)\n        # Interpolate for timestamps in the testing set:\n        t = step_positions[:,0]\n        x = step_positions[:,1]\n        y = step_positions[:,2]\n        fx = interp1d(t, x, kind='linear', fill_value=(x[0],x[-1]), bounds_error=False) #fill_value=\"extrapolate\"\n        fy = interp1d(t, y, kind='linear', fill_value=(y[0],y[-1]), bounds_error=False)\n        # Output result:\n        t0 = np.array(samples.loc[(building,path_id),'timestamp'],dtype=np.float64)\n        sample_submission.loc[(sample_submission.building==building)&(sample_submission.path_id==path_id),'x'] = fx(t0)\n        sample_submission.loc[(sample_submission.building==building)&(sample_submission.path_id==path_id),'y'] = fy(t0)\n        sample_submission.loc[(sample_submission.building==building)&(sample_submission.path_id==path_id),'floor'] = floors.loc[building+'_'+path_id,'floor']\n#         break\n#     break\n\nsample_submission[['site_path_timestamp','floor','x','y']].to_csv('submission.csv',index=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Thank you for reading! \nI'm new to kaggle and RNN, and this is my first published notebook.\nSo please let me know if you notice any mistakes or have suggestions."}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}