Skip to content

Instantly share code, notes, and snippets.

@ruoyu0088

ruoyu0088/numba-skeletonize.ipynb

Created February 9, 2021 01:39
Show Gist options
  • Save ruoyu0088/5678aaacdc198e901cd74d572bdfc998 to your computer and use it in GitHub Desktop.
Save ruoyu0088/5678aaacdc198e901cd74d572bdfc998 to your computer and use it in GitHub Desktop.
Download ZIP
numba-skeletonize.ipynb
Display the source blob
Display the rendered blob
Raw
numba-skeletonize.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"from numba.typed import List, Dict\n",
"from numba import njit, int64, int32\n",
"%matplotlib inline\n",
"from matplotlib import pyplot as plt\n",
"from PIL import Image\n",
"import numba as nb"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Wall time: 396 ms\n"
]
},
{
"data": {
"text/plain": [
"array([[ 0. , 0. , 0. , ..., 0. ,\n",
" 0. , 0. ],\n",
" [ 0. , -0.61544942, -0.36143147, ..., -0.19687525,\n",
" -0.19397574, 0. ],\n",
" [ 0. , -0.06983371, -0.17894958, ..., 0.41112256,\n",
" -0.69122705, 0. ],\n",
" ...,\n",
" [ 0. , -0.15352685, -0.887787 , ..., -0.08003802,\n",
" -0.59011623, 0. ],\n",
" [ 0. , -0.40265092, -0.14167998, ..., -0.39880908,\n",
" -0.4883661 , 0. ],\n",
" [ 0. , 0. , 0. , ..., 0. ,\n",
" 0. , 0. ]])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@nb.stencil\n",
"def kernel(a):\n",
" return 0.2 * (a[0, 0] + a[0, 1] + a[1, 0] + a[0, -1] + a[-1, 0])\n",
"\n",
"arr = np.random.randn(1000, 1000)\n",
"%time kernel(arr)"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [],
"source": [
"@nb.stencil\n",
"def kernel(a):\n",
" a00 = a[-1, -1]\n",
" a01 = a[-1, 0]\n",
" a02 = a[-1, 1]\n",
" a10 = a[0, -1]\n",
" a11 = a[0, 0]\n",
" a12 = a[0, 1]\n",
" a20 = a[1, -1]\n",
" a21 = a[1, 0]\n",
" a22 = a[1, 1]\n",
" if a00 == 0 and a01 == 0 and a02 == 0 and a11 == 1 and a20 == 1 and a21 == 1 and a22 == 1: return 1\n",
" if a01 == 0 and a02 == 0 and a10 == 1 and a11 == 1 and a12 == 0 and a21 == 1: return 1\n",
" if a00 == 0 and a02 == 1 and a10 == 0 and a11 == 1 and a12 == 1 and a20 == 0 and a22 == 1: return 1\n",
" if a00 == 0 and a01 == 0 and a10 == 0 and a11 == 1 and a12 == 1 and a21 == 1: return 1\n",
" if a00 == 1 and a01 == 1 and a02 == 1 and a11 == 1 and a20 == 0 and a21 == 0 and a22 == 0: return 1\n",
" if a01 == 1 and a10 == 0 and a11 == 1 and a12 == 1 and a20 == 0 and a21 == 0: return 1\n",
" if a00 == 1 and a02 == 0 and a10 == 1 and a11 == 1 and a12 == 0 and a20 == 1 and a22 == 0: return 1\n",
" if a01 == 1 and a10 == 1 and a11 == 1 and a12 == 0 and a21 == 0 and a22 == 0: return 1\n",
" return 0"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x2a2379fc250>"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeOElEQVR4nO3deXxU9b3/8ddndrJBQhIIEELYwr6ZClhFLWrdbrG3pVfrXlra29al1gXae3/drNf+2lp7+/NXy9VaXNrrRtW6wqVqsSKLiAgEDEuAsCUkBEgCycw53/tHpgFNMDPJrDmf5+ORRzhnzpnzyTDvOed853u+R4wxKKV6P1eyC1BKJYaGXSmH0LAr5RAadqUcQsOulENo2JVyiB6FXUQuFpGtIrJNRBbEqiilVOxJd79nFxE38CFwIVANrAGuMsZsjl15SqlY8fRg3TOBbcaYHQAi8t/AHOC0YfeJ3wTI7MEmlVKf5ARNtJoW6eyxnoR9MLDnlOlqYPrHFxKR+cB8gAAZTJfZPdikUuqTrDLLT/tYT87ZO/v06HBOYIxZZIwpN8aUe/H3YHNKqZ7oSdirgeJTpocA+3pWjlIqXnoS9jXAKBEpFREfcCXwQmzKUkrFWrfP2Y0xIRH5NvAa4AZ+b4zZFLPKlFIx1ZMGOowxLwMvx6gWpVQcaQ86pRxCw66UQ2jYlXIIDbtSDqFhV8ohNOxKOYSGXSmH0LAr5RAadqUcQsOulENo2JVyCA27Ug6hYVfKITTsSjmEhl0ph9CwK+UQGnalHELDrpRDaNiVcggNu1IOoWFXyiE07Eo5hIZdKYfQsCvlEBp2pRxCw66UQ2jYlXIIDbtSDqFhV8ohNOxKOYSGXSmH0LAr5RAadqUcQsOulENo2JVyiC7DLiLFIvK6iFSIyCYRuSU8P09ElolIZfh3bvzLVUp1VyR79hDwXWPMWGAG8C0RGQcsAJYbY0YBy8PTSqkU1WXYjTH7jTHrwv8+BlQAg4E5wOLwYouBK+JUo1IqBqI6ZxeRYcBUYBUwwBizH9o+EIDC06wzX0TWisjaIC09LFcp1V0Rh11EsoBngVuNMUcjXc8Ys8gYU26MKffi706NSqkYiCjsIuKlLehPGGOWhGcfFJGi8ONFQE18SlRKxUIkrfECPAxUGGPuO+WhF4Drw/++Hng+9uUppWLFE8EynwauBT4QkfXhed8D7gWeEpF5wG5gblwqVErFRJdhN8a8BchpHp4d23KUUvGiPeiUcggNu1IOEck5u1LOJac7g/0YY+JbRwxo2JXqhCs7mwPXT+ToSLvLZb1HheF/OoRVUZmAyrpPw65UJ1xZmRR+fjfvjX2xy2Vfag7w81XX4q9IQGE9oGFX6hTuggIOXTqSxmLhmoLlEa1T7Glg9yUucoafRdHr9dgbt8S5yu7RsCt1CntoIefcvIqFhSvo6woA7i7XGe/18e6cX7Ej5GF+6y3kb4x/nd2hYVfqVCLkepvJd2dGvIpbXOS6MxhgGrHdETboJYF+9aaUQ+iePQ14hgzGGhi/gYDcew8R2n8gbs+vUoOGPdWJsPPGEr48969xefqgcfPcH85l4K807L2dhj3VuNy4c/siPl942sXxklb+LT8+LbxBY/H40HPwFA3s8JgJWVh19WBbcdm2SiwNe4rxDBrI5n8bzOhR+wAQMfxo0HNx255X3Nzx2b/w50lTOzxWuXcQY37aF2vrtrhtXyWOhj3ZRBCPF1xtrbh2v2wuO+N9/t/gVQkr4Rv99vKNfns7zL9/4DBeKjgPd1V4hCHbYELBtOgaqjrSsCeZPWsK277sQfq0HSr3yWzhrrw1Sa6qzWcyt/DondNpaBgPgAm6KHlG8L+SGvWp6GjYk+zwqACvXfwLRnsj/143USb5Aqwrf7J9+rDVzKyK2xn0ShKLUt2mYU+S0Owz2DvLj2fiEfLSpLeDXzxkfqaGqr4zKXonhP8l3cOnEw17kuyd5WfVvF+SIT68knp79c5kuHysmPwk9mSbsTk3Mepl0fP3NKJhTwBXIMDx8yfSVHTy5XZPOBIOetd9r1NJW71uispqqJs3g5yqVnxvvI8JhZJdmuqChj0BXP360vztBp6Y8If2eXku0maP3pkXxz9O/VibS9/5JiNWZ2AdjfhWAipJNOxx5MrMxJo0koYhAaYWbIh5I9whq4nHj45nf2u/Lpcd6q/jmpwP6evqE5Nt57ozyHVDZp+W9q8NVWrTsMfTiGJcPz3EHcVvMt1/AMiK6dO/3FTCk3dfTL9NDV0uu/zMXEYteIiLMoIxrUGlDw17HBmfh7Pzt3NFZiOxCHp1qJGNrf2xwiN7v1o3kb6Vjdgbuu5K26//NF5smEKQDz4y342hzFtHqTe2H0Qq9WjY08i8yitpfmAw7pa2cdE8TRa+D7cRSc9134YqNt01iQ2BKR+Zb1xCw7xjvH/mn2JfsEopGvY4EI8HV0YGrVk+/K7uHzYftpppMicHPKzcW0jZqx9gNze3z4v0EhWrrh7P8vqO/+EuN8cuLKd6WmP7rIBIVIM3qPSgYY8DUz6OLd/0MGpwDZdmbQSibxQ7bDVzxuvfpu+qQPu8km2t2CdifNtr22LEs63809Y722c1TAzx10vu00P7XkbDHmsiNA3pwx/O/h2zAtCdoFvGpsG26bfST8Fv3455iR/nevM9Ct88OR24agYHL+rDUM/Jowq3dOzmZxkbY7QlPl1o2GPIOn8aOz/nI6v0CMM83WuU+zDYxOVvfxPZkUHJ+uauV4iDfpuPcMNjN2H1aesdZ2Xa/OSCZ7g6u659mRt3n8OKFRPI3Qz28Z1JqVNFR8MeQ7VTAqz84s8pdGfS3db3ymB/Bj3uw//SytgWFwX7/QpK3j857Skp5vnJU7g6++TQym+uGceou1aDbaEdZtODhj3GuntNy+qWIN/e/GVq9+QyZm8jXd+HJHFMYxMblo+hbO/gtmkD/de5wKRSlaorGvYU8VzDGeT9MED/DzZgt8S4Ea6HrLp6ht3zLuI+2Y/ftLZi9CKYtKJhTxFB40ZagtgnTiS7lE6ZlpaUOFx3ZWdz4qwyWvtG/tb1HLfJWrOL0IGDcaws9WnYVXopGUzBv+/ktkFLI17lr43jWHrXLPyvaNhVD3mKhxAqyqW5yOAmuq+idgYbWdY8mjf2jaLghPZb/zh3QQH2sIGY8J1WGkozuDR3JTMC0VwavJkny2Yz8PCktknb4NlXT6i647h7vVnEYRcRN7AW2GuMuVxE8oAngWFAFfAlY8zheBSZ0kSounYo1121jO/02U2OK9D1Oqe4a88cau4ZTv7eRszu3XEqMn3VXTKSWbe+Q6Gv7RLavu7jXJFVCUTew2+Sz2Lhv/6JPV/NA8A2Lh574kIG36thP51bgAogJzy9AFhujLlXRBaEp++KcX1p4UShzV39/3Fv7uja4/c29qXvmiqs2trYF5aG3Lm5SGZG+3TjEOGOgrfCX2f+Q3RdeTNcPq7MPgy07YssY/Ngyfl4ioe0j7RjWlrbxsjvxSIKu4gMAS4DfgrcFp49Bzgv/O/FwBs4NOwqNlyBADtuHcvIc0920rmu/zJyozxa6opbXPzg/Od4ZuwZ7fM2by5l3I9dEV9rkI4i3bPfD9wJZJ8yb4AxZj+AMWa/iBR2tqKIzAfmAwTI6GwRRwoai2bTSmvIo99XA+L1IdnZeCYe4cXRHx++NvZDd92QU8MNOSe3c51/FofyBmL7PXild0a+y7CLyOVAjTHmXRE5L9oNGGMWAYsAciQvFb69SQlX77yITS+X0Xe7jWl09rm6O78/O79dhjWmkdvH/k9SarimYCW33f0l/N5mPpu1EfAnpY54imTP/mngcyJyKRAAckTkceCgiBSF9+pFQE08C+1t1mweTtn/XY0JhVKqt1wySE42Uy+s4I+lryethosygmyc8UR4qvcFHSIIuzFmIbAQILxnv90Yc42I/By4Hrg3/Pv5+JWZesTro+FfplE/XvjUp7Ymu5y05CkpZveXimkaavGT/BXJLqfX68n37PcCT4nIPGA3MDc2JaUHV58Axz5/jE0zHsWF0P1e8c7VOjSf+Te8xDf67Ui7IbXTUVRhN8a8QVurO8aYOmB27EtKHy6X0TdpN3iGD6N2VhENo2FcoFpfwwTRHnQq4eqnD+R733+Mcv8BBrj7EI/WdtWRhl0lnO0RhnnqGOLp/rBXLSbI4qMlbG4eFPE6XrG4MncVZ/h93d5uOtOwq7RUHWrhN49cQfFLdV0vHGZl+lj541LemrQkjpWlLg17gm1oPUFlayGew8576d3982BAPscLBX83O640262sbgmwqnkqOVU21qbIvwlxZWezd88YnhuRxTjfwZS8TXY8Oe8dl0TNditzlt7EsCUwancdlsNuhljz+TLOmL+eK7LfosTTvbfe305kc9cD8+j/QQv9Nu4kmlfQbmpmzANN3L/kKg7ceIItZz/WrRrSlYY9gWxs+uzy4nv17V7dB/t0mgcIvxz0OlmuABDdeXOLCVJrtbCuuYzCtcdxrXgvqqADYFvY6zfjf1+wzp/BzumNZLucM0a+hl2lhUUNI3nw8cvI3m3TvzK6PXoHxlD652a+uP0OGs45wfvn/Tb8AdS7adhVygsai1VHSil5ej/Wth4GPUxWvk/+Sgj1OYtjs0L0EbvTsfF7Ew27Smk/qB3PE0tnkb1TKKqriPnzD1jdxPl/uAN7VDOvnfVAr74LjoZdpbQlOyYz+iebsI4ejUs7h7z9PiVvQ9MXprPrzBxKvb33siQNu4ovEVouLqd2ipfMmYfwizfZFXUqa1cTNy79KnlDGnh4wmNM8fe+K9807CquxO1m1xxh/eW/JCAevCkadvPuJso+8GGXj+W1BycwxV/Z9UppRsOu4s9l6OuK/gaXCWUMpqUFV0uIoOmdffV7d/OjUqqd7tlVSlrdEuStpjKa9mVjLCd2QYo9DbtKSVf9fT4jHrAZU1uD1ZycW1f3Nhp2lZoO+ZF3VmHpzSNjRs/ZlXIIDbtSDqFhV8ohNOxKOYQ20CWQX7wM/swedmTPZMAam8wlq9tvLKhUvGnYE8grbl4d8zyMgZG5X6fseTfGYaPVqOTRsCdY+zXTktw6lPPoObtSDqFhV8ohNOxKOYSGXSmH0LAr5RAa9m4yxnD8uI8Pg00ctqK/Kkv8Fu6CfNw5OXGoTqmONOzdZI4fp+RhN1+++3Yu+eC6qNf/P9Nf5MgjGWy/c7wGXiWEhr2bTCiEZ/m79H/oHWo/zI96/Rtyanhr0hLyyw9CLxzcUKUeDbtSDqFhV8ohNOxJdsmgzWy/ZSR1X5uJO79/sstJGWOn7GLP92fSOHc64o3uJpCqcxGFXUT6icgzIrJFRCpEZKaI5InIMhGpDP/OjXexvdHC/pt594ZfMWP+OszAgmSXkzL+POpFVn/jPuq/1IyrT++/6WIiRLpn/zXwqjFmDDAZqAAWAMuNMaOA5eFp5zGGnG0uPltxOfccKiNoohsJ1S0uslwB+rhbU/84y+XGPnsKR66ZgZRPiOumvOImyxXA40nMyLKewYNonDudPRfmMDqwPyHbTLQu314ikgPMAh4GMMa0GmMagDnA4vBii4Er4lNi6hu4+AM814R4/KnZNNotyS4nblw+L9u+4uHBn/ya7XOzQXrPpXuN04Zw5Q9f5Zmv/4I5mYeSXU5cRHKJ63CgFnhERCYD7wK3AAOMMfsBjDH7RaSws5VFZD4wHyBARkyKTjX2sWPYx47hbRqORfcGoyjyHWHFpH70zZiEe8surIYjMa4yeq7MTMy44VgZbW+ToN9N8eA6pvj9SEkT1rlT8dU0Ym3ZDnbne2BjG/wHPNx9aAxnZmznooxgIv+EiFk+YWJgD2N9vfM9CpGF3QNMA24yxqwSkV8TxSG7MWYRsAggR/J0WJbTuK7vBwz698O8XD+Rfd8bifv1dckuCTN6GPysnrlF7wHgxuacjO1AJk9N/y9WTy7lP1ZcxpjvZmAfO9b5k9gWI/5rD2+9UM7vrz2PjV/4TzJc2uCWDJGEvRqoNsasCk8/Q1vYD4pIUXivXgTUxKvIdOE+bnjzeBFjfAcZ7fXhlcjvGZbvzuTK7MMUeP7OwqHjKSgtaX/MHG5I6J7eFQjgGlDA0WFZ3FC0nHl9D5zyaCYAU/x+pvj38eDAo4j7k88GQ3uqYU81gc+ehU10t0TOCrRghg7CU5dB6GDtaY8gouXOzUX6ney5eCLPRUCCQO+8zxtEEHZjzAER2SMiZcaYrcBsYHP453rg3vDv5+NaaRooWrqf+w9eRc00F09f+ysm+aJ/40zzHeOS2/5G5fy2syIbYesfx1P4wNuxLve0WmZNIPTdOmYUrObSzG1AVsK2/XE/K3uWRx46hzc2ljHuh25C1Xt7/qQi7L1hLOO/WNE+68ysSsb6WoEUvwFlD0Q6LNVNwBMi4gN2ADfS1rj3lIjMA3YDc+NTYvoI7agiY0cVhUyn8l8KKXYfIMcVODkUVQRy3Rn8qGATFGwCwDI2ZcNHMzA7G+zwXtGysE+c6HG94vcjno5vgcZBXn47+knO8PtIZtABZgVg1tAVfE1sqvOKcdVnYB8/3u2BOsXvR3w+jo2w+GPp6x97tPcGHSIMuzFmPVDeyUOzY1pNL5Hz/kF+et/VLBgG93/xES7L6H4w3eLi6gtW8HTJ1PZ51pZsRvxqK1Zdfbef15WZye6bJ2POONrhsXEDKhnuCQGpc259bcHfuf0/5lK7ewJjf9OAtfnDqJ/DFQhQfdM0rOlH+eroN2JfZIrTASfjILSjioIHqygsn8B7lw7jsowtPXq+HxVsatvbh31p4GyaHsqCHoRd+gTIPLuW1VOfPs0SqdUqPSsAq6c+zaMj8ln89Bw8m6N/Dgn44awGNk//Y+wLTAMa9jT0zwXrWLhgLu6mQd1+DttvuHnYazGsqmsD17Qyuf+t9BtVzytTHqHQnRn1c0wJVPPjeTZcNqNthi0MXRbEu3Rth2U9pSXsvHowrbltpz+21/DVEW/04C9Ibxr2NHRl9mHm/tODPX6eaNoSYsG7dC0jlwkN187gwEQ3hd1o+J7kC7D1vIfbp4+bVj517DaGLu24bMuw/tx5zTNcm33y24RE/82pRMMeR+76Rn7/xrm8MGIid5c9F9MOJcl+0/68fgSPVk7HXt0P07on8hWNAQO26X7vu1P/dj9efNMOc/Dmszos1zjUZpTvQNJfq1ShYY+j0M5dlC2sRQYP5D8XXcBFo19Jdkkx8//fms3YhZWY4zti8s1Ad3nFzd/LH+HYtI531vGKkOvq3S3s0dCPvHgyBrupCeob2FRRzNerZ7L8eO/otCFBwWpo6FbQ+xwKcdeOL/CD2vHUWE09riXLFaDIk9XhJ9+dGfVevcZq4ge147l536fY0Jq8D7F40LAngFXfwLi7q9lz41Bueu+qZJeTdH3+thnP1zwsu+ccljUPTXY5H7GseSjL7jmHLTeP4+7qy5JdTkzpYXwi2BahvfuQGh8tu6fx6ISTY9YN8x1ipt+KqmttorWYIG+dCLA3eHLIAl999+u1m5qwdzSRnZ/NswenEZBVnNNnf7da52OlxmpixfEinj04jeydTXj21nGkpXedAmjYE8iEgpT97hCP/uVz7fP2XOBn2bU/Z6gnuT3VPsnWoMXND32dondOHtYO37Mfq4e3m3ZVVNF0Zym/LPkyqxe8w88GrO9hpd33y0OfZsW9M8je1YyrYifkZCetlnjRsCeSMVhbt+HeenJWv+KZrG8ppMluu4baJzZDPH784k1SkdBst1JtBbHCLebrTgwjb4v1kSvxYnE5in3sGLyzgdwDQ1lbN5SK3JP9/wvchvw47ukPWU3UWie/EVhbN5TcVfsIVe3GBlxRhN0yNnutZra05uJuTd0LOzXsSZb/1n7u+eF12N62N97xfOEbX/kL3+oXxddZMfazuqk8//C5+I62vXHdLYa8tdXE607yds0hrPvGc33hd9vnhf65nnXlT8Zpi3DRe1/BsySvfTqjJoRds+kT1ji99a0hrnryDvI2GgpXHYzJB2E8aNiTLLSjir47qtqn80ePYN0XS2jO2d5hWa+4Y3puHzRWp8NovXOolMFLdn3kCrN4BR3Abm7G//IaTh09f/v4mTROO4Er3IbsFunR0U6LCbafdtjYHN2ax4g/rPxoHR+ZsGm1vTTbrfjFc9pW/RYTpCqYz6C/hfC/vCZlgw4a9tRTW8/6RyZyRtGkDg9llR/izSmPx2Twh8NWM2ev/hr2B307bmeXIf/Ixh5voyeK/yfIpxpvg/CRdkvpCV497zeM9kZ/aP9hsImL37gJ/87wwJUGit9p/cR17KPHaHq8lKkjb+G8z67nd0NWdljmhaYMbn3lWjJ3uyn+cH9KBx007CnHOnyY/N91fGMB7Lv9LFomh8iIwdVoDbZNYGkO+b/r/Dr56IaYiD3fa2sZekrX/eNXnEnVp/sx2ht9L8SqYD+G/NlDn+ciHxPAbmoid/FK8vx+lpWMh07CvuJYGWUPHcHesCXlgw4a9rSSvzHItFduAW8MotjiZuSW9Ok0klnVyL++OA+7b/Rhdx3xMrrqaPc+wCyLvBV+Sq15HR7y7/Exom5X0j8YIyWmh1+fRCNH8sx00Uvgu00E8cSuld6Egt0eBCLhRBC3G7rTz93YGMvq/t/qcrdtu7PnDcWzNSN6q8xyjpr6Ti880D17OjEGE/zkc81ey5jkBcu2MDEa+y6ZtLusUg6hYVfKITTsSjmEhl0ph9CwK+UQGnalHELDrpRDaNiVcggNu1IOoWFXyiE07Eo5hIZdKYfQsCvlEBp2pRxCw66UQ2jYlXIIDbtSDhFR2EXkOyKySUQ2isifRCQgInkiskxEKsO/c7t+JqVUsnQZdhEZDNwMlBtjJgBu4EpgAbDcGDMKWB6eVkqlqEgP4z1AHxHxABnAPmAOsDj8+GLgiphXp5SKmS7DbozZC/wC2A3sB44YY5YCA4wx+8PL7AcKO1tfROaLyFoRWRukJXaVK6WiEslhfC5te/FSYBCQKSLXRLoBY8wiY0y5Mabc+5Eb/CilEimSw/gLgJ3GmFpjTBBYApwFHBSRIoDw75r4lamU6qlIwr4bmCEiGSIiwGygAngBuD68zPXA8/EpUSkVC13eJMIYs0pEngHW0XYzz/eARUAW8JSIzKPtA2FuPAtVSvWM3v5JqV7kk27/pD3olHIIDbtSDqFhV8ohNOxKOYSGXSmH0LAr5RAadqUcQsOulENo2JVyCA27Ug6hYVfKITTsSjmEhl0ph9CwK+UQGnalHELDrpRDaNiVcggNu1IOoWFXyiE07Eo5hIZdKYfQsCvlEBp2pRxCw66UQ2jYlXIIDbtSDqFhV8ohNOxKOYSGXSmH0LAr5RAadqUcQsOulENo2JVyCA27Ug6hYVfKITTsSjmEhl0phxBjTOI2JlILNAGHErbRnssnfepNp1ohvepNl1pLjDEFnT2Q0LADiMhaY0x5QjfaA+lUbzrVCulVbzrVejp6GK+UQ2jYlXKIZIR9URK22RPpVG861QrpVW861dqphJ+zK6WSQw/jlXIIDbtSDpGwsIvIxSKyVUS2iciCRG03UiJSLCKvi0iFiGwSkVvC8/NEZJmIVIZ/5ya71n8QEbeIvCciL4anU7nWfiLyjIhsCb/GM1O1XhH5Tvg9sFFE/iQigVStNRoJCbuIuIEHgEuAccBVIjIuEduOQgj4rjFmLDAD+Fa4xgXAcmPMKGB5eDpV3AJUnDKdyrX+GnjVGDMGmExb3SlXr4gMBm4Gyo0xEwA3cCUpWGvUjDFx/wFmAq+dMr0QWJiIbfeg5ueBC4GtQFF4XhGwNdm1hWsZQtub7jPAi+F5qVprDrCTcIPwKfNTrl5gMLAHyAM8wIvARalYa7Q/iTqM/8cL+A/V4XkpSUSGAVOBVcAAY8x+gPDvwiSWdqr7gTsB+5R5qVrrcKAWeCR82vGQiGSSgvUaY/YCvwB2A/uBI8aYpaRgrdFKVNilk3kp+Z2fiGQBzwK3GmOOJruezojI5UCNMebdZNcSIQ8wDfitMWYqbddHpORhcPhcfA5QCgwCMkXkmuRWFRuJCns1UHzK9BBgX4K2HTER8dIW9CeMMUvCsw+KSFH48SKgJln1neLTwOdEpAr4b+AzIvI4qVkrtP3/VxtjVoWnn6Et/KlY7wXATmNMrTEmCCwBziI1a41KosK+BhglIqUi4qOtweOFBG07IiIiwMNAhTHmvlMeegG4Pvzv62k7l08qY8xCY8wQY8ww2l7LvxpjriEFawUwxhwA9ohIWXjWbGAzqVnvbmCGiGSE3xOzaWtMTMVao5PAho9LgQ+B7cD3k91Y0Ul9Z9N2arEBWB/+uRToT1tDWGX4d16ya/1Y3edxsoEuZWsFpgBrw6/vc0BuqtYL/AjYAmwEHgP8qVprND/aXVYph9AedEo5hIZdKYfQsCvlEBp2pRxCw66UQ2jYlXIIDbtSDvG/PxI4ytBq66IAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"img = Image.open('bin.png')\n",
"arr = np.asarray(img)\n",
"arr = (arr > 0).astype(np.uint8)\n",
"\n",
"plt.imshow(arr)"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x2a2351d84f0>"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD7CAYAAACscuKmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAagklEQVR4nO3de3yUhZ3v8c9vZpIJScgVcoMAUZCAWhdFAW+1UlutVqwet7a2S6u7rBaL2lYXrbseT0/Py9W21vZoj1i0WG+1rqtUrTfUqlURUFCQOwGSkIQQAgm5z8zv/DFjCJKQSebO83u/XryS5/4l5MtzmWeeEVXFGHP0cyU6gDEmPqzsxjiEld0Yh7CyG+MQVnZjHMLKboxDRFR2ETlfRDaKyBYRWRitUMaY6JPhvs4uIm5gE3AeUAOsAL6lqp9GL54xJlo8ESx7GrBFVbcBiMiTwBxgwLKni1czyIpgk8aYI+mkjW7tkv6mRVL2MUB1n+EaYMbnZxKRecA8gAwymSGzI9ikMeZIluuyAadFcs7e3/8eh50TqOoiVZ2uqtPT8EawOWNMJCIpew1Q3md4LLArsjjGmFiJpOwrgEkiUiEi6cAVwNLoxDLGRNuwz9lV1Sci1wEvA27gIVVdF7VkxpioiuQCHar6IvBilLIYY2LI7qAzxiGs7MY4hJXdGIewshvjEFZ2YxzCym6MQ1jZjXEIK7sxDmFlN8YhrOzGOISV3RiHsLIb4xBWdmMcwspujENY2Y1xCCu7MQ5hZTfGIazsxjiEld0Yh7CyG+MQVnZjHMLKboxDWNmNcQgruzEOYWU3xiGs7MY4hJXdGIewshvjEFZ2YxzCym6MQ1jZjXEIK7sxDmFlN8YhrOzGOISV3RiHGLTsIlIuIm+IyHoRWSci14fGF4jIqyKyOfQ1P/ZxjTHDFc6e3Qf8WFWnADOB+SIyFVgILFPVScCy0LAxJkkNWnZVrVPVD0PftwLrgTHAHGBJaLYlwCUxymiMiYIhnbOLyARgGrAcKFbVOgj+hwAUDbDMPBFZKSIre+iKMK4xZrjCLruIZAP/Bdygqi3hLqeqi1R1uqpOT8M7nIzGmCgIq+wikkaw6I+p6jOh0Q0iUhqaXgrsjk1EY0w0hHM1XoDFwHpV/VWfSUuBuaHv5wLPRT+eMSZaPGHMcwbwXeATEVkdGncrcCfwlIhcDewELo9JQmNMVAxadlV9B5ABJs+ObhxjTKzYHXTGOISV3RiHsLIb4xBWdmMcwspujEOE89KbMY7hHj2aPV+bSE+WUPpaA/5NWwec1zN2DA3nj0NdUPJyDb4d1XFMOnRWdmP6CIwr4qwFy/lq7ifcvu9qco5Q9s7JJcy98UVGujv4Q/UcvFZ2Y1KH60AXf9l0IuuKS2kd58LzjRmMXNuIf/O23nk85WNpPbmMlnEe/lR9Cl0+D/n7uxOYOjx2zm5MH4Gt25l00x74twIyz2rkprv+SM3XSw6Zp+mLY5l751KO++ZGvHfnUzS/A9eqDQlKHD7bsxvTh/p8+Gpqcff0sHffGLZ1F9GdCzL9hN552otc7OwuZGdLPgXV+5L+XP0zoqpx21iOFOgMsTtsTfITjwepnEh3URbbrnBx9cy3e6c98ulplC3xklHfjqzfSqCzM4FJD7Vcl9Gie/u9vd327Mb0Q30+dO0G0rxe+MeTuDhnNSVuP0XuLP6aN5XMlXX4GxuJ364ycnbObswRaHc3kx7p4l/+/UbOfu8aetSf6EjDZnv2FCdp6eAa6E2Jffj9qM8X+0BHG1Xk3TXkvQv7jptF4IxAohMNm5U9hXkqxrPpmjL8pYM/2y9th5eJ91fhq6uPQzKTjKzsKcxXnMu1F77Mjwq2DTrvDXXT2fTksWBldywrewrovOg06me4DxuvHvjt21/mt+7BLxNJjwvPlS7EP+vQ8Qplb3WT9tqqqOU1ycnKnuxEqP6qsP7S3xw26b7myfz1mrNxv7920NV0fmUa19zzFN/IPvS5oD3q5xTXjUx4LWqJTZKysicZ18iRtJ07hY7C4J5cBQgosz/5Ji45dA9evauAKU0t+MO48JbR0M5tq+Zwf+H+Q8YHVPB7oenqWcgABwjuLqXwndqUuXnE9M/KnmSkZDQFP9rBXROCT+z2q3Dxn35MzndbD5t3im8f/v1hPsJ/zUaO+2EOuA49HRC3i9rb/Tz+H/cMuOj7HeN58JbLyLSypzQre5Jw5+XiO76CfWMz6Gndy+K9ZwDBPW9mveBvbIxo/erz4W/ae/gEEUbUHtu7vf5Ut+fTWu4m/dxTyNhUj6+mNqIsJjGs7Emi58RjGHv3FgrS2njn3hl8vLoyOCEAYxq3ELNbOVSpWLyNj5dWDpytMJMRNzdw+TXv8+C9FzPqASt7KrKyJ5hr5EikZDStY7wUpLXhEiW7tpvAx/F7F5Wvrv6IL8l5x5TR1J1DnrudzkLBPXki7N0f8dGGiS8re4Id+PJUCm/YTlfLXt65dwbZtd14V1fFbk8+DP49TYy+czT3F19O+9e7mfBoDW//+WTK7t4DcXwjlYmMlT2WRHBlZyPugd+C0Fbs4mdjX2VJ4xk0rBL86zYmVdEBtKsLeW8NWZmZyKWVLCx+jZfKT8JdkI92dBJob090RBMGK3sMeSrGs/7HxeSO3T/gPPsafFz/22vIqg+Qt2tjHNMNXaCzi/GPuZiz/GZcxwSoXVxMz4p8xt+z2gqfAqzsMeQvyGbumW9z26i1uCW4d/froW+kOPPjy8n/aQv+xsak26MfJuAn/eWVFL0MB/5zFh+c+ggz5btIehpY15OelT2G3HV7efrRc1gy/iwWzv4L38+p5qINc9i2qrx3ntxNoG11CUw5PGPe8nECC8jbAIGOqkTHMWGwsseQr3YXZXfvwjOmjGemTuPbI7dR/fo4jvn5u4fMl4pvmvS+sIJjXgh+b5foUoM9vCJOVMN4z7kxMWRlN8YhrOzGOISV3RiHsAt0MeTKyIDjJtBRkk1jSzv3NZ+ItznRqYxThb1nFxG3iHwkIs+HhgtE5FUR2Rz6mh+7mKlJjh1P850+Tvg/a3C9mceyfz2dsue2JzqWcaihHMZfD6zvM7wQWKaqk4BloWHTR2BEGheNXctVhe+QsTf4lFJf7a5ExzIOFVbZRWQscCHw+z6j5wBLQt8vAS6JajJjTFSFe87+a+BmYGSfccWqWgegqnUiUtTfgiIyD5gHkEHm8JMaEyJeLyLDu28h0N0DgaS/MTkmBi27iFwE7FbVVSJyzlA3oKqLgEUQ/Ky3oS5vTF+e8rFsml+Olg/989X83S4mPCmkv7wyBsmSXzh79jOAi0Xka0AGkCMijwINIlIa2quXAruPuBZjoiBQmMO3L3iLO0avG/Kydb4DXPDRzRS/HINgKWBIn+Ia2rP/RFUvEpG7gSZVvVNEFgIFqnrzkZZ3yqe4uqdMYsc3RtOdo/izA0i3MOmJNnTFJ4mOlrLckyey47Ki3p+peoZxkKjgbnPh7hJK/+7D++KKwZdxuTlw2XQap7nwZyiBzAC56zyUPrSGQFvb0DPEWKw+xfVO4CkRuRrYCVwewbqOKq1TCrjnqgfZ3F3Cf//gPNxvrUEdep4YLQcqC7jrqoeo78nj8esuxPPG6iGvw52bw87FZXxw2sN8wbuAY18cfBlJ81B3oY+N5y3iuBeuoXLBJwS6ugik4BN6hlR2VX0TeDP0fRNw9O+mh0EF0sVPmvgRVcdeEIoG99Tj2D2rkPYS4eY1l9HekkFlQxuBYfxMtaODwKpczvFeibpgz7/OGnSZgEf4QsVW0iT0HH+fL2UfxWV30JmkVv/FQn570338suartF9fhKtqC/6WA8NaV6Czk/G/XoM8kEHTbYU8devdYb32PNrtIXi5KrVZ2U3ycLlxnXgcnaXZvaM6RwkP7T6bNTvGUrmnEd++gR/xFY5AWxu0t5NZ6+Lexi8d9ik7n+cmwKX5Kzkj9btuZTfJw5XhZcMPsvnFl/7UO27hykupvW4ClU178UfrE2hVKX9kC1tePGbwWdPcvPwfU1g787HobDuBrOwm4cTjwT22DH/hSNzZh35una8tDdm4FV/r4R9/FQl/w25oGPzVYvF6aT9wfFS3nShWdpNw7uIi1t8xitmVG9nxwkn83z9+s3daZX0bgTZ7mmU0WNlN4nnTmTmxip+XvcJ5rSeR/tLB179T8fl8ycoeXmGMQ1jZjRlMQOjR0Ov6bjcM8004iWaH8cYcgfb4KH3Bw/G7rsOVoWy896Skvl32SGzPbsyRBPxk/3k5Fbe+h6dN2PD1+8i6sB7JTL23a9ue3SSMZ0wZdRePp3OU0LF2NDM3L2Dchp5ExzpqWdlNwnRXFPGt+a8w2VvHf/70n8hZ+knw4RImJqzsJnHcQq67gxxXJ+JXAp1DfyCFCZ+dsxvjEFZ2YxzCym6MQ1jZjXEIK7sxDmFlN8YhrOzGOISV3RiHsLIb4xBWdmMcwspujENY2Y1xCCu7MQ5hZTfGIazsxjiEld0Yh7Cyx4DLp2ztLqK2O5+eTA/uwgIkLT3RsYzD2ZNqYmDkh7t48OeX0F7konPuATxXl1F6vxfP66sSHc04mO3ZY8BXXUPuo+9T9lYL/zT1Ax4/ZTEHxtie3SSWld0Yh7CyG+MQVnZjHCKssotInog8LSIbRGS9iMwSkQIReVVENoe+5sc6rDFm+MLds98LvKSqlcBJwHpgIbBMVScBy0LDpg/X/nYe+vAM5m/8FgfKhf1XzsR9/ORExzIONWjZRSQHOBtYDKCq3aq6D5gDLAnNtgS4JDYRU5d/6w4qb9pJzr95yTu7nrt+9v+ovqAw0bGMQ4XzOvsxQCPwsIicBKwCrgeKVbUOQFXrRKSov4VFZB4wDyCD1PswvIgE/PgbG/Gkp5Gd5uLk9E783kSHioAI7uOOpadoJOk1e/FV7Yhode62Hh6vPo2KnCbait2MOPMfSN9aj6+uPkqBTV/hHMZ7gJOB36nqNKCNIRyyq+oiVZ2uqtPTSOXfdOPyetkwv5DLHniFbXPHRPw55fLpNkb8KIOqOyrxXrSbL93/HrsvqIhSWvN54ZS9BqhR1eWh4acJlr9BREoBQl93xybiUcDvZ1dLDu905qIe8FSMx11YkOhUQ+dykVbcwTV5tXQV+iNeXaC9ncDaDWRtbOSY3Cb+Of9DunMj+w8kFtyjCvFUjAcX/K0jkz37syEQ+d8/3gYtu6rWA9Ui8tmVpdnAp8BSYG5o3FzguZgkPAr4m5op+WU6P7vle3SW+sh9tIUd8yrB5U50NDMISUtn+7WTyX20hZ68AP9r4VWMv1cI7G9JdLQhC/fe+B8Cj4lIOrAN+D7B/yieEpGrgZ3A5bGJmPq0pxvXO6sZ6fXC+Sfym3HPc+qEibhzsgl0dKJdXYmOOChXRgaSlYUGhN3+NsSXfHvgmHAJnZM6ebziDSrWTST72Y9Qnw9NdK5hCKvsqroamN7PpNlRTXOU0x4fE/4szF57E+4yperBccjqkYz/zScEWlsTHW9A7sICtt44GTnuALotk/M+uImK1Z2gqfgr71x2B108Bfykv7SCknvexd0Jq09/iPyz6pGMjEQnOyLJzuLEL27mvVkPkLFHgvnf+DAm29IkO2CQCC9CJhN7i2uCXVG+inv+5/m4Oo8FIK3VRcUTu/Fv3JLgZPGjzftZ9/QUZo2fDGMDbP3lTMreDjDi2Q/ilsF9/GS2/WMhvuzAodnccNkJywdYKrVY2RPsh/k7+MElv+sdfqE9m1+tuBLvxgSGijN/czMlv34P98iRVC0ez8en/4FK13wmPidxO1VonZzHfd95gHMyeg6b5paj4wDYyp5gS9syuWPDRXT7gv8UXZ1peE5OI/3Y03vnKVrZhry7Jm6Z3BMraJhdgt8bPITtzoF5BS/iFQ89p7XSsOB0Rq/pwPW3j6K3UVVUFRElTdyUH19P/YJZSKjrmQ0Bcl9cF7VrG4GzptE4bUTvcFc+XLviSrxe34DLSKeLuvmnkbvDR+Zf16TEhdW+rOwJ9quqr1CyoAttCt41Fji2nFH3fcrdY58HoAc475GbmPBu/DLtO6WY23+yhFO9wVsn3CLkuzJIk3SWn76I9pl+znz6J0x8K3Z73uenPklr5cHiXVt1Gd2rRkE0yi5C1cUZvPnNu/jsxc9rqy6ja34+snPXgIttvfl4nv3RXVy6+p/Jficbv5XdhCOrFv5l52x2ri+hsnk9/pbg67aePft5d8sx3OY6H4AAgt8LHXNO6102Y083rpXrI9qziMeDnjyFjtIRh01rLXfx+11n85z3wIDLZ9bF9tD29Y4C/rJ3GoHQFbstTaPImpVD1oQCMlZtw9/cHPa6XBkZ+KZX0lV48GlB6lFu33V+7/CaLeVMbd6Fr+UIr58rlHm8ZHq7IQUP7a3sCVLypw3sea2QKW3b8fW5QcNfV8+UW5W6EaG3GogQ+GGA2375cO88t6z/BiXXjsJXUzvs7btGjmTT9R7un/HwYdNu33wxnQuLqdsz8HsZyps34o/h+fT1f/s2U+5q7j1yyJqVw9dveoNdXXlsuaES+fsQyl5aTPttLdw+8S8ABHBx3XPfp+67B9/OMbWjDl9dQ3T/EknGyp4g/uZm6GfvpD7foSUWIW1fEfW+3N5Rae4AnZNLSB+Ve9jy4erJTseT5jtkvZ/Z25JFQVU9vvo4//L7/XQ0ZvKHliIyqtPxb952sOwVBezqyqPL76F1XAZ5bVPDXm17cSZuaer9uwZwkd7sCq5/GEaPaOPA8eWk7RsNgPT40e01BNrahrW+eBGN440ROVKgM8Tuwxkq98QKussPPhtkzwkZfOHba6nIbBr2Opu6s3nz6VMo+aDzsGlpzZ3oui1oT/ew1z8sIrinTKK7OBvvjr34tm3vneTOz8c3dTyt4zJI+14DXyrZFPZqNx8oYsOjlRRsOHja493eNOR37W3/2Sw+vuo3bOvp4b9bptGlwX3lquZxdP17Ca63o3jBcpiW6zJadG+/NwfYnj0F+LdU4d5S1Ts8Mm8G5+ZvYEbG9mGvc5uvgOW7T+735phAP/PHhSr+Tzfh/hQ+f03c39yM/L2ZvNZKSnL3cEXuirBX+3raZHbtnHjI33Xga+6DK/e4+B+5H+IPXU/I97TxXOYYkv35wVb2FJSzspYH7riU+73Dv7vL1aOMXl5Pqr13S3bWsfXOKczNOyHsZTwdSv7qnREVvK9b68/i7YdOJa1dD67/4+1RW3+sWNlTkK+6hpwnaiJeT6oVHcC/bz8jnv2Aw19DOLKIiuhyI+7gi3TtgR7eb5hA6VMb8e85eBqV7EUHK7sxRyQeD3u+dypNs3pwtSqnL/4JuVsC6IHqREcbstR7sdCYeHK72Xt2F1UX/J5App/x//sDch99n0Dn4Rc2k53t2Y05Er+fgre9VPivZvS7HtSfiic/QVZ2Y45AfT4KF3/AqD+4g0VP4ffwW9mNGUzAj6bgM+c+z87ZjXEIK7sxDmFlN8YhrOzGOISV3RiHsLIb4xBWdmMcwspujENY2Y1xCCu7MQ5hZTfGIazsxjiEld0Yh7CyG+MQVnZjHMLKboxDWNmNcYiwyi4iN4rIOhFZKyJPiEiGiBSIyKsisjn0NX/wNRljEmXQsovIGGABMF1VTwDcwBXAQmCZqk4CloWGjTFJKtzDeA8wQkQ8QCawC5gDLAlNXwJcEvV0xpioGbTsqloL/ALYCdQB+1X1FaBYVetC89QBRf0tLyLzRGSliKzsIbU+vN6Yo0k4h/H5BPfiFUAZkCUi3wl3A6q6SFWnq+r0NLzDT2qMiUg4h/FfBqpUtVFVe4BngNOBBhEpBQh93R27mMaYSIVT9p3ATBHJFBEBZgPrgaXA3NA8c4HnYhPRGBMNg35IhKouF5GngQ8JfljlR8AiIBt4SkSuJvgfwuWxDGqMiYxoHD/OJkcKdIbMjtv2jHGa5bqMFt0r/U2zO+iMcQgruzEOYWU3xiGs7MY4hJXdGIewshvjEFZ2YxzCym6MQ1jZjXEIK7sxDmFlN8YhrOzGOISV3RiHsLIb4xBWdmMcwspujENY2Y1xCCu7MQ5hZTfGIazsxjiEld0Yh7CyG+MQVnZjHMLKboxDWNmNcQgruzEOYWU3xiGs7MY4hJXdGIewshvjEFZ2YxzCym6MQ1jZjXEIK7sxDmFlN8YhrOzGOISV3RiHEFWN38ZEGoE2YE/cNhq5UaRO3lTKCqmVN1WyjlfV0f1NiGvZAURkpapOj+tGI5BKeVMpK6RW3lTKOhA7jDfGIazsxjhEIsq+KAHbjEQq5U2lrJBaeVMpa7/ifs5ujEkMO4w3xiGs7MY4RNzKLiLni8hGEdkiIgvjtd1wiUi5iLwhIutFZJ2IXB8aXyAir4rI5tDX/ERn/YyIuEXkIxF5PjSczFnzRORpEdkQ+hnPSta8InJj6HdgrYg8ISIZyZp1KOJSdhFxA/cBFwBTgW+JyNR4bHsIfMCPVXUKMBOYH8q4EFimqpOAZaHhZHE9sL7PcDJnvRd4SVUrgZMI5k66vCIyBlgATFfVEwA3cAVJmHXIVDXmf4BZwMt9hm8BbonHtiPI/BxwHrARKA2NKwU2JjpbKMtYgr905wLPh8Yla9YcoIrQBeE+45MuLzAGqAYKAA/wPPCVZMw61D/xOoz/7Af4mZrQuKQkIhOAacByoFhV6wBCX4sSGK2vXwM3A4E+45I16zFAI/Bw6LTj9yKSRRLmVdVa4BfATqAO2K+qr5CEWYcqXmWXfsYl5Wt+IpIN/Bdwg6q2JDpPf0TkImC3qq5KdJYweYCTgd+p6jSC749IysPg0Ln4HKACKAOyROQ7iU0VHfEqew1Q3md4LLArTtsOm4ikESz6Y6r6TGh0g4iUhqaXArsTla+PM4CLRWQ78CRwrog8SnJmheC/f42qLg8NP02w/MmY98tAlao2qmoP8AxwOsmZdUjiVfYVwCQRqRCRdIIXPJbGadthEREBFgPrVfVXfSYtBeaGvp9L8Fw+oVT1FlUdq6oTCP4sX1fV75CEWQFUtR6oFpHJoVGzgU9Jzrw7gZkikhn6nZhN8GJiMmYdmjhe+PgasAnYCvw00Rcr+sl3JsFTi4+B1aE/XwMKCV4I2xz6WpDorJ/LfQ4HL9AlbVbgH4CVoZ/vs0B+suYF7gA2AGuBPwLeZM06lD92u6wxDmF30BnjEFZ2YxzCym6MQ1jZjXEIK7sxDmFlN8YhrOzGOMT/B7nUcJawRKIuAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import scipy.ndimage.morphology as m\n",
"\n",
"def skeletonize(img):\n",
" h1 = np.array([[0, 0, 0],[0, 1, 0],[1, 1, 1]]) \n",
" m1 = np.array([[1, 1, 1],[0, 0, 0],[0, 0, 0]]) \n",
" h2 = np.array([[0, 0, 0],[1, 1, 0],[0, 1, 0]]) \n",
" m2 = np.array([[0, 1, 1],[0, 0, 1],[0, 0, 0]]) \n",
" hit_list = [] \n",
" miss_list = []\n",
" for k in range(4): \n",
" hit_list.append(np.rot90(h1, k))\n",
" hit_list.append(np.rot90(h2, k))\n",
" miss_list.append(np.rot90(m1, k))\n",
" miss_list.append(np.rot90(m2, k)) \n",
" img = img.copy()\n",
" while True:\n",
" last = img\n",
" for hit, miss in zip(hit_list, miss_list): \n",
" hm = m.binary_hit_or_miss(img, hit, miss) \n",
" img = np.logical_and(img, np.logical_not(hm)) \n",
" if np.all(img == last): \n",
" break\n",
" return img\n",
"\n",
"plt.imshow(skeletonize(arr))"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [],
"source": [
"from numba.core import config"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {},
"outputs": [],
"source": [
"config.DEBUG_ARRAY_OPT = 0"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Wall time: 5.67 s\n"
]
}
],
"source": [
"%%time\n",
"from numba.core import typing\n",
"args = (arr,)\n",
"result = None\n",
"array_types = tuple([typing.typeof.typeof(x) for x in args])\n",
"array_types_full = array_types\n",
"\n",
"(real_ret, typemap, calltypes) = kernel.get_return_type(array_types)\n",
"new_func = kernel._stencil_wrapper(result, None, real_ret, typemap,\n",
" calltypes, *array_types_full)\n"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"67.6 碌s 卤 1.26 碌s per loop (mean 卤 std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit new_func.entry_point(arr)"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x2a236f78970>"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAq4AAAKsCAYAAAAtNz8NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX9klEQVR4nO3dX6jn913n8dd7M/ljUoLN7jbkT3dTIailYCqD7dpFlo1LXC0mN3ErVILbJTfuWsVFond7sdALEXuxCCFVAhY11ELCIsYyKuzehMY2oG22ttRums6YZLdqpLBJip+9OD/I1J4zc2bm/HvN7/GAcM7v3/w+c+Y7Z5588v2+z6y1AgAAJ90/Ou4FAADAfghXAAAqCFcAACoIVwAAKghXAAAqCFcAACqcupIXz8yPJPlokmuSPLbW+siFnn/dXL9uyE1X8pYAAFzF/l++kdfXa7PbY3O5c1xn5pokf5Hk3yR5Mcmnk/zkWuvze73m5rllvWfuvaz3AwDg6vfMOpNX19d3DdcrOVXgB5J8aa315bXW60l+J8n9V/DrAQDAnq4kXO9I8tXzbr+4uQ8AAA7clZzjutsW7reddzAzDyd5OEluyI1X8HYAAGyzK9lxfTHJ28+7fWeSs//wSWutR9dap9dap6/N9VfwdgAAbLMrCddPJ7l7Zt4xM9cl+UCSpw5mWQAA8K0u+1SBtdY3Z+Y/Jnk6O+OwfmOt9bkDWxkAAJzniua4rrV+P8nvH9BaAABgT35yFgAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVLhquM/P2mfnjmXl+Zj43Mx/e3H/LzHxqZr64+fjWw18uAADbaj87rt9M8gtrre9N8t4kPzMz70zySJIza627k5zZ3AYAgENx0XBda51ba31m8/nfJXk+yR1J7k/y+OZpjyd54JDWCAAAl3aO68zcleTdSZ5Jcuta61yyE7dJ3rbHax6emWdn5tk38toVLhcAgG2173Cdmbck+b0kP7fWenW/r1trPbrWOr3WOn1trr+cNQIAwP7CdWauzU60fnyt9cnN3S/NzG2bx29L8vLhLBEAAPY3VWCSfCzJ82utXz3voaeSPLT5/KEkTx788gAAYMepfTznfUl+Ksmfzcxzm/t+OclHkjwxMx9K8kKSBw9lhQAAkH2E61rrfyaZPR6+92CXAwAAu/OTswAAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqHDquBcAAHAxT599btf777v9niNdB8fLjisAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVTh33AoDt9fTZ5w7017vv9nsO9NcD4GSx4woAQAXhCgBABeEKAEAF4QoAQAXhCgBABVMFgEtykJMAjmoKgOkFAFcHO64AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUMA4LrnJGQR38mhtHggFcDey4AgBQQbgCAFBBuAIAUEG4AgBQQbgCAFDBVAG4ClzoKndXrR88X1OA42HHFQCACsIVAIAKwhUAgArCFQCACsIVAIAKwhUAgArGYcFVwHgmALaBHVcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqGIcFJ8zTZ5/b8zFjrwDYZnZcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKhgqgAck72mB5gcAAC7s+MKAEAF4QoAQAXhCgBABeEKAEAF4QoAQAXhCgBABeOw4JgYewUAl8aOKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVTBWAQ/T02ef2fMxUAQC4NHZcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKggXAEAqGAcFgBwIhghyMXYcQUAoIJwBQCggnAFAKCCcAUAoIJwBQCgwr6nCszMNUmeTfK1tdb7Z+aWJL+b5K4kX0nyE2utvz6MRcJJt9eVsK6CBYCDcyk7rh9O8vx5tx9JcmatdXeSM5vbAABwKPYVrjNzZ5IfS/LYeXffn+TxzeePJ3ngQFcGAADn2e+O668l+cUkf3/efbeutc4lyebj2w52aQAA8KaLhuvMvD/Jy2utP72cN5iZh2fm2Zl59o28djm/BAAA7OvirPcl+fGZ+dEkNyS5eWZ+K8lLM3PbWuvczNyW5OXdXrzWejTJo0ly89yyDmjdAABsmYvuuK61fmmtdeda664kH0jyR2utDyZ5KslDm6c9lOTJQ1slAABbb9/jsHbxkSRPzMyHkryQ5MGDWRIAdNtrRN5JYEwfzS4pXNdaf5LkTzaf/98k9x78kgAA4Nv5yVkAAFQQrgAAVBCuAABUEK4AAFS4kqkCAHDVu5wJASf5yv2DnnhwOb/XvdZwkr9unAx2XAEAqCBcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKhgHBbs04VGyBjhAlevq+3v90H/fq62cWGcbHZcAQCoIFwBAKggXAEAqCBcAQCoIFwBAKhgqgAAcNkuZ0LAXpMITBvgYuy4AgBQQbgCAFBBuAIAUEG4AgBQQbgCAFBBuAIAUME4LACuKnuNWkqMW4J2dlwBAKggXAEAqCBcAQCoIFwBAKggXAEAqGCqAACV9poeYHIAXL3suAIAUEG4AgBQQbgCAFBBuAIAUEG4AgBQQbgCAFDBOCwA4EjtNbJsrxFnF3oN28WOKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVTBWAfXJFKwAcLzuuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVDAOC4AT6+mzz+35mBF1sH3suAIAUEG4AgBQQbgCAFBBuAIAUEG4AgBQwVQB2CdXNwMcrgt9L93re7Dvv9vFjisAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABVOHfcCAGAv991+z3EvAThB7LgCAFBBuAIAUEG4AgBQQbgCAFBBuAIAUEG4AgBQQbgCAFBBuAIAUEG4AgBQQbgCAFBBuAIAUEG4AgBQ4dRxLwBa3Hf7PXs+9vTZ5y75NQDApbHjCgBABeEKAEAF4QoAQAXhCgBABeEKAEAF4QoAQAXjsAA4sfYaNZcYNwfbyI4rAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFYQrAAAVhCsAABWEKwAAFU4d9wLganDf7ffsev/TZ5+75NcAbCvfM7kYO64AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUOLWfJ83MdyZ5LMm7kqwk/z7JF5L8bpK7knwlyU+stf76MBYJre67/Z7jXgJUu9DfoafPPnfJrwG67XfH9aNJ/mCt9T1Jvi/J80keSXJmrXV3kjOb2wAAcCguGq4zc3OSH0rysSRZa72+1vqbJPcneXzztMeTPHA4SwQAgP3tuH5XkleS/ObMfHZmHpuZm5LcutY6lySbj2/b7cUz8/DMPDszz76R1w5s4QAAbJf9hOupJN+f5NfXWu9O8o1cwmkBa61H11qn11qnr831l7lMAAC23X7C9cUkL661ntnc/kR2QvalmbktSTYfXz6cJQIAwD6mCqy1/mpmvjoz373W+kKSe5N8fvPfQ0k+svn45KGuFK4yrogGtpXvf1yufY3DSvKfknx8Zq5L8uUkP52d3donZuZDSV5I8uDhLBEAAPYZrmut55Kc3uWhew90NQAAsAc/OQsAgArCFQCACsIVAIAKwhUAgAr7nSoAHLC9xr7sNSbmQq/h6uRYAPhWdlwBAKggXAEAqCBcAQCoIFwBAKggXAEAqGCqAJwwF7pafK+rzF1hfnXy53phJnPA9rHjCgBABeEKAEAF4QoAQAXhCgBABeEKAEAF4QoAQAXjsKDI5Yz/uZxfj6NjdNPB83U7eJf7PQYOmh1XAAAqCFcAACoIVwAAKghXAAAqCFcAACqYKgBXgcu9ivogrxTepiu5L+frttfXZ5u+bifZQf6ZnnQn+fdqygYXY8cVAIAKwhUAgArCFQCACsIVAIAKwhUAgArCFQCACrPWOrI3u3luWe+Ze4/s/YCjc5CjtZLLG31z0GvYi7E8JEd7zBtdxzZ5Zp3Jq+vrs9tjdlwBAKggXAEAqCBcAQCoIFwBAKggXAEAqGCqAAAAJ4apAgAA1BOuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFQQrgAAVBCuAABUEK4AAFTYV7jOzM/PzOdm5s9n5rdn5oaZuWVmPjUzX9x8fOthLxYAgO110XCdmTuS/GyS02utdyW5JskHkjyS5Mxa6+4kZza3AQDgUOz3VIFTSb5jZk4luTHJ2ST3J3l88/jjSR448NUBAMDGRcN1rfW1JL+S5IUk55L87VrrD5PcutY6t3nOuSRvO8yFAgCw3fZzqsBbs7O7+o4ktye5aWY+uN83mJmHZ+bZmXn2jbx2+SsFAGCr7edUgR9O8pdrrVfWWm8k+WSSH0zy0szcliSbjy/v9uK11qNrrdNrrdPX5vqDWjcAAFtmP+H6QpL3zsyNMzNJ7k3yfJKnkjy0ec5DSZ48nCUCAMDORVcXtNZ6ZmY+keQzSb6Z5LNJHk3yliRPzMyHshO3Dx7mQgEA2G6z1jqyN7t5blnvmXuP7P0AAOjyzDqTV9fXZ7fH/OQsAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACoIVwAAKghXAAAqCFcAACrMWuvo3mzmlST/e3PznyT5P0f25pxkjgUSxwFvciyQOA622T9fa/3T3R440nD9ljeeeXatdfpY3pwTxbFA4jjgTY4FEscBu3OqAAAAFYQrAAAVjjNcHz3G9+ZkcSyQOA54k2OBxHHALo7tHFcAALgUThUAAKDCsYTrzPzIzHxhZr40M48cxxo4ejPz9pn545l5fmY+NzMf3tx/y8x8ama+uPn41uNeK4dvZq6Zmc/OzH/f3HYcbKGZ+c6Z+cTM/K/N94Z/4VjYPjPz85t/F/58Zn57Zm5wHLCbIw/XmbkmyX9L8m+TvDPJT87MO496HRyLbyb5hbXW9yZ5b5Kf2fzZP5LkzFrr7iRnNre5+n04yfPn3XYcbKePJvmDtdb3JPm+7BwTjoUtMjN3JPnZJKfXWu9Kck2SD8RxwC6OY8f1B5J8aa315bXW60l+J8n9x7AOjtha69xa6zObz/8uO/9A3ZGdP//HN097PMkDx7JAjszM3Jnkx5I8dt7djoMtMzM3J/mhJB9LkrXW62utv4ljYRudSvIdM3MqyY1JzsZxwC6OI1zvSPLV826/uLmPLTIzdyV5d5Jnkty61jqX7MRtkrcd49I4Gr+W5BeT/P159zkOts93JXklyW9uTht5bGZuimNhq6y1vpbkV5K8kORckr9da/1hHAfs4jjCdXa5z2iDLTIzb0nye0l+bq316nGvh6M1M+9P8vJa60+Pey0cu1NJvj/Jr6+13p3kG/G/g7fO5tzV+5O8I8ntSW6amQ8e76o4qY4jXF9M8vbzbt+Znf8lwBaYmWuzE60fX2t9cnP3SzNz2+bx25K8fFzr40i8L8mPz8xXsnOq0L+emd+K42AbvZjkxbXWM5vbn8hOyDoWtssPJ/nLtdYra603knwyyQ/GccAujiNcP53k7pl5x8xcl50TsJ86hnVwxGZmsnMu2/NrrV8976Gnkjy0+fyhJE8e9do4OmutX1pr3bnWuis7f///aK31wTgOts5a66+SfHVmvntz171JPh/HwrZ5Icl7Z+bGzb8T92bnGgjHAd/mWH4Awcz8aHbOcbsmyW+stf7rkS+CIzcz/zLJ/0jyZ3nz3MZfzs55rk8k+WfZ+Qb24Frr68eySI7UzPyrJP95rfX+mfnHcRxsnZm5JzsX6V2X5MtJfjo7myqOhS0yM/8lyb/LzvSZzyb5D0neEscB/4CfnAUAQAU/OQsAgArCFQCACsIVAIAKwhUAgArCFQCACsIVAIAKwhUAgArCFQCACv8fWoekdhhS3IQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 864x864 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"arr2 = arr.copy()\n",
"for i in range(10):\n",
" hm = new_func.entry_point(arr2)\n",
" arr2 = np.logical_and(arr2, np.logical_not(hm)) \n",
"# arr2 = (arr - new_func.entry_point(arr)).astype(np.uint8)\n",
"# arr3 = (arr2 - new_func.entry_point(arr2)).astype(np.uint8)\n",
"# arr4 = (arr3 - new_func.entry_point(arr3)).astype(np.uint8)\n",
"# arr5 = (arr4 - new_func.entry_point(arr4)).astype(np.uint8)\n",
"# arr6 = (arr5 - new_func.entry_point(arr5)).astype(np.uint8)\n",
"plt.figure(figsize=(12, 12))\n",
"plt.imshow(arr2, interpolation='none')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"h1 = np.array([[0, 0, 0],\n",
" [0, 1, 0],\n",
" [1, 1, 1]]) \n",
"\n",
"m1 = np.array([[1, 1, 1],\n",
" [0, 0, 0],\n",
" [0, 0, 0]]) \n",
"\n",
"h2 = np.array([[0, 0, 0],\n",
" [1, 1, 0],\n",
" [0, 1, 0]]) \n",
"\n",
"m2 = np.array([[0, 1, 1],\n",
" [0, 0, 1],\n",
" [0, 0, 0]]) "
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"hm1 = np.array([[2, 2, 2],\n",
" [0, 1, 0],\n",
" [1, 1, 1]]) \n",
"\n",
"hm2 = np.array([[0, 2, 2],\n",
" [1, 1, 2],\n",
" [0, 1, 0]]) \n",
"\n",
"hm_list = []\n",
"\n",
"for k in range(4): \n",
" hm_list.append(np.rot90(hm1, k))\n",
" hm_list.append(np.rot90(hm2, k))"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"a00 = a[-1, -1]\n",
"a01 = a[-1, 0]\n",
"a02 = a[-1, 1]\n",
"a10 = a[0, -1]\n",
"a11 = a[0, 0]\n",
"a12 = a[0, 1]\n",
"a20 = a[1, -1]\n",
"a21 = a[1, 0]\n",
"a22 = a[1, 1]\n",
"if a00 == 0 and a01 == 0 and a02 == 0 and a11 == 1 and a20 == 1 and a21 == 1 and a22 == 1: return 1\n",
"if a01 == 0 and a02 == 0 and a10 == 1 and a11 == 1 and a12 == 0 and a21 == 1: return 1\n",
"if a00 == 0 and a02 == 1 and a10 == 0 and a11 == 1 and a12 == 1 and a20 == 0 and a22 == 1: return 1\n",
"if a00 == 0 and a01 == 0 and a10 == 0 and a11 == 1 and a12 == 1 and a21 == 1: return 1\n",
"if a00 == 1 and a01 == 1 and a02 == 1 and a11 == 1 and a20 == 0 and a21 == 0 and a22 == 0: return 1\n",
"if a01 == 1 and a10 == 0 and a11 == 1 and a12 == 1 and a20 == 0 and a21 == 0: return 1\n",
"if a00 == 1 and a02 == 0 and a10 == 1 and a11 == 1 and a12 == 0 and a20 == 1 and a22 == 0: return 1\n",
"if a01 == 1 and a10 == 1 and a11 == 1 and a12 == 0 and a21 == 0 and a22 == 0: return 1\n"
]
}
],
"source": [
"for i, j in np.ndindex(3, 3):\n",
" print(f'a{i}{j} = a[{i-1}, {j-1}]')\n",
" \n",
"for hm in hm_list:\n",
" and_terms = []\n",
"\n",
" for (i, j), v in np.ndenumerate(hm):\n",
" i2 = i - 1\n",
" j2 = j - 1\n",
" if v == 0:\n",
" continue\n",
" elif v == 1:\n",
" term = f'a{i}{j} == 1'\n",
" elif v == 2:\n",
" term = f'a{i}{j} == 0'\n",
" and_terms.append(term)\n",
"\n",
" cond = ' and '.join(and_terms)\n",
" print(f'if {cond}: return 1')"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['a[-1, -1] == 0',\n",
" 'a[-1, 0] == 0',\n",
" 'a[-1, 1] == 0',\n",
" 'a[0, 0] == 1',\n",
" 'a[1, -1] == 1',\n",
" 'a[1, 0] == 1',\n",
" 'a[1, 1] == 1']"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"and_terms"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'ab'"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@nb.njit\n",
"def f(a, b):\n",
" return a + b\n",
"\n",
"f(\"a\", \"b\")"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"True\n",
"@overload(operator.add)\n",
"@overload(operator.iadd)\n",
"def unicode_concat(a, b):\n",
" if isinstance(a, types.UnicodeType) and isinstance(b, types.UnicodeType):\n",
" def concat_impl(a, b):\n",
" new_length = a._length + b._length\n",
" new_kind = _pick_kind(a._kind, b._kind)\n",
" new_ascii = _pick_ascii(a._is_ascii, b._is_ascii)\n",
" result = _empty_string(new_kind, new_length, new_ascii)\n",
" for i in range(len(a)):\n",
" _set_code_point(result, i, _get_code_point(a, i))\n",
" for j in range(len(b)):\n",
" _set_code_point(result, len(a) + j, _get_code_point(b, j))\n",
" return result\n",
" return concat_impl\n",
"\n",
" if isinstance(a, types.UnicodeType) and isinstance(b, types.UnicodeCharSeq):\n",
" def concat_impl(a, b):\n",
" return a + str(b)\n",
" return concat_impl\n",
"\n",
"True\n",
"@overload(operator.add)\n",
"@overload(operator.iadd)\n",
"def charseq_concat(a, b):\n",
" if not _same_kind(a, b):\n",
" return\n",
" if (isinstance(a, types.UnicodeCharSeq) and\n",
" isinstance(b, types.UnicodeType)):\n",
" def impl(a, b):\n",
" return str(a) + b\n",
" return impl\n",
" if (isinstance(b, types.UnicodeCharSeq) and\n",
" isinstance(a, types.UnicodeType)):\n",
" def impl(a, b):\n",
" return a + str(b)\n",
" return impl\n",
" if (isinstance(a, types.UnicodeCharSeq) and\n",
" isinstance(b, types.UnicodeCharSeq)):\n",
" def impl(a, b):\n",
" return str(a) + str(b)\n",
" return impl\n",
" if (isinstance(a, (types.CharSeq, types.Bytes)) and\n",
" isinstance(b, (types.CharSeq, types.Bytes))):\n",
" def impl(a, b):\n",
" return (a._to_str() + b._to_str())._to_bytes()\n",
" return impl\n",
"\n"
]
}
],
"source": [
"import operator\n",
"import inspect\n",
"\n",
"res = []\n",
"for r in f.typingctx._registries:\n",
" for func, F in r.globals:\n",
" try:\n",
" if func.__name__ == 'add':\n",
" overload_func = F.templates[0]._overload_func\n",
" source = inspect.getsource(overload_func)\n",
" print(func == operator.add) \n",
" res.append(func)\n",
" print(source)\n",
" except:\n",
" pass"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"func = res[0]"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<function numba.cpython.unicode.unicode_concat(a, b)>"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"func"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'templates': (numba.core.typing.templates.Registry.register_global.<locals>.decorate.<locals>.Template,),\n",
" 'typing_key': <function _operator.add(a, b, /)>,\n",
" '_impl_keys': {},\n",
" '_depth': 0,\n",
" 'name': 'Function(<built-in function add>)',\n",
" '_code': 881}"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f.__dict__"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"group = np.random.randint(0, 20, 10000)\n",
"values = np.random.randn(10000)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"@njit\n",
"def filter_values(arr, vmin, vmax):\n",
" #r = List()\n",
" s = 0.0\n",
" for i in range(len(arr)):\n",
" v = arr[i]\n",
" if vmin < v < vmax:\n",
" s += v\n",
" #r.append(v)\n",
" return s\n",
" #return r"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"27.5 碌s 卤 721 ns per loop (mean 卤 std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit filter_values(values, 0.4, 0.41)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"27.1 碌s 卤 89.6 ns per loop (mean 卤 std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit filter_values(values, 1.4, 1.41)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[(array(float64, 1d, C), float64, float64)]"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"filter_values.signatures"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"26.6 碌s 卤 101 ns per loop (mean 卤 std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"%timeit filter_values(values, 1.4, 1.41)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.04 ms 卤 96.9 碌s per loop (mean 卤 std. dev. of 7 runs, 1000 loops each)\n"
]
}
],
"source": [
"%timeit filter_values.py_func(values, 0.4, 0.41)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"15.5 碌s 卤 327 ns per loop (mean 卤 std. dev. of 7 runs, 100000 loops each)\n"
]
}
],
"source": [
"vmin, vmax = 0.4, 0.41\n",
"%timeit values[(vmin < values) & (values < vmax)].sum()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"@njit\n",
"def group_mean(group, arr):\n",
" d = {}\n",
" for i in range(len(group)):\n",
" g = group[i]\n",
" v = arr[i]\n",
" if g not in d:\n",
" d[g] = 1, v\n",
" else:\n",
" n, s = d[g]\n",
" d[g] = n + 1, s + v\n",
"\n",
" r = {}\n",
" for key, val in d.items():\n",
" r[key] = val[1] / val[0]\n",
" return r\n",
"\n",
"s = pd.Series(group_mean(group, values))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 0.032395\n",
"1 -0.019159\n",
"2 0.074304\n",
"3 0.038414\n",
"4 0.037896\n",
"5 -0.006833\n",
"6 0.028384\n",
"7 0.012377\n",
"8 0.019447\n",
"9 0.039177\n",
"10 0.063982\n",
"11 -0.007606\n",
"12 0.011451\n",
"13 0.049213\n",
"14 -0.014121\n",
"15 0.049763\n",
"16 -0.011851\n",
"17 0.024308\n",
"18 -0.012959\n",
"19 0.067572\n",
"dtype: float64"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s.sort_index()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"459 碌s 卤 48.5 碌s per loop (mean 卤 std. dev. of 7 runs, 1000 loops each)\n"
]
}
],
"source": [
"%timeit group_mean(group, values)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"7.55 ms 卤 620 碌s per loop (mean 卤 std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%timeit group_mean.py_func(group, values)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 0.032395\n",
"1 -0.019159\n",
"2 0.074304\n",
"3 0.038414\n",
"4 0.037896\n",
"5 -0.006833\n",
"6 0.028384\n",
"7 0.012377\n",
"8 0.019447\n",
"9 0.039177\n",
"10 0.063982\n",
"11 -0.007606\n",
"12 0.011451\n",
"13 0.049213\n",
"14 -0.014121\n",
"15 0.049763\n",
"16 -0.011851\n",
"17 0.024308\n",
"18 -0.012959\n",
"19 0.067572\n",
"dtype: float64"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vs = pd.Series(values)\n",
"vs.groupby(group).mean()"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numba as nb"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.52.0'"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nb.__version__"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"buf = ffi.cast(\"unsigned long long *\", id(mylist._opaque.data))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"buf1 = ffi.cast(\"unsigned long long *\", buf[0])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"buf2 = ffi.cast(\"unsigned long long *\", buf1[1])"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"140730398914576"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"buf2[1]"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [],
"source": [
"buf1 = ffi.cast(\"unsigned long long *\", buf[2])"
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"buf1[4]"
]
},
{
"cell_type": "code",
"execution_count": 98,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2246896846616"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mylist._opaque.data"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"140731033875760"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"id(type(mylist._opaque))"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {},
"outputs": [],
"source": [
"buf2 = ffi.cast(\"unsigned char *\", buf[1])"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 60,\n",
" 6,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 48,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 41,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 1,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 16,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 190,\n",
" 252,\n",
" 189,\n",
" 235,\n",
" 98,\n",
" 16,\n",
" 0,\n",
" 144,\n",
" 128,\n",
" 77,\n",
" 26,\n",
" 37,\n",
" 11,\n",
" 2,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 0,\n",
" 59,\n",
" 0,\n",
" 8,\n",
" 0]"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"list(buf2[0:100])"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1.11011454e-311])"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.frombuffer(mylist._opaque)"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(memoryview(mylist._opaque))"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"['__class__',\n",
" '__delattr__',\n",
" '__dir__',\n",
" '__doc__',\n",
" '__eq__',\n",
" '__format__',\n",
" '__ge__',\n",
" '__getattribute__',\n",
" '__gt__',\n",
" '__hash__',\n",
" '__init__',\n",
" '__init_subclass__',\n",
" '__le__',\n",
" '__lt__',\n",
" '__ne__',\n",
" '__new__',\n",
" '__reduce__',\n",
" '__reduce_ex__',\n",
" '__repr__',\n",
" '__setattr__',\n",
" '__sizeof__',\n",
" '__str__',\n",
" '__subclasshook__',\n",
" 'acquire',\n",
" 'data',\n",
" 'refcount',\n",
" 'release']"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dir(type(mylist._opaque))"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"21"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"0x15"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "object of type '_nrt_python._MemInfo' has no len()",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-43-d0d773010e39>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmylist\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_opaque\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;31mTypeError\u001b[0m: object of type '_nrt_python._MemInfo' has no len()"
]
}
],
"source": [
"len(mylist._opaque)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"25.2 碌s 卤 960 ns per loop (mean 卤 std. dev. of 7 runs, 10000 loops each)\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x214ed33d820>]"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deXQc13Xmv9c7GjsaAEFiJUCA+w6JABdRomSZlCXLsiVbsi07sWL5zCQeOxkndsYzGWcSTzJ2kuMl3uRYE49jW5YiKbZsyZQocSfAfScIYuECklgbxI7eqt78UV2NBthAN9DVXfUa93eOjohGo/pVV73v3freffcxzjkIgiAIcTHp3QCCIAgiPkjICYIgBIeEnCAIQnBIyAmCIASHhJwgCEJwLHp8aH5+Pq+oqNDjowmCIITl5MmTfZzzgqmv6yLkFRUVOHHihB4fTRAEISyMseuRXidrhSAIQnBIyAmCIASHhJwgCEJwSMgJgiAEh4ScIAhCcEjICYIgBIeEnCAIQnBIyIPIMsfLJzrgC8h6N4UgCGJWkJAHuXB7EH/x7+dwuLVP76YQBuD3Fzrx23O39W4GQcSELis7jcioVwIAeCkiJwC8eOgaJM7x6JpFejeFIKJCEXkQb0ARckmmHZMIYNQXgEy7ZxGCQEIeRI3EAzJF5AQw7pNAYzohCiTkQTx+isiJCcZ8EmS6FwhBICEPMhGRU+clyFohxIKEPIiXInIiDLJWCJEgIQ9CETmh4gvICMicrBVCGEjIg4Q8cokmO+c7Y74AAJC1Mk/wBiT0Dnv1bkZckJAHoYicUBnzKYM6Cfn84MVD1/Dodw/q3Yy4ICEPogo5eeTEhJDr3BAiKbhHvOgZ9oILPHCTkAdRrRWKyIlxisjnFRLn4BzwS+JebxLyIF4/ReSEwmjQI6d7YX6gTmqrq7tFJG4hZ4yVMsb2MsaaGGMXGWNf0KJhycYToIicUFAjcgrI5weBkJCLm+igRdGsAID/yjk/xRjLBHCSMfYO5/ySBsdOGhMRubgXk9AGmuycX6jXWWQhjzsi55x3cs5PBf89DKAJQHG8x002FJETKmStzC/U6yzyXgSaeuSMsQoA6wEcjfC75xljJxhjJ3p7e7X8WE0IReQCT3gQ2jBOWSvzCnXpyLz2yFUYYxkAXgXwRc750NTfc85f4JzXcs5rCwoKtPpYzfBSRE4EGQt55HQvzAdC1op/nkfkjDErFBH/Oef8NS2OmWw8lLVCBFFXdkok5PMCKQUmO7XIWmEAfgKgiXP+T/E3SR8oIidUQpOddC/MC8gjV9gC4FkAOxhjZ4L/PaLBcZOKh7JWiCBjlH44r5BSII887vRDzvkhAEyDtugK1VohVMhamV9IlH6YOlA9ckKF8sjnF7SyM4WgiJxQofTD+YUakc93j1x4ZJnDJ1EeOaGgLgiiyc75AWWtpAjhF5AicoKqH84vQkI+3/PIRSfcG6OsFYLqkc8vUiFrhYQcE6mHAEXkxETWCkD2ynxAJo88NZgckadux+Wco2fIo3czDI8akQNkr8wHyCNPEeZLRN7Q5kb937+HWwPjejfFsMgyx7hfgsWkLI1I4duBCKLmN5CQC858ici7hjyQZI6uQRLy6fAEJHAOZDiUtXIUkac+lEeeIqgRucXEUjoiVyOOofFAlHfOX1RbJd1GQh4LX3zpNP76jYt6NyMuAimQtaLFDkHCo47ETps5pbNW1MmcIY9f55YYFzX1MDMUkevZGuNzpXsEGXaxZSQUkUvi9n2KyDExEmfYLQjMsCDoxLV+jHjFjWbVAWvII+45JBp1MZAqTqlstWmBJHPhAwOJ6pGnBuo2b067ZdqOO+6T8LEXGvHKiY5kNk1TQhH5uNgdL5GErJWgkNPmEjMTkGUMCx4YkEeeIqgjcfoMQu7xS5BkLvRN6yVrJSqqtZJB1kpMpFRETlkrYqNG5Ok2M/zTeOTq6yKP2mpELvJglGhGg9ZZJlkrMRGQOUa8AaEXTtHGEimCGpE7bZZpi2ap3rnIPpqXrJWojPvJWpkNkszBOTDiEzc4oAVBKYIakWfYzdOmH6pC7hN4ZnvCWhG30yUa1SNXJzsFDjSTgtpfRH7Ko1orKUIoIp/BIw9ZKwJH5BPWCkXk0xGyVoIeOe0SNDNqfxH5KU+mrJXUwBuQYbOYYDObokbkIo/aofRDgTtdohmfkrUisvebDKQUishFftomIYeSkWK3mGA2sekjckmd7BT3YvvIWonKmF+C1cxgMytdg1Z2zkwqROQT9cjFDdJIyKGIs8NqDi7RnyZrJQWE3EvWSlTGvAE4bRaYqWhWTKj9Zdgr7j2lXmOR+zYJOZSROFpEHkiBCRE1Ivf4ZaHPI5GM+SQ4bWYwRccpIo/CREQu7lNeQJ4I0kTNUiIhx9SIPIq1IvCESLh4i+xpJpIxv4Q0mxmmoJKTRz4zE1krAkfkYV3aL+ievSTkCPfITeA8cuedmOwUV8jDJ3NIyCMz5g0gnayVmJCDOeSA2PMuEuewWxQpFPVJlYQcijjbLSZYzErnjRSVB1JgZafXr2TnAGJPTiWSMZ8akSs/k7UyPeH9ROSIXJI5nDYzAHEDNRJyKOLssJpDUVgkn1x95BJ5Ga9PklGQYQdA9VamY9yveuTT3wuEQvh3I2pErj59p1lJyIXH4w9G5CY1Ir/7YqaCteL1y8jPsAGIz1q51jeK7d/ci+4U3P9zzCcp1kpQyCkgn57wfiLqE5664CtNjcgFTUEkIUdsEXn4zLao+CQZ+WpEHkfHa+kZwXX3GNp6R7RqmmEY8wYUayXYM8hamR5pkrUiZkSunoMzuCOUqIuCSMgRKSKf3loRdcQGlLYXZMZvrajzBGNecb+L6RgLWitq1got0Z+ewCRrRdCIXJ4akZOQC8tERK58HREj8hRYEOSTZOQ4bTCx+PJ+1XmCMYEHtelQ8sgtISEXNa84Gaj9hDGBI3KuRuTkkQtPTBF58LWAzEOiLhKyzOGXOBxWEzId1riyDEJCLvC2d5EISDJ8AXlSRE5zndOj9pPsNKuwHrksTxVyMYMTEnIoF88e7pFHWBTgDxupRfTR1DbbLCZkpVniyjJQjzXqE/Omnw71CcMZ5pFT1sr0qP0k12mDNyALmdEVslasQY9cwHMASMjBOVdWdk7KI4+QtRL2mog+mtpmu8WMTHt8EVSqRuRq5cNwa4UmO6dH7RO5TisAMXPJyVpJEXySDM4xOSKfYbITEPNieyVFpNSIPB5PUz3/VIvI1Vrk4dYK6fj0qP0k16mktIqYS37XZOd8tlYYYy8yxnoYYxe0OF4yUUUpmkceCBNyER+/JiJyE7Ic1jizVoIRucDbe0VC3R0ozWaGmayVqKj9JDddXZsgYEQ+1SMX8Gkb0C4i/1cAOzU6VlLxBH1Re7SslXBrRcBRW/W17RZlslMLa2U0xdIPQ/t12iyhlZ1krUzPRESuWCsiVkBUu7Uq5CLOfwEaCTnn/ACAfi2OlWhkmU+KJNUR2BFjHjkgqLUSHpFPsVZmm2LnS9GIXLVW0shauYtI94jaT3KcAkfkoZWdymRnLBE55xwjBpsfsujdgGQx4g3gH99uxpvnO9Ez7MXyoiz8/UdWh0biyR55pCX6qRGR24LWyrA3gJ5hD/5x9xW8fLID71u+AF96/1LULMiM4VjK+aeaRz4YfEpx2swhK226ssbzAc45Gtv78f19rWhoc+MzWxfjTx+qCfnJaj/JS1c9cgGF/K5aKzPf02+cvY2v/64JXUMerC3Jxkc2luATm8pD2qEXSZvsZIw9zxg7wRg70dvbm6yPBaBcrC/88jT+X8N1rCvNwecfWILBcT8++S9HcfL6HQBTPPII6YfhHVpEH01dkWq3mEMbC+/4h/149dRNfGD1QjS0ubHzWwfw56+cDWVvTEcqZq1cd4/ib3/XhKIsB8pdTizMcSDNasY/v9cS9ftIRTr6x/CRHxzBMz9uRFPnMB5cXogXDrRj17cPoKHNDWCin0xkrYh3P6jWmcNqAmMzP22/evImvvDSaSzIduC/7FiCgMzxV7++iG/ubk5Wc6claRE55/wFAC8AQG1tbVLDnG/ubsa7l3vwN4+vxLP1FQCAp2pL8fQLjfjL184DQAzVD8MjcvGEPDwiX5STBgBYV5qDr31wJZYUZuDOqA/f39eKHx+8CleGHV/ZtWzaY01MdqaGwPUMefDsT47BL8n4xefq4bRZ4LRZ8O2n1+Fz/3YSf/byGXzv4xtg0jnqShayzPGlV86ipXsEf/OhVXhqYwkcVjOOtPbhK6+dxzM/bsS/PbcpVOo3O80GxsTMWlEHIzNjsFtM0/btl4934MuvncOWqnz8+FO1SLOZ8WcPL8VXXz+PH+5vw9qSbOxavTCZTZ9Eyqcfvn76Jn64vw2f2FQWEnEAKM1z4pefrUNRlgMAotcjn+SRiydg4R75zpVFeOdP78PPnrsXSwozACiZB1/9wAp8tLYEPz7YjqbOoWmPlUoe+ag3gE+9eAx9I1786x/ei+owa+nhlUX46iPL8daFLvyf3Zd1bGVyeflEB45e7cdXP7Acz9aVwxG0HTYvycfvv7gN5S4nvvbGRXgCakorQ4bNIuTqTjUiN5uUDbcjZaTtvdyDv3j1HLZVF+BfPl0bspYA4K8eW4H1ZTnBgW84ae2eilbph78E0ABgKWPsJmPsOS2OGy8Xbw/iy6+eR11lHr72wZV3/b7M5cQvn6/DxzeVYeWiLFhmyFpJpYjcZGKoXpAZyswI5789shw5aVb85Wvnp93mzJdCeeTHrvXjctcwvvHkGqwrzbnr989tXYxP1pXhR/vbsedStw4tTC49Qx58/c0m1FXm4WP3lN71e6fNgv/xgRVo7RnB/z18DQBgNpmQlWYV0lpR+7rZxGC3mu8K0mSZ4xu7m7E4Px0vPLsxNKip2C1m/OATG5Fms+BzPzup2ySoVlkrz3DOF3LOrZzzEs75T7Q4brx8e08LnDYzvv+JjbCaI59quSsd//uJ1ch0WEPWynS1VmzBYwjpkQcmPPKZyHHa8N8fXY4zHQP4+bEbEd+jDgqp4JF7goOR+mQyFcYYvvbYSpTlOfG9fa0pX0Trf/7mIrwBGX/34TURB3oAeHB5IbbXFOBgSx8AwGJiyHRYxJzsDF5PkylorUzp229f6kZT5xA+v2PJXSKuUpTtwHefWY/2vlH84uj1hLc5EilrrVztG8U7Td14tq48NKseDdVamS5rJd0enNkWMNdUjaLVrd5m4kPrirF1ST6+8dZl9A577/q9N6z6oeibE6v2gGOGAc5iNuGPti3G6RsDOBGcHE9F9jb34K0LXfjCg9VYnJ8+7fsYY/irx1aEkgPMJqZkQgko5Or9G8kj55zjO++2YHF+Oj64dtGMx6mvcqGuMg8/PXJdl6J6KSvkLx66CqvJhGfry2P+m2h55Ol2NddUPEtBTa1zxCDkjDH89eMrMewN4FfH747K1Zud8wkhFJXQ3IF15u/lqY2lyHVa8aP97cloli786+FrKMpy4Pn7KqO+t6ogA5/ZuhiAkrqX6bAIuSAo3FqxWcyThPztS9241DmEP3lgCSzTPNGH85kti3FrYBy7LybfgktJIR8Y8+GVkx14fN0iFGY6Yv67aCs7M1QhF9Aj/935Liwryoz56aSqIAObq1z41YmOu6Lu8Akh0Vd3TpRomNlySrOZ8Wx9BfY0dafkzkg374zhQEsvPlpbMq0NOZX/+nANfvqZe1GRn6545F7xInK1r5tCEblyP3PO8e09LahwOfH4upmjcZUHly9AucuJnxxK/mCfkkL+86M34PHL+KNt0SOLcGbMI5d4WGEdsYS8pXsYZzsG8OTGkml9z0h87J5SdPSP40gwb1jFFxaFi565opZocESJyAHg0/XlsFtM+JeDqReVv3y8AwDw0QgTnNNht5ixvaYAAMSNyIMeucU82Vp573IPLnUO4fM7qmOKxgElqv+DzRU4dWMAp28k14JLOSH3BiT865FruK+mAEuLoq9SDEdd5dk9fPemwn5JhtVsmjRqi8IrJ2/CYmJ4Yn3xrP7u/SuLkOO04qUp9opPkkMr4USMyBvb3fjUi8cQkORQx7XF0FldGXY8ubEEr566FXHuQFQCkoyXT9zEfdUFKMl1zukYWQ4rRrwB4SaDJ0Xk1glr5d9P3kR+hj3maFzlqdpSZNoteDGY0ZMsUk7ID1zpQ++wF5/ZUjHrv3Vl2LGmJBu/v9B11+8CMofVzGCLMLNtZPySjNdO3cKOZYVwBTdejhWH1Ywn1hfj7Yvd6B/1hV73BeSQRTPuFy8K+8mhqzhwpRfDngC8AQkWE4s56vqjbZXwBWS8crIjwa1MHvuv9KJryINn7o09Gp9KpsMCSebCLRKLlEc+7PHjvcs9eHTNwpjvC5UMuwVP31uKN8934vbAeCKaHJGUE/KuQeXLW7Eoa05//+iahTh3cxDX3aOTXg+EInKzUNbK/uZe9I148VTt3Drpx+4phU+S8dqpm6HXvAEZOcFl2aJF5MMeP/Y3KyUiAjKHxy9Pm1YWicX56VhXmoO3zt892IvKL491ID/DjgeXL5jzMTIdwQqIgmWuqAkmZsZgtypP2+9c6oY3IOOxtXNbqflsXQUkmeONs7c1bOnMpJyQ3xlTbqSctNgm9abySHCZ7W/PdU563S9xWEziWSuvnOxAfoYN9y8tmNPfLyvKwrrSHPzqeEfosdkXJuSieeR7mrpDefCSzJVt/mLI5AnnkdVFOH9rEB39Y4loYlLpHvJgb3MPntwY+yRnJLLSlEQA0RYFhawVE0J55G+cvY3inDRsKMud0zHLXE6sKs7CWxGe7BNFCgq5D5l2S0z50pEoyXVifVnOXUIekGVYzeqoLUZE7h7x4t2mHjyxvjiuTvr0PaVo6RnB6Y4BAKqQKwOlaBH578Kuq1+S4Q1uvD0bdq1SBvu3LnRGeafxef30LUgyx9OzmOSMRCgiF2yZ/qSVnRYT7oz5cLClD4+uXTirxICp7Fq1EGc6BtA5mBx7JeWEfGDMj5x0a1zHeHTNIjR1Dk1KM/NLHJagtSLKDkFvnL2NgMznbKuoPLJmISwmht0XuyDLHAGZhyreiRSRD477ceBKHwoylbkCSebwBGZnrQBKnZ7Vxdl4MwXslXcudWNVcRYqZlgAFAtZwYqawlkrataKicFuMWPMJyEgczy2ZnaTnFN5/8oiAMDuJEXlKSfk/aO+0B6Cc+UDqxeCsbujN6tp5gppRmNPUw+qCzNiqjE+E1kOKzZV5uHdpp6QLaF+xyLVW9lzSbFV1E4akDm8fmlOT2+7VhfhTMcAbiVxQktr3CNenLpxBw8um7s3rlLuSgdjwNmOQQ1aljzkKXnkAFBZkI6Vc5xjU1lSmIHqwgz8/iIJ+ZwYGPOFHvvnSlG2A/eU5+G35yYmKwISn8g1FWBl57DHj6NX3dixvFCT4z20fAFae0bQ3KVUeMtOs4Ixseqt/O58J4pz0rCxXPE+FY9chn2WETkAPKLaK+fFtVfeu9wDzoH3rYhfyPPSbVhfmoO9zT0atCx5TLVWAOCxNYvislVUdq4qwrGr/XCPJD5VNeWE/M6YH3nO+KwVAHh07UJc6R7BlWBpyoAsK9aKVYyslUMtffBLHDuWaifkAEITOHarGU6rWZiIfHDMj4MtvXhkdVFYuWIZHv/sJzsBoCI/HSsWJndCS2vebepBUZYj7uhTZceyQpy7OYieCOswjEqoaBZjcAZXbj8Wpa5KrOxcVQSZK/ZVoklBIY8/IgeUyQrGEMop90tcKGvl3cs9yHJYQtFnvJTmOVGzIANvBiNQu9kEp90ijEe+p6kbfonjkdULQyt41Yh8th65yiOri3Dy+h10DYojXCoev4QDLb14cHmhJtEnAOwIWjT7Lid3B7B4kMMi8ic3luDHn6qdthLmbFmxMAuleWlJsVdSQshlmcMvyfBLMoY9gbg9cgAoyLRjdXE2DlwJ5hxLwYhcgPRDWebY19yD+5cWznpBw0w8tHwBbgRT7mwWE9JtZmEWgOy/0ov8DBvWluRMKlfsDcw+a0VF3RFmd5J8UC1paHdjzCfhIQ1sFZXlCzOxMNuB9y6LY68EwoQ8P8Ouic2kwhjDzpVFONzaFypalyhSQsj/z+7LeOy7hzAQzCHPizNrReW+6gKc7hjA4LgfflnxyEVY2Xnu1iD6Rnx4UCN/XCV8wYjNYoLTZhEi/VCWOQ619mFbdQFMJhbaQCQgKZOdcxXyqoIMlLucocFeJN5t6obTZkZ9pUuzYzLGcP/SQhxs6RUmsyt8ZWci2LlqIfySElglkpQQ8rfOd+Fy13BoIk4LawUAti8tgCRzHGntU1Z2msRY2fleUzdMDKGCRlqxrjQH+RnKd2szm+C0mYWwVi7cHkT/qA/31eQDQFhELsdlrQDAtup8NLS7hREuQKns925TD7ZV58d17pF4cFkhRn0Sjl3t1/S4iSI02amRvTSVdaU5yE6zhjbhSBTCC/l192jocV8d9bSwVgBgfWkOMh0W7G3ugcwnKqT5DG6tvNfcg43luZoNaCpmE8MDwclTm0XxyEWY7FQj5m3VysBmNYd75HOPyNVjjvkknEpytbt4uHh7CJ2DntAEtpZsXuKCzWISxl6ZWNmZGCE3mxi2LHHhYEtvQguKCS/k4SPdvmCHzdEgawVQdobZUpWP94KTN1azyfArO7sGPbhwayg08aQ1O1cpCx1ynTbFIxcg/fDAlT6sXJSF/GDRsHCP3OOXo9Yin4nNVS6YTQwHW8SxV9673APGgAeWaWu9AcqenpurXMKkISbaWgGUwb57yIvWnsTVsRdeyA+19KE4Jw0rF2WFvqhYN0+Ihe1LC9AXzAO1mlnIWjFquU61A2ntj6vsWFaI335+K1YVZ8Fpsxh+snPY48epG3cm2UyhTbYlJSKPpRb5dGQ6rNhQlpPwR2ctOdLWhxULJwY2rdmxrBBX+0bRLsAGHOFFsxLF1iWKpZfIe0RoIQ9IMo609WFbdf6kHdC1slYA4L4pAqA+hvsMum/nodY+LMx2oFqjFKqpMMawqjgbjDGk280YNbhHfqTNjYDMJ11HNfrySTL8Eo8rIgeUiOv8rcFJpX6Niscv4dT1AWyu0m6Scyqq/bZfgElgdX9eUwKVsDTPicX56Ql9ahNayM/dGsSQJ4Ct1flYGxRyu8UU2slHC4pz0kJ5pVbzxOovI9ornHM0trlRX+nSLDd4Jpw2C8YMnrVy4Eov0m3mSZXs1AVBo0FbKNp+ndHYVp0PzoHDrcaPyk9dvwOfJKM+gUJemudEWZ7zrp2ljIgaj1kSqeRQovLG9v6EpS4LLeSHWvrAGLClKh/rg0KuZTSucl9wkkxd2QnAkCmIV7pH4B71oS6BnTScdJsZPkk2bMYG5xwHWnpRX5U/qZ6KGpGrtlAsG1LPxJqSHGQ5LEL45A3tbphNDPdU5CX0czZXudDY7o64/62RmFjZmdjP2Vadj/Hg01AiEFrID7b0YnVxNnLTbagsyECG3YJcDf1xle3BWt4WU3hEbrxItKFNiQgT+dgcjrqkedygPvk19xg6+sexPZh2qGINRl8joYg8vic4s4lha3U+Drb0GXbuRKWhzY1VxdmhsrOJYvOSfAx7Arh429hFtMa8AVjNLOFPsHXBSfFDrYkZ7IUV8hFvAKdvDIQmEswmhvtq8lGzQHtveHOVC3/8QBW21xQY2lppaHejNC9tzvsuzpb0oIU1ZtDt3tQI+b4p+fRmsxqRB4U8zogcUHzyzkHPpNLHRmPUG8CZjsT64yrqQiOj2ysnb9zB6uLshH9OlsOK9aWJmxQXVsgbg5NYW6snoq3vPrMB3/rYOs0/y2o24c/fvwyFWY4JITeYtSLLHI3t/Zqu1IuGOhdh1NWdje1uFOekodw1uda2WmtFbXe8k52A8ugMKKmORuXE9TsIyDwp90hBph01CzIMPW8w4g3g3M3BhM4XhLO1Oh/nbw3iTgImxYUV8kOtfXBYTZOKQplNiX9EUju90bJWLnUOYXDcn7SbEgDSbYq1YsTVnZwrA1tdBNEymyZPdsaTfqhSkutEucuJhnbjRqANbW5YzQy1FdoUUovG5qp8HL/Wb9g5lOPX+iHJHPWV+dHfrAHqpHgi7hFhhbyx3Y17KvI0iaZmw0REbqwotDF4cyTrpgQAp924EfmV7hH0j/pQV3n3pJ5lymSnVvdQ3WIXjl3tD1XUMxoN7W6sLcmBMzgAJ5r6Khc8fhlnOhIzwRcvjcGBTasKodFYW5KDN/5kK3YGdw/SEiGFvH/Uh8tdwxGjrUSjpqoZzSNvaHNjcX46irIdSftMI0fk6sA2U0Q+olH6oUp9lQuD4340dQ1pcjwtGfL4cf5mcvxxlbpKF0xMWYBkRBra3VhfmqtpuvJMWMwmrC7JTkg5ACGF/NhVtZMmNoUqEmr0ZiQhD0gyjl3tT6qtAgDpakRuwKwV1R8vzbt74lfNGVYHIIdGEfmm4P3Y2G68glHHr/ZD5khaaiqg7CK1qjgbR1qNYTf92ctn8P19rQCUge3CrcGkfh+JREghb2zvR5rVjNXFOdHfrDFGTD+8cHsIw95AUic6ASDDrqSwJWLyJh5kmePo1cj+ODCxIGhEnezUKCJfmJ2GCpcz9DRgJBra3LBZTJMWRiWD+ioXTnfcMcRT29H2fpy8phQ3O9auDGzJ7jOJQlAhd6O2IndOm+bGSygiN1DWSkPb9DZCIlmQZUdBpt1wlf9aehR/fLonFLWuRmhlp4b3UV2lMX3yhnY3NpTlaF62Nhqbq/LhlziOX9P/Hhn1BTDsUa55Q7sysK0vS34wmAiEE3I9/XHAmB55Q7sb1YUZKMhMTBGk6WCMoa7ShYY2t6EWwqgLozYtjmy9mUwMJhZmrWgobnWVxvPJB8Z8uNQ5hM1VyZsIV7mnIhcWEzPEU8qYV8JwcPBuaNNnYEsUwgm5nv44YDxrxReQceJa8v1xlfpKF3qGvbjaN6rL50eisb0fJbmR/XEVi8kUlkeuXTdQffIGAy2EaWzvB+fQ5R5x2ixYW5qj+/fhC8jwSTJGvH4MjPnQ1DWU1AyvRCOckDe0uZFmNWNNiT6PRD/7wH4AABy3SURBVDaDrew8d3MAYz4pqdkI4agDqlHypxV/3B31ic1sYhj3a5t+CIT75MaZ8GxsV/rMWp36TH2lC+dvDYayhPRALSMx4gnoOrAlCk2EnDG2kzHWzBhrZYx9RYtjTkdjez9qK3Jh1XBT4dlgC36uURY5NLS5wRiwabE+N+Xi/HQsyLIbRriu9Azjzpg/qpBbwlLAtIzIAdUnN07BqIY2/eaUAEUwJZnjuI7bv6nllke8AZzuuAOb2YS1pYlfmp8s4r6yjDEzgO8B2AVgBYBnGGMr4j1uJNwjXjR36+ePA0ouaLrNjN5hr25tCKeh3Y1lRVkJKRYWC6pP3thuDJ+8sS02602tt2IzmzTP662vcmHIE0BTp/4+eV+wz+gZfW4oy4XNbNLVJ1fnQ/wSR9+wD5kOS9IXEyYSLYboewG0cs7bOec+AC8BeFyD496FuqGrnkIOABvKc3H8mv4RqDcg4eT1O7qnUNVVutA77EVbr/4+eayFw9Rccq2jcWDi6cgIE3wTK371u0fSbGasK83R1X4LX33sHvWmzCSnihZ3cTGAjrCfbwZf0xzV61tTou8jUV2lC5e7hnXfEeb0jQF4A4ndJCAWVJHQW7hC+eMx2EyqtaJVDnk4RdkOLM5PN4Td1NDmRobdkpQKfzNRV+XChVuDGPL4dfn88J2s+ka8CbnueqLF2UR6Lr3rGZsx9jxj7ARj7ERv79xq8j5VW4r//eFVuvnjKupj+zEdPT9AKRFqYsC906TZJYtylxNFWQ7dJzybu4cxEIM/Dkws00/U43VdZZ4hfPKGdreSAqhzn6mvdEHm0M0nD9/JqnfYq9lqXqOgxdW9CaA07OcSALenvolz/gLnvJZzXltQUDD11zGxqjgbT6wvmVsrNWR1cQ4cVn09P0Dxg1cVZyM7LbGbBESDMYb6KheO6uyTq9djUwypqerqzkRFZnWV+vvk3UMetPeO6pI/PpX1ZTmwWUy6pSGGR+TuEZ8mFS+NhBZncxxANWNsMWPMBuBpAL/R4LiGxWYxobY8T1chH/dJON2hvz+uUleZh74Rn64bKzTOYmONREfkRvDJVdHU23oDlEVXG8r088nDPfKAzMkjnwrnPADgTwDsBtAE4GXO+cV4j2t06irzcLlrWLc6Iyev34Ff4oYp+qPaGXpFXLPxx4EJjzxRkdmET66vkGc5LFi+MEu3NoRTX5mv1M0fS75PPrXWCwl5BDjnb3LOazjnVZzzr2txTKOzKShcx3TKXjnS1peUTXRjpSzPiUXZDt0m+C53Kf54rNGnOYFZKyp1lS4cvdqvm0/e0O7GpkpX6OlDb+oq88A5cPRq8ge3qTXzyVohAABrSrJ19cmVTQKykWFPziYB0dA7n3zCH49NyK3mxForgCJcwzr55DfvjOFG/5huK34jsa4sB3aLSRd7ZcwXQPh4RpOdBABFADaW5+oSgSZ7r8FYqatywT3qQ0tP8n3yxnY3yvKcKM5Ji+n95gRbK8CE3aTHYG8kf1zFbjGjtkKfPjPqCyA7zToxgJO1QqjULXbhctcQBsaS65Mne6/BWNErn3yi/njsNpMlwZOdALAgy4FKnXzyhnY38tJtqCnMTPpnz0R9pQtNnUNJn1sa80pIt1tCT7BkrRAhNlW6wHny88mTvddgrJTkpqE4Jy3pE56Xu4YxOB5b/rjKRNZKYrvAJh18cs45GtvcqKvMS8i2YvGgXqNk++SjvgDSbRZkOpRUXZrsJEKsLc3WxfM70pbcvQZjRfXJjyZ5Y4WGWfrjwMQS/UR3aNUnv3Q7eT75dfcYbg96UG+A/PGprCnJQZrVnHR75Ub/OPIzbRMROXnkhIrdYsa9i/OSuifh4LgfF28bzx9XqavMQ/+oD1d6hpP2mUda+1Duit0fB5IXkYfSMtuTtwFxgwHqq0yHzWJCbUVuUp/a3CNeNHUOob7ShQwHWStEBLYsyUdz9zB6hjxJ+bxjwU10jSvkQZ88SR3VL8lobHdj65LZRZ+JrLUSzoIsByoL0nEkicLV0OZGQaYdVQXpSfvM2VBX6UJz9zDcI8mpIKoObFuW5CMz5JFTRE6EoQrI4bbkRFwNbW7YDbzXYGmeEyW5aUl7dD7TMYBRn4Rt1bMU8iSkH6psW5KPo+39SdlVinOOhnY36itdYMxY/riKGoQcjXNuKSDFtifA4dY+ZAYLh6kReaKfxJJNap2NDqxYmIUcpxWHWpITcR1p68PG8lxD11Kuq3Sh8ao7KT75oZY+MIZZZ/BMeOSJ7wJbqwsw7pdw6vpAwj+rrXcEvcNeQ+WPT2V1cTbSbea47JWrfaNY/le/R3NXdAvvUGsf6qpcsJhNYVkrxu0/c4GEPE5MJoYtVfk43NqX8IUw6sbTRvQ+w6mvdGFgzI/m7sT75Idb+7CmOBvZztkVDkt0rZVw6irzYDYxHGqdW9XP2WDE/PGpWM0m1FbkxZUk0DkwDr/Eoy62uuEeQ0f/eOjJmTxyYlq2LMlH15An4RsrHG03ficFEKr/kmhfeNjjx+mOAWyZpT8OhOeRJ74LZDqsWF+ag0MtibffDrT0oTgnDWUzbDxtBOqrXGjtGZnzTlv+4NNed5S5KdXyVO8R1SOnBUHEXYR88tbEdtSGdjecNv02no6V4pw0VOan48CVxEagx4L52Vtn6Y8D4Ss7k9Oht1bn49ytwYQuhPEFZBxp7cP9SwsM64+rxLt4TPXHu6II+aHWPizImpj4DeWRG9ianAsk5BpQ5nKiNC8NhxIo5Jxz7GvuxabFebptojsbti8tQEO7O7R7eSI42NIHh9WEDWWzXxg1MdmZnO9yW3UBOE/sU8rJ63cw6pOwvWZu9f6TycpFWchyWHCwZW6DvV9SIvKeoekjelnmaGhzY8uS/NDARis7iRnZuiQfjW3umGfSZ0t73yhu9I9hx7LChBxfa+5fWghfQE7o8vTDrX24pyJvTlF1aM/OJHXotSXZyHRYEuqT77/SC6uZYfMcrKZkYzGbsK2mAPuae+c0t6SulJ0pIm/qGkL/qG9SampVYQZsFhMWzWLNgQiQkGvE1iUFGPYGcO7WYEKOv/dyDwBFIEVg0+I8OKwm7GvuScjxu4c8aOkZmXX+uErIWknSI7bFbEJ9pQsHWxI3Kb6vuQe15XmGqYgZjR1LC9Ez7MXFOax6DchKwDSTR65aneFzKOtKc3D5f+3EgizHrD/TyJCQa0R9lQuMIWETWvuae1FdmIFSg09iqTisZmyuysfeOUZc0VA76Vz8cSB5C4LC2Vadj5t3xnHdPab5sbuHPLjcNYztS41vq6hsX1oAxiaClNkQbq1Md3/tv6L0mamibbT6M1pAQq4Reek2rCnOxntzuCmjMeoN4OhVNx4QxFZRuX9pAW70j+Fqn/bZPPuae+FKt2F50dx2vzEncUGQytZqRWQPJmAuZX9wYlkEf1wlP8OONSU5eG8OT22qhemTZNyJsOPQ4LgfR9v78eDyBXG3UwRIyDXkfSsW4EzHQNSUqNlyuLUPfonjfoGiLQC4v0YZePY1a+sL+wIy9l7uwUPLF8w5ukr0Vm+RqHApq173JWCw33+lFwuy7FhWZKyytdF4YGkBznQMoH+W2Tz+sMVmXYN397d9zT0IyBzvW0FCTsySh1cWAQD2NHVrety9zT3IsFtQW26Mbd1ipczlRGVBOvZpnIbY0O7GsDeAh1fOvZNObPWWvIicMYaHVxThYGsfRryB6H8QIwFJxsErvdheY/y0w6nsWFYIzoH9V2Y3uElhSQXdw3cL+TuXupGfYcO6UmOn6moFCbmGVBdmoMLlxNsXtRNyzjn2Xu7F1iX5QqQdTuX+mkI0apyG+PbFLjht5jktBFJJ5oKgcHauKoIvIGs6CXz25gCGPAFsrxHLegOAVYuykZ9hx97LsxvsA2ERefeUiNwXkLG/uRcPLltgmP1KE414ymBgGGN434oFONLWh2GPNjuFX+4aRteQBw8sE8tWUbl/aQF8AVmzMq6yzPHOpW5srymIazFPMotmhbOxPBf5GTbs1nCw33u5FyaGOWfw6InJxHD/0gLsv9I7q9RddbITuDsF8ehV5YltvtgqAAm55jy8sgh+iWvmC+9tFivtcCr3Ls5DmtWs2STw2ZsD6Bn2xmWrAMDakhxsrpqoT50szCaG960owntN3fD4439K4Zzjd+c7UVfpmnW9GaPwwNJCDI77caYj9qJiquhnp1nRPWVR0J5L3XBYTXE9sYkGCbnGbCjLhSvdhncuaRNx7b7QhZWLsoTNe3VYzdixrBBvne+CX4PFUm9f6obFxLBjaXxCvmVJPn7x2TpdHr13rirCqE/CEQ1KH1+4NYSrfaP44NpFGrRMH7bV5MNiYrPqM+pkZ3FO2qTkAs6VJ7Zt1QWG20ErkZCQa4zZxPDg8kLsvdwDXyA+4WrrHcHZm4P40LpijVqnD4+vWwT3qE+THPu3L3YJHX0CSp2RTIcFv7/QFfexfnP2Fqxmhl2rFmrQMn3IclhxX00BfnP2dsx7mwYkGRYTQ1G2Y5KQX7w9hNuDHrxvnqQdqpCQJ4CHVxRh2BuIe3n666duwcQUIRSZ+5cWIsdpxX+cuRXXcVp7RtDWOxq3raI3NosJDy4rxDuXuuMq6SDLHL8914ntNQVCD2wA8OENxegc9MRco1ySOSxmhgVZ9klCvqepG4wBO5aLaUXOFRLyBLC1Oh9pVjN+f3HuEZcsc7x++ha2LMlHoaC2iorNYsIHVi/E7otdcaXd7Q5+nw+lQLS1c1UR7oz5ceza3HfJOX6tH52DHjwmsK2i8tDyBch0WPDaqZsxvd8vcVhMJizIcqBvxAe/JINzjt+cuY17yvOQn2FPcIuNBQl5AnBYzdi1qgi/OXMbo3MUrmPX+nFrYBwf2VCicev04Yn1xfD4Zbw9x8FNljleOdGBjeW5KVHw6L6aAjisJrx1fu6D/a/P3kaa1ZwS2RkOqxmPrlmEty50xdRnArIcjMiVIKdn2IuGNjfa+0bx9L2liW6u4SAhTxCfrC/HiDcwZzvh9VO34LSZhbcRVDaW56IkNw2vn57b93GwtQ/X3GN4tq5c45bpg9Nmwc6VRXj99K05par6JRlvne/EQysWwGkTo0hWND6yoRjjfimmuQM1Ii8KCnn3kAc/P3YD2WlWPLJa3PmCuUJCniDWl+Zg5aIs/Kzh+qyLRnn8Et4834ldqxamTCdljOFD64pxuLUPPRFW4kXjZw3X4Uq3YdfqogS0Th/+YMtijHgD+PeTsdkJ4Rxq6cOdMb/Q2SpT2Viei3KXE6/GYK8EJBlWM0NhlmKhXLw1iN0XuvDkxpKU248zFkjIEwRjDM/WleNy1zBOXr8zq79951I3hr0BfHiD2NkqU/nQ+kWQOfDG2c5Z/d3NO2N473I3nr631NCbTs+WdaU52FCWg58euTbrjapfOq5En/fVpE6uNGMMH15fgoZ2N24NjM/4XnWyU43If7i/HQGZ45l7y5LRVMNBQp5APrhuETIdFvys8fqs/u6l4zdQlOVAncE3WZ4tSwozsbZUEa7ZpGb+/OgNAMDHN6WGrRLOH25ZjGvusdDCr1i4eHsQuy9249ObK1JqYAOUuRTOgVejPKX4ZQ6ryYS8dBusZoZbA+Ooq8zDksKMJLXUWJCQJxCnzYKPbCjBm+c7Y95k9nBrHw63uvHc1sUpWSfiiw9V40b/GF46fiOm93v8En51vAMPLV+A4hSY5JzKzlVFWJjtwIuHr8b8N9/e04JMhwXPbV2cwJbpQ5nLiQeWFuDHB9tnrIgYkJTJTsYYCjOVqDwVB/pYISFPMJ+sK4df4njpWHThkmWOv3urCcU5aXi2PjVvyvtrCnDv4jx8593WmLIT3jzfif5RHz5VX5H4xumA1WzCs/XlONzqRnPXcNT3X7g1iLcvdeOPtlYiO03s3PHp+MtHlmPUG8B33m2Z9j1+iYcqWC7MdsCVbsP7UyQxYC6QkCeYJYUZeGBpAX6wvw3X3TNvsPDGudu4cGsIX3p/TcpO2DDG8OWdS9E34sX/jSEK3X+lF4WZdmxZklo2UzjP3FMGh9WE7+9rjfreb+25giyHBX+4tSLxDdOJmgWZePreMvxb43W0945EfE9AViY7AeAru5bhux9fn3I202wgIU8CX39iNcwmhj/91ZlpV/J5AxK+ubsZKxZm4fG1qTXJOZWN5Xl4aPkC/Gh/O+5E2VBgzCchL90mXJ3t2ZCbbsNnt1Xi12duz5jBcu7mAPY09eCz2yqR5UjNaFzlTx+qgd1iwt+9dTni7yWZh0oR11bkYXNV6kz6zoW4hJwx9hRj7CJjTGaM1WrVqFRjUU4a/vZDq3DqxgB+sK8t4nte2N+Om3fG8ZVdy1JyT8Gp/MXOpRjxBfB3bzXNmLHh8Usp+3QSzhcerEZ9pQtfff08LkXYjLhnyIMvvXIWOU4r/mBLRfIbmGQKMu34zw8swTuXuiMWF/NLMixmikNV4v0mLgD4MIADGrQlpXl8XTE+uHYRvv1uC86Glev0BWT89/84j3985wp2rSrCfQLtuRgPNQsy8fy2Srx84ib++BenMOab8MvD8+49fglp80DILWYTvvPMeuQ4rfhPPz+JwfGJRUI33GN48ocNuHlnHN//+AZkpng0rvLc1sUoyU3DZ396Ar8OW1jHOYfHP2GtEEBcq004500AUvqxV0v+5vFVOHGtH0/9sAEfWr8IH1izCN95twUnr9/B57ZX4s8fXqp3E5PKV3YtQ0GmHV9/swkdPxrDl3cuw6GWPrxy8iYeW7MQf/34Koz7JRRkpMaiqGgUZNrxvY9vwNMvNOJjP2rAJ+vKUZKbhr/493PwSTJ+8dm6ebN1GaAs23/5c/X4wkun8YWXzuBwax+21xTiXw6140zHAD68PrUtyNnAZrvqMOJBGNsH4Euc8xOxvL+2tpafOBHTW1OOm3fG8KP97XjlZAc8fhlOmxnfeHINHl2TOiv0Zsu7Td34L788jVGfBLOJwWJiuHdxHn723CY89E/7UbMgA9//xEa9m5k0fneuE99+9wqudCsTfYWZdvzsuU1YKtjGyloRkGR8a08LvrevFZwDJblp+E/3V+HJjSXzboKTMXaSc36XjR1VyBljewBEWhf9Vc75r4Pv2YcoQs4Yex7A8wBQVla28fr12S2SSTXujPrQ2O7G6pJslOQ69W6O7nQNenCmYwDry3Lw+V+chskEvPR8Pbb8/XvYVJmHf/roOr2bmFQ457h4ewi9I15sLM9N+cnNWLh4exDdQx5sqy6AdZ7649MJedRnVs75Q1o0gHP+AoAXACUi1+KYIpObbsOueVjcZzqKsh3Yma3ECxYzC6389Abmx2TnVBhjWFWcrXczDMXKRdlYuYi+k0jMz2GNMDRWsym0Ldy4b35MdhJEPMSbfvgEY+wmgHoAv2OM7damWcR8xmpm8EtcyU4IyHBYKd4giJmIN2vldQCva9QWggAwEZH7JQ5J5hSRE0QUKNQhDIfFbEJA5vAEJACYlx45QcwGEnLCcCjWigyPj4ScIGKBhJwwHFaTYq14/MqEJwk5QcwMCTlhOCxmhoDEMe5XInLyyAliZkjICcNhNZvgk2R4/Kq1QrcpQcwE9RDCcFgpIieIWUFCThgONf1QjcjtJOQEMSMk5IThCKUfUkROEDFBQk4YDluwzvSwR6lRTh45QcwM9RDCcKg7v6hCnmajiJwgZoKEnDAc6l6MoYh8ntWcJojZQkJOGA6bRY3Ile3OKCIniJkhIScMh8U02VqxW+g2JYiZoB5CGA51U91hrx8Oq4n2hCWIKJCQE4bDGjbZSXVWCCI6JOSE4VCFfMgToBxygogBEnLCcFiC1sqIx08ROUHEAAk5YTisYQuCSMgJIjok5IThmOyR0y1KENGgXkIYDjX9cNwvkUdOEDFAQk4YDptlIt2QhJwgokNCThgONSIHaJs3gogFEnLCcKgeOUBCThCxQEJOGA41awWgErYEEQvUSwjDYQmLyMkjJ4jokJAThmNyRE5CThDRICEnDEe4R04lbAkiOiTkhOEIF3IqYUsQ0aFeQhgOS5i1QhE5QUSHhJwwHNbwPHLa5o0gokJCThgOK0XkBDErSMgJw2E2UR45QcwG6iWE4WCMwRac8KT0Q4KIDgk5YUjUCU8ScoKIDgk5YUjUFERa2UkQ0YlLyBlj32SMXWaMnWOMvc4Yy9GqYcT8xkoROUHETLwR+TsAVnHO1wC4AuAv428SQUyUsqWInCCiE5eQc87f5pwHgj82AiiJv0kEAVgtakRO7h9BREPLXvIZAG9peDxiHqMuCiJrhSCiY4n2BsbYHgBFEX71Vc75r4Pv+SqAAICfz3Cc5wE8DwBlZWVzaiwxf7CaTWCMaq0QRCxEFXLO+UMz/Z4x9mkAjwJ4kHPOZzjOCwBeAIDa2tpp30cQgJJ+6LCYwRiL/maCmOdEFfKZYIztBPBlANs552PaNIkglIic/HGCiI14e8o/A8gE8A5j7Axj7IcatIkgYDUzylghiBiJKyLnnC/RqiEEEY7FZIKDCmYRREzQsythSKwWE5WwJYgYiSsiJ4hE8YdbKjDqDUR/I0EQJOSEMXlgaaHeTSAIYSBrhSAIQnBIyAmCIASHhJwgCEJwSMgJgiAEh4ScIAhCcEjICYIgBIeEnCAIQnBIyAmCIASHzVB5NnEfylgvgOtz/PN8AH0aNsdIpOq50XmJR6qem+jnVc45L5j6oi5CHg+MsROc81q925EIUvXc6LzEI1XPLVXPi6wVgiAIwSEhJwiCEBwRhfwFvRuQQFL13Oi8xCNVzy0lz0s4j5wgCIKYjIgROUEQBBEGCTlBEITgCCXkjLGdjLFmxlgrY+wrerdHKxhj1xhj54MbWJ/Quz3xwBh7kTHWwxi7EPZaHmPsHcZYS/D/uXq2cS5Mc15fY4zdCl63M4yxR/Rs41xgjJUyxvYyxpoYYxcZY18Ivi70NZvhvIS/ZpEQxiNnjJkBXAHwPgA3ARwH8Azn/JKuDdMAxtg1ALWcc5EXKgAAGGP3ARgB8P8456uCr30DQD/n/O+DA3Au5/zLerZztkxzXl8DMMI5/wc92xYPjLGFABZyzk8xxjIBnATwIQB/AIGv2Qzn9VEIfs0iIVJEfi+AVs55O+fcB+AlAI/r3CZiCpzzAwD6p7z8OICfBv/9UygdSiimOS/h4Zx3cs5PBf89DKAJQDEEv2YznFdKIpKQFwPoCPv5JlLnwnAAbzPGTjLGnte7MQlgAee8E1A6GIBU2pDzTxhj54LWi1D2w1QYYxUA1gM4ihS6ZlPOC0iha6YikpCzCK+J4QtFZwvnfAOAXQD+OPgYTxifHwCoArAOQCeAf9S3OXOHMZYB4FUAX+ScD+ndHq2IcF4pc83CEUnIbwIoDfu5BMBtndqiKZzz28H/9wB4HYqNlEp0Bz1L1bvs0bk9msA57+acS5xzGcCPIeh1Y4xZoYjdzznnrwVfFv6aRTqvVLlmUxFJyI8DqGaMLWaM2QA8DeA3Orcpbhhj6cHJGDDG0gE8DODCzH8lHL8B8Ongvz8N4Nc6tkUzVKEL8gQEvG6MMQbgJwCaOOf/FPYroa/ZdOeVCtcsEsJkrQBAMFXoWwDMAF7knH9d5ybFDWOsEkoUDgAWAL8Q+bwYY78EcD+UcqHdAP4ngP8A8DKAMgA3ADzFORdq4nCa87ofyiM6B3ANwOdUX1kUGGNbARwEcB6AHHz5v0Hxk4W9ZjOc1zMQ/JpFQighJwiCIO5GJGuFIAiCiAAJOUEQhOCQkBMEQQgOCTlBEITgkJATBEEIDgk5QRCE4JCQEwRBCM7/B4ggytGOcVCXAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"@njit\n",
"def get_peaks_y(x, y, xstart, xend, count):\n",
" index0 = np.searchsorted(x, xstart)\n",
" index1 = np.searchsorted(x, xend) \n",
" index0 = max(index0 - 1, 0)\n",
" \n",
" if index1 - index0 < count * 2:\n",
" return x[index0:index1], y[index0:index1]\n",
" \n",
" xstart = x[index0]\n",
" xend = x[index1]\n",
" xstep = (xend - xstart) / count\n",
" res_index = np.zeros(count*2, np.int32)\n",
" \n",
" pos = 0\n",
" xnext = xstart\n",
" i = index0\n",
" \n",
" while pos < count * 2:\n",
" xnext += xstep\n",
" min_val = max_val = y[i]\n",
" min_idx = max_idx = i\n",
" \n",
" while x[i] < xnext:\n",
" val = y[i]\n",
" if min_val > val:\n",
" min_val = val\n",
" min_idx = i\n",
" if max_val < val:\n",
" max_val = val\n",
" max_idx = i\n",
" i += 1\n",
" \n",
" if min_idx > max_idx:\n",
" min_idx, max_idx = max_idx, min_idx\n",
" res_index[pos] = min_idx\n",
" pos += 1\n",
" res_index[pos] = max_idx\n",
" pos += 1\n",
"\n",
" return x[res_index], y[res_index]\n",
"\n",
"n = 10000\n",
"x = np.linspace(0, 100, n)\n",
"y = np.sin(x)\n",
"y[np.random.randint(0, n, 50)] += np.random.randn(50)\n",
"xp, yp = get_peaks_y(x, y, -10, 28, 100)\n",
"%timeit xp, yp = get_peaks_y(x, y, -10, 28, 500)\n",
"pl.plot(xp, yp)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"4.58 ms 卤 271 碌s per loop (mean 卤 std. dev. of 7 runs, 100 loops each)\n"
]
}
],
"source": [
"%timeit xp, yp = get_peaks_y.py_func(x, y, -10, 28, 500)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.5"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment