1027 lines
217 KiB
Plaintext
1027 lines
217 KiB
Plaintext
|
{
|
|||
|
|
"cells": [
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 91,
|
|||
|
|
"id": "54634e7b",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"import matplotlib.pyplot as plt\n",
|
|||
|
|
"import numpy as np\n",
|
|||
|
|
"from numba import jit\n",
|
|||
|
|
"import itertools\n",
|
|||
|
|
"\n",
|
|||
|
|
"fontsize=14"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "8e090825",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__В pdf-файле задачи 1-4__"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 374,
|
|||
|
|
"id": "705fbec9",
|
|||
|
|
"metadata": {
|
|||
|
|
"scrolled": true
|
|||
|
|
},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"text/html": [
|
|||
|
|
"\n",
|
|||
|
|
" <iframe\n",
|
|||
|
|
" width=\"1000\"\n",
|
|||
|
|
" height=\"800\"\n",
|
|||
|
|
" src=\"HW1_14.pdf\"\n",
|
|||
|
|
" frameborder=\"0\"\n",
|
|||
|
|
" allowfullscreen\n",
|
|||
|
|
" \n",
|
|||
|
|
" ></iframe>\n",
|
|||
|
|
" "
|
|||
|
|
],
|
|||
|
|
"text/plain": [
|
|||
|
|
"<IPython.lib.display.IFrame at 0x7ff362a4af20>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"execution_count": 374,
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "execute_result"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"from IPython.display import IFrame\n",
|
|||
|
|
"IFrame(\"HW1_14.pdf\", width=1000, height=800)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "d4960afd",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"### Задача 5 \n",
|
|||
|
|
"\n",
|
|||
|
|
"- сгенерируем s точек равномерно в квадрате стороной 2, оставим только те, для которых $|| x ||_p \\le 1, p=1,2,3$\n",
|
|||
|
|
"- сравним реализацию с jit и без"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 2,
|
|||
|
|
"id": "55fc75ad",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"@jit(nopython=True)\n",
|
|||
|
|
"def numba_unit_disk(p=1, s=100000):\n",
|
|||
|
|
" X_old, Y_old = np.random.uniform(-1.0, 1.0, s), np.random.uniform(-1.0, 1.0, s)\n",
|
|||
|
|
" X_new, Y_new = [], []\n",
|
|||
|
|
" \n",
|
|||
|
|
" for x, y in zip(X_old, Y_old):\n",
|
|||
|
|
" if (abs(x)**p + abs(y)**p)**(1/p) <= 1:\n",
|
|||
|
|
" X_new.append(x)\n",
|
|||
|
|
" Y_new.append(y)\n",
|
|||
|
|
" return X_new, Y_new\n",
|
|||
|
|
"\n",
|
|||
|
|
"def unit_disk(p=1, s=100000):\n",
|
|||
|
|
" X_old, Y_old = np.random.uniform(-1.0, 1.0, s), np.random.uniform(-1.0, 1.0, s)\n",
|
|||
|
|
" \n",
|
|||
|
|
" X_new, Y_new = [], []\n",
|
|||
|
|
" for x, y in zip(X_old, Y_old):\n",
|
|||
|
|
" if np.linalg.norm([x, y], ord=p) <= 1:\n",
|
|||
|
|
" X_new.append(x)\n",
|
|||
|
|
" Y_new.append(y)\n",
|
|||
|
|
" return X_new, Y_new"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 3,
|
|||
|
|
"id": "108af16e",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAF2CAYAAABj+Z+GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoR0lEQVR4nO3db5Bc1X3m8e8jCUkOcoKQVAQElnBCYkuVLZydEGddlWyssY3zApENSeS0QGBcsmbihF1XssbFC285ca2drVpMaj0CLcYSqMuYkHVZqcTFMmPYVLmMw1CLDSOCkTGUEdgIAV4kGf0Z/fbFvW2a0fzpmb7d99/zqZqa7tu3u391Z6afOfece44iAjMzq69FeRdgZmb5chCYmdWcg8DMrOYcBGZmNecgMDOrOQeBmVnNOQjMzGrOQWCWI0k3SfqmpKOSfFGP5cJBYJavZcD/Aj6fcx1WY/KVxWYg6UHgX4HjwDXp5tuBT0TE6T68/1XA30WEev1eZlO5RWD2hgbJ38RvAR8FtgP/caadJb1N0pE5vm7tT+lmC7ck7wLMCuQF4M8jaSb/q6RfAT4O/PcZ9n8euHSO1/x/2ZVn1hsOArM3PBRvPlf6LeCvJP18RJzxgR4Rp4ADfavOrEd8ashsgXxqyKrCLQKzN/ymJLW1Ct4NPD9dayDlU0NWCQ4CszdcAHxe0gjwa8BfAn89085ZnBqS9DbgXGB9ev/S9KEDEXGkm9c265SHj5rxpuGjp4CtQAB3AP85IiZ7+L67gW3TPPS7EfFgr97XrJ2DwIyfBcHjEfGxvGsx6zd3FpuZ1ZyDwMys5nxqyMys5twiMDOrOQeBmVnNlfI6gtWrV8f69evzLsPMrFQeeeSRlyJizdTtpQyC9evXMz4+nncZZmalIunZ6bb71JCZWc05CMzMas5BYGZWcw4CM7OacxCYmdWcg8DMrOYcBGZmNZdJEEi6Q9KLkh6f4XFJ+ltJByR9V9Kvtz22TdJT6dd087KbmVkPZdUi2A1cPsvjHwQuSb+2AzsBJJ0LfAr4TeAy4FOSVmZUk5mZdSCTIIiIfwZenmWXzcCdkXgIOEfS+cAHgPsj4uWIeAW4n9kDxaywmk1Yvx4WLUq+N5t5V2TWmX71EawFfth2/7l020zbzUqh2YQVK0CCrVvh2WchIvm+dSsMDuZdodncStNZLGm7pHFJ44cOHcq7HDOGh5MP+6NHZ95nbAw2buxfTWYL0a8gOAhc1Hb/wnTbTNvPEBG7ImIgIgbWrDlj8jyzvhoehp07O9t3//5kf7Oi6lcQ7AOuSUcPvRv4SUS8ANwHvF/SyrST+P3pNrPCajY7D4GW+e5v1k9ZDR/9MvAt4FclPSfpekk7JO1Id/kn4GngAPA/gWGAiHgZ+Cvg4fTr0+k2s8LasWPufaazerU7kK2YMlmPICI+NMfjAfzpDI/dAdyRRR1mvdJswg03wOHDC3+Nw4dh+/bkdqORTV1mWSjlwjRm/dRswrXXwqlT3b/WsWNwzTXJbYeBFUVpRg2Z5eWGG7IJgZbTpz201IrFQWA2i8HB7k4HzWZszKOJrBgcBGbTaDZhyZLkw7qXdu50y8Dy5yAwm6LZTM7jT0725/3GxhwGli8HgdkUH/1och6/n3rd8jCbjYPALNVsJmP9Z5syopeWLfN1BpYPB4EZyQfwtm296xjuxIkTyWgiz01k/eYgsNprNpMP4H71Ccxl/35Y6zl4rY8cBFZrrYvFiub5590ysP5xEFitZX2xWJb274fly91vYL3nILDaGh7Ot0+gE8ePJ30XDgPrJQeB1dJCppLOy+RkMqTVrFccBFZLW7fmXcH8HD3qi86sdxwEVivDw8n6wmU0NgYrV+ZdhVWRg8BqodmEt761PKeDZvLqqx5NZNlzEFjltS4WO3Ik70qysX+/O48tWw4Cq7Th4WJdLJaVrVthxQoHgmXDQWCVNTxc/lNBszl6NJkl1WFg3XIQWGVVOQRaWqudeYEb64aDwCqn2UxOm9TJzp0OA1s4B4FVSrMJ27fnN5V0nm69Ne8KrKwcBFYZrdFBx47lXUk+ItwqsIVxEFglVHV00Hz5FJEthIPASq/qo4Pmy2Fg8+UgsFLbuNEhMJ2dO2HRIgeCdcZBYKU1OJhcZWvTi0gCwZPV2VwcBFZKzWYyCZvNbWzMLQObnYPASmdwsHzTSOdt505fgWwzyyQIJF0u6UlJByTdOM3jN0t6NP36nqRX2x6bbHtsXxb1WHUND7slsFCejsJmsqTbF5C0GPgC8D7gOeBhSfsi4mdnbyPiP7Xt/2fAu9pe4qcRcWm3dVg97NqVdwXl1ZqO4pvfhJGRvKuxIsmiRXAZcCAino6IE8DdwOZZ9v8Q8OUM3tdqZHgYlizxdQJZ8GkimyqLIFgL/LDt/nPptjNIWgdcDHyjbfNySeOSHpJ05UxvIml7ut/4oUOHMijbyqJ1nYBDIDteA9na9buzeAtwb0S0/0mvi4gB4E+Az0v6pemeGBG7ImIgIgbWrFnTj1qtAHyxWG8cPer1DOwNWQTBQeCitvsXptums4Upp4Ui4mD6/WngQd7cf2A1VZWlJYvs6NGkz8BhYFkEwcPAJZIulrSU5MP+jNE/kt4BrAS+1bZtpaRl6e3VwHsAXyJUc60ZRKuytGTROQys6yCIiFPAx4D7gCeAeyJiQtKnJV3RtusW4O6IiLZt7wTGJX0HeAD4bPtoI6unm26q7wyiedm61Vcg15ne/LlcDgMDAzE+Pp53GdYD7hPI19CQh5ZWmaRH0j7ZN/GVxVYYDoH8+fjXk4PACsMrbBWD5LmJ6sZBYLlqNmH16uTDp4RnKSvLs5bWi4PActNsJp2Uhw/nXYlNZ2zMo4nqwkFgufHVrcV39dUOgzpwEFgums3kgiYrtgjPWloHDgLru2YTrrsu7yqsU6dPw44deVdhveQgsL5q9QucPJl3JTYfR454DeQq63o9ArNObdzoNYbLrLUG8ve+B6OjeVdjWXKLwHqu2YSzznIIVIVHE1WPWwTWU61TQVYtrZ9po5FvHZYNtwisp66+Ou8KrFe2bnWfQVU4CKxnBgd9tXDV+QrkanAQWE8MDyfnkq36xsbcMig79xFY5jw6qH527oT3vMd9BmXlFoFlptmEZcscAnXlK5DLy0FgmWiNDjpxIu9KLC+nT3uls7JyEFjXPETU2o2NOQzKxkFgXWk2Ydu2vKuwohkbS9aZ8KmicnAQWFduuAEmJ/Ouworo8GH48IcdBmXgILCueFEZm82JE3DTTXlXYXNxENi8tS8vaTaXZ5/1dQZF5yCwjrUCwMtL2nzt3Alr1+Zdhc3EF5RZRzwyyLr1/POwciW88kreldhUbhFYRzx5nGXh1Vdh+XJ3IBeNg8Dm5MnjLEvHjydLlToMisNBYLNqNj15nGXv5MlkaKkVg4PAZuWLxaxXTpyAt7zFLYMiyCQIJF0u6UlJByTdOM3j10o6JOnR9OsjbY9tk/RU+uWPnQIZHPTFYtZbr78O117rMMhb10EgaTHwBeCDwAbgQ5I2TLPrVyLi0vTr9vS55wKfAn4TuAz4lKSV3dZk3WkNE/UpIeuHU6dgx468q6i3LFoElwEHIuLpiDgB3A1s7vC5HwDuj4iXI+IV4H7g8gxqsgUaHvZ1AtZ/R47AihVuGeQliyBYC/yw7f5z6bap/kDSdyXdK+mieT4XSdsljUsaP3ToUAZl21TNZnLhj1kejh71aKK89Kuz+B+A9RHxb0j+698z3xeIiF0RMRARA2vWrMm8QHPHsOXv5MlkIkPrryyC4CBwUdv9C9NtPxMRhyPieHr3duDfdvpc649m0x3DVgyHD8Nb3+qWQT9lEQQPA5dIuljSUmALsK99B0nnt929AngivX0f8H5JK9NO4ven26yPPH2EFc2RI17trJ+6DoKIOAV8jOQD/AngnoiYkPRpSVeku/25pAlJ3wH+HLg2fe7LwF+RhMnDwKfTbdYHrTWGHQJWVF7trD8UJZw7YGBgIMbHx/Muo9SazWSx8dOn867EbG5790KjkXcV5SfpkYgYmLrdVxbXkEPAymbrVvcZ9JKDoEb
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 432x432 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAF2CAYAAABj+Z+GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmZklEQVR4nO3df5Ac5X3n8fdHK68wkhPrVxF+eSUnutgiuRK5CfGdq5KLWWPsPxC5kERkBcKQEtLGia9SSRkXV+UrEi44qTrs3HkROiwQaMvY4cplpRIXhwRc/jGEVQUDEoclyxiQsVkL7ETCSNbu9/7oHtQazf6cnumZ7s+rampnunt6vppdzWf6eZ5+WhGBmZlV14KiCzAzs2I5CMzMKs5BYGZWcQ4CM7OKcxCYmVWcg8DMrOIcBGZmFecgMCuIpFWSvijpsKSfpD//UtI7i67NqmVh0QWYVdj7gD5gK3AQeD+wHVgObC6wLqsY+cxiM5D0OPD/gBPA9enie4BPRcRkB+sYBv48IpZ36jXN3DRkdtoQyf+Jfw/cTPKt/D9PtbGk90g6NsNt2xxr+Bngjfn+A8zmw0cEZrx9RHAB8IuR/qeQ9F+ALRFx0RTPWQismmHX/xIRr82yhgFgH/DfIuK/z7J0s5Y5CMx4OwheiojrM8suB/YAPxsR/9Lm1z8PeBx4Gvj98H9M6yA3DZnNU15NQ5J+DngMeA64ziFgneZRQ2an/ZokZT6IPwB8b5qjge8B62bY57RHEpLOJwmB/cC1EXFqDvWa5cJBYHbaBcDnJI0Avwz8GfAXU22cfmgfmu+LSbqApDnoeySd0isk1VePR8TEfPdtNhcOArPTRknG9T8JBPBF4M42vt4VwJr09lLDutXAi218bbO3ubPYjLc7i5+LiE8UXYtZp7mz2Mys4hwEZmYV56YhM7OK8xGBmVnFOQjMzCquJ4ePrlixIlatWlV0GWZmPWXfvn0/jIiVjct7MghWrVrF2NhY0WWYmfUUSd9tttxNQ2ZmFecgMDOrOAeBmVnFOQjMzCrOQWBmVnEOAjOzinMQmJlVnIPAzKziHARmZhWXSxBI2iHpNUnPTbFekv5G0iFJz0j6lcy6TZIOprdNedRjVrTRUVi1ChYsSH6OjhZdkdnU8joiuA+4cpr1H+X0Jfk2A3cBSFoGfAb4NeAy4DOSluZUk9mM6h/YEvT1JT8lWLEiWTc8nHyY15fP9rZxI3z3uxCR/Ny4cebnODCsKLkEQUT8I/D6NJusB+6PxBPAuyWdD3wEeCQiXo+IN4BHmD5QzOZtdDT5gG/2gQ0wOXl626NHk3V33ZV8mHfCbAKjry8JJ7M8daqP4ELg5czjV9JlUy0/i6TNksYkjY2Pj7etUOtd2eaYFSvO/ia/cWPyAd/LJieTcMr+uxYudDhYa3qmszgitkdELSJqK1eeNYuqVdzw8JnNMUePdu6bfNEmJk6Hw6JFZzdvmc2kU0FwBLg48/iidNlUy82ayn7rX7LkdLv+XXcVXVl3OHny9P1685aDwWbSqSDYDVyfjh76APDjiHgVeBi4QtLStJP4inSZ2RkGB8/uhD1+/Mx2fZtePRiWLDnddOYOaoP8ho9+CfgG8IuSXpF0k6Qtkrakm/wDcBg4BPwvYBggIl4H/hx4Kr3dli6zihseTtq+699m9+4tuqLyOH78dNNZvYN60SIHQpUperAhtVarha9QVi6XXAIHDhRdhUFypHDzzTAyUnQlljdJ+yKi1ri8ZzqLrVyybf3veIdDoJs0jkxy81H5OQis40ZH4brrTrf1nzpVdEU2nXrz0YIFHqZaVg4Ca7vGE7k2bqzO0M4yiUiOFC65pOhKLG8OAmurwcFynMhlpx044BPaysZBYG0zOOjRPlVQP6HNRwq9y0FgucpO4uYQqJb6kYIDofc4CKwlo6Pwrnc1n8TNqqkeCO5c7h0OApu30VG4/no4dqzoSqwb1TuXzznHw0+7nYPA5iQ7AmjjRk/xYDM7cQI2bXIYdDMHgc2aRwDZfE1MJOeOeH6j7uQgsBnVr9Llzl9rRcTZV2wbHCy6KgMHgU1jdDSZ/qGTV+myatm712HQDRwE1lR9GghP/2DttnevRxcVzUFgb8teqN3TQFgn1Se588VziuEgMCAJATcBWdHqF8/x9RE6y0FQcfXhoL7Uo3WTkyeTQHAYdIaDoKLqAeDhoNbNNm70cNNOcBBU0PBw0hHsALBeUB9u6tFF7eMgqJjRUfcFWG/auzeZ18pHB/lzEFREfVbQjRuLrsRs/o4dg49/3GGQNwdByY2Owjvf6VlBrTx++lP3HeTNQVBiw8PJf5i33iq6ErP8ffe7SV+XT0ZrnYOgpOp9AWZlVp/q2h3JrXEQlEz9QjHuC7Aq2bs3OSveTUXz4yAoEV8oxqosIvkC5KaiuXMQlES9P8AXirGqu+suWLjQRwdz4SAogfo8QWaWqF8Ix2EwOw6CHudOYbPmInyJzNnKJQgkXSnpBUmHJN3SZP2dkp5Ob9+S9KPMuonMut151FMF2bmCzKy5iQnYvNlhMBNFi3MNSOoDvgV8GHgFeAq4NiIOTLH9HwGXRsSN6eNjEbFkLq9Zq9VibGyspbp7mZuCzOZm8WIPogCQtC8iao3L8zgiuAw4FBGHI+Ik8CCwfprtrwW+lMPrVpKbgszm7vhxX/hmOnkEwYXAy5nHr6TLziJpAFgNPJpZfI6kMUlPSLp6qheRtDndbmx8fDyHsnvPJZe4KcisFUePwo03OgwadbqzeAPwUERMZJYNpIcqvw98TtLPN3tiRGyPiFpE1FauXNmJWrvG6GgyHO5A08Y2M5uLkyfh1luLrqK75BEER4CLM48vSpc1s4GGZqGIOJL+PAw8DlyaQ02lUT9JbGJi5m3NbHY8AeOZ8giCp4A1klZL6if5sD9r9I+k9wFLgW9kli2VtCi9vwL4IODvvRk33OCTxMzawc1Dp7UcBBFxCvgE8DDwPPCViNgv6TZJV2U23QA8GGcOU3o/MCbpm8BjwB1TjTaqouFhOHWq6CrMymnjRncg17U8fLQIVRg+6iGiZp2zYAHcfz8MDRVdSXu1c/io5cwhYNZZk5Nw881FV1EcB0EX2rat6ArMquf4cejrq+bspQ6CLjI6CkuW+MLyZkWZnEyOxqsWBg6CLlEfJnr8eNGVmFnVmmYdBF1gdDSZJdHDRM26R5Uuf+kgKNjwcDJvuk8YM+sue/cWXUHnOAgKNDqadAy7T8CsO0nV6C9wEBTo5psdAmbdrgqdxw6CggwOumPYrFeUfUi3g6DDRkdh0aJqtT+a9bqI5Ozjsh4ZLCy6gCoZHfX1BMx6VcTpYaUjI8XWkjcfEXTQJz9ZdAVm1qoyNhM5CDro6NGiKzCzVkWUb7ZSB0EHDA8nc5iYWTlcf325wsBB0GaDg0m7os8aNiuPycmkv68s1zJwELTR6KhHB5mV2dGjsHlz74eBg6CNfIFss/J7883e/7/uIGgjXyDbrBpeeqnoClrjIGiTXj9UNLPZi+jtk80cBG1S5cvemVVRL89J5CBog+FhzyNkVkW9ekEbB0GORkdh4cLe/WMws9ZdcknRFcydgyAn9XmEfIEZs2o7cKD3+ggdBDm54YaiKzCzbnHjjUVXMDcOghwMDsKpU0VXYWbd4uTJ3joqcBC0yGcPm1kzGzf2Thg4CFrUa4eAZtY5vTKM3EHQgtHR5BDQzKyZXhlGnksQSLpS0guSDkm6pcn6GySNS3o6vf1BZt0mSQfT26Y86umUXp9fxMwMcggCSX3AF4CPAmuBayWtbbLplyNiXXq7J33uMuAzwK8BlwGfkbS01Zo6YXTUcwmZ2cze+c7u7yvI44jgMuBQRByOiJPAg8D6WT73I8AjEfF6RLwBPAJcmUNNbTU87GsPm9nsvPUWXHddd4dBHkFwIfBy5vEr6bJGvy3pGUkPSbp4js/tGqOjPnPYzOYmArZsKbqKqXWqs/jvgFUR8W9JvvXvnOsOJG2WNCZpbHx8PPcCZ8sXoDez+Th2rOgKppZHEBwBLs48vihd9raIOBoRJ9KH9wD/brbPzexje0TUIqK2cuX
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 432x432 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAF2CAYAAABj+Z+GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAldklEQVR4nO3dfZAc9X3n8fdHq5MwkmP0VJgnrSDmHIvLlWxPcc75KnEs2cZUHSIXYsu3MnIgJViFxHeppIxLV5UrEtXh3B/YubOEVIRYoC2wQ8pl5WIfhwRc7qoMYVWRAYkTWuMIJLARAhwjWZKRvvdH91it1czuzk7PdM/051U1tTO/fpjv9s72Z7p//aCIwMzMqmtG0QWYmVmxHARmZhXnIDAzqzgHgZlZxTkIzMwqzkFgZlZxDgIzs4pzEJgVRNIMSdslvSjpuKRXJG2TdEnRtVm1OAjMivUo8CngvcBvAlcA3yy0Iqsc+cxiM5D0OPD/gBPAjWnzPcAXIuJ0F+u4DvgW8I6ION6t97Vq8xaB2RlDJP8TvwLcAqwF/kOzkSUtlvTWJI+7p/rmkuanNTzpELBu8haBGT/fIrgYeG+k/xSS/hNwa0Rc2mSamcCSSWb9TxHx6iTv/SXgNuB84Ang30bEay39AmZtcBCY8fMgeDEibsy0LQd2AO+KiH/q4HsvBOYDg8AfA28Bnwz/c1qXzCy6ALNeJWkxsHeS0bZFxK0TjZB++38NeF7Sc8BLwL8B/k8uhZpNwkFgdsa/kqTMN/EPAS9PsDXwMrBsknm2uiVR77eb3eJ0ZtPmIDA742Lgy5I2Ar8M/BHwp81Gjoi3gbHpvpmkXwE+APxf4E3gF4E/Af4xbTPrCgeB2RkjwADwJBDAXwB3dfD9fgrcANwBzAFeAf4n8GkfNWTd5M5iM37eWfxsRNxWdC1m3ebzCMzMKs5BYGZWcd41ZGZWcd4iMDOrOAeBmVnF9eThowsXLowlS5YUXYaZWU/ZtWvXaxGxaHx7TwbBkiVLGB0dLboMM7OeIulAo3bvGjIzqzgHgZlZxTkIzMwqzkFgZlZxDgIzs4pzEJiZVZyDwMys4hwEZmYV5yAwM6u4XIJA0r2SXpX0bJPhkvTnksYkPS3pA5lhayTtTx9r8qjHrFNGRmDJEpgxI/k5MjL9cdetg5kzQUp+rlvXubrNJhQRbT+AXyW59+qzTYZfC3wHEMkNwZ9M2+cDL6Q/56XP5032fh/84AfDrG7btoi5cyMgeUgRy5dHDA4mz+fMOTMs+xgePnseg4NJ+4wZjcdv9njHO1obfzqPek0LFiQPKal3eDh53WiaBQuS32v871ifNjvMqgEYjQbr1NzuRyBpCfA/IuJfNBi2GXg8Ih5IX+8DPlJ/RMQtjcZrplarha811H/WrYNNm4quotoGBuC974V9++DUqeT12rWwcWPRlVkeJO2KiNr49m71EVwCvJR5fTBta9Z+DklrJY1KGj18+HDHCrXOGL+b5Kqrkl0i2YdDoHinTsHevcnP+utNm5K/T3bXVSu7yKz8eqazOCK2REQtImqLFp1zFVUr2MgILFx49op94UJYsSJZWaxeDQcOJDstDhxIVjbWW+qBIJ3791y9+sww93f0nm4FwSHgsszrS9O2Zu1WcuNX/KtXw5EjZ49z5Ajs3JmsLKw6slsR2S8F3moor24FwXbgxvTooQ8BP46IV4CHgY9LmidpHvDxtM1KZiorfrNmjhxJPjMLF3p3UhnlcmMaSQ+QdPwulHQQ+GPgnwFExN3At0mOHBoDjgG/nQ57XdKfAE+ls7ojIl7PoyZrz7p1sGXLmX3FZnmof3k4cABuuil5PjRUXD2WyO2ooW7yUUOdMTKS/HOePFl0JVZ1CxbAV77ikMhb0UcNWUmNjMDcuWd29zgErAzqu5Kk5PPp3Uid5SCokPohf9lOvNWr4ejRoisza+7o0eRz6iOROsdBUBErVpw55M+sF40/EmnGDIdDXhwEfSx7LZudO4uuxixfEUk4rFhRdCW9z0HQZ7L7/Ddt8lE/1v927vS5Cu1yEPSR+u4f7/O3Kqp3MHsLoXUOgj5x1VXe/WMGyf/BVVcVXUVvcRD0qHXrkitD1jvOfO0eszP27j3zv+GzmCfnIOhB9cs1nz5ddCVm5Ve/KJ63EppzEPSI7DkAvlyzWevqWwnuVD6Xg6AHjIwkNwfxOQBm7at3KvuM5TMcBCW3bl3yoT12rOhKzPqLz1g+w0FQUiMjyclg3g1k1lmbNnnrwEFQMiMj8M53Jt9UfDKYWXccPQqf/Wx1wyCX+xFYPkZGkg9jD14Z3KznRSRfwKB6l7/2FkGJ3HKLQ8CsaKtXJ1vlVdo6cBCUxMiILw1hVhZvvVWtXUUOgoLVO4Xrm6RmVg4RZ26n2e/cR1CQSy6Bl18uugozm8jJk8n/6qFDRVfSWd4iKIBDwKx3vPxy/1+ewkHQZSMjDgGzXrN3b39fvM5B0EX1s4TNrPfUL17Xj2ciOwi6pH7FUDPrbZs29V8YOAi6YGTEIWDWTzZt6q/dRA6CLlizpugKzCxvq1f3Txg4CDpsxQpfM8isX61endzjoNfvk+wg6JD6iWK+j7BZ/9u5s7fDwCeUdcCKFQ4As6rp5f/5XLYIJF0jaZ+kMUm3Nxh+l6Td6eN5SW9mhp3KDNueRz1FWreutz8QZjZ9vXo0UdtbBJIGgK8CHwMOAk9J2h4Re+vjRMR/zIz/e8D7M7P4aUQsa7eOstiypegKzKwo9aMDN24sto5W5bFFcDUwFhEvRMRJ4EFg5QTjfwZ4IIf3LSV3DJtVWy8eWppHEFwCvJR5fTBtO4ekQeBy4NFM83mSRiU9Ien6HOopTK/98c2sM9avL7qC1nT7qKFVwEMRkf3ePBgRNeDfA1+W9IuNJpS0Ng2M0cOHD3ej1paMjMDatUVXYWZlcOBAb30xzCMIDgGXZV5fmrY1sopxu4Ui4lD68wXgcc7uP8iOtyUiahFRW7RoUbs15+7WW+HYsaKrMLOyWL26dw4pzSMIngKulHS5pFkkK/tzjv6R9EvAPOC7mbZ5kmanzxcCHwb2jp+2rEZGkisSSskdjczMsnrl/IK2gyAi3gZuAx4GngO+ERF7JN0h6brMqKuAByPOuivv+4BRSd8DHgPuzB5tVGYjI3DjjckmoJlZMzt3ln83kaIH75Zeq9VidHS00Bre8Q44frzQEsysR8yYAffdB0NDxdYhaVfaJ3sWX2JimhwCZjZVp08nB5OUdcvAQTANvXr2oJkV59ix8h5W6iCYBt9bwMymo6x9ig6CFvX7TazNrLPKuA5xELRg3brkJtZmZtNVxnWIg2CKfM9hM8tL2TqNHQRT4BAwszx9/vNFV3A2B8EUOATMLE9HjhRdwdkcBJMo2yacmfWHuXPLs35xEEzilluKrsDM+tHRo/DZz5YjDBwEE1i3LvljmZl1QkQ5+gscBBPYvLnoCsys35Whv8BB0MTISHJ9EDOzfucgaMJ9A2ZWFQ6CJtw3YGbdUvSFLB0EDRT9RzGzatm0qdg7mTkIxhkZ8QlkZtZ9Rd7W0kEwThkO5TKzatq5s5j3dRCMU4ZDuczMuslBYGZWIkX0UToIzMxK5O67u/+eDoKMMlzzw8yqLaL77+kgyPBJZGZWRQ6C1IoVPonMzMqh23snHAQkC72ow7bMzMZbs6a77+cgANavL7oCM7MzTp3q7tFDDgLgxReLrsDM7GybNnVvF5GDAJg/v+gKzMzO1a0DWHIJAknXSNonaUzS7Q2Gf07SYUm708fvZIatkbQ/fXR5z1jiJz8p4l3NzCbWrQNYZrY7A0kDwFeBjwEHgackbY+IveNG/XpE3DZu2vnAHwM1IIBd6bRvtFtXK06e7Oa7mZmVSx5bBFcDYxHxQkScBB4EVk5x2k8Aj0TE6+nK/xHgmhxqmrKrrurmu5mZlU8eQXAJ8FLm9cG0bbzflPS0pIckXdbitB2zd/x2i5lZiXTj0tTd6iz+G2BJRPxLkm/9W1udgaS1kkYljR4+fDj3As3Myqgb5zjlEQSHgMsyry9N234uIo5ExIn05T3AB6c6bWYeWyK
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 432x432 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"for n in range(1, 4):\n",
|
|||
|
|
" x, y = numba_unit_disk(n)\n",
|
|||
|
|
" \n",
|
|||
|
|
" plt.figure(figsize=(6,6))\n",
|
|||
|
|
" plt.scatter(x, y, color=\"blue\")\n",
|
|||
|
|
" plt.title(\"p = {}\".format(n), fontsize=fontsize)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 4,
|
|||
|
|
"id": "43ba0e12",
|
|||
|
|
"metadata": {
|
|||
|
|
"scrolled": true
|
|||
|
|
},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"759 ms ± 74.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
|
|||
|
|
"6.3 ms ± 1.44 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"%timeit unit_disk()\n",
|
|||
|
|
"%timeit numba_unit_disk()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "0b618b90",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__Вывод:__\n",
|
|||
|
|
"- JIT крут, вычисления в 150-200 раз быстрее"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "9aff9c18",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"### Задача 6"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 92,
|
|||
|
|
"id": "74c4965f",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"c = 3"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 93,
|
|||
|
|
"id": "41852c3b",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"# Generate some data\n",
|
|||
|
|
"np.random.seed(42)\n",
|
|||
|
|
"lambda1 = np.random.normal(size=(c, c))\n",
|
|||
|
|
"lambda2 = np.random.normal(size=(c, c))\n",
|
|||
|
|
"lambda3 = np.random.normal(size=(c, c))\n",
|
|||
|
|
"G1 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
"G2 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
"U = np.random.normal(size=(c, c, c, c))"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 94,
|
|||
|
|
"id": "c72140e1",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def Z_naive(lambda1, lambda2, lambda3, G1, G2, U):\n",
|
|||
|
|
" c = lambda1.shape[0]\n",
|
|||
|
|
" Z = np.zeros(shape=(c, c, c, c))\n",
|
|||
|
|
" for a, b, c, d, e, f, g, h, i, j in itertools.product(*([range(c)]*10)):\n",
|
|||
|
|
" Z[a, h, i, j] += lambda1[a, b]*lambda2[d, e]*lambda3[g, h]*G1[c, b, d]*G2[f, e, g]*U[i, j, c, f]\n",
|
|||
|
|
" return Z\n",
|
|||
|
|
"\n",
|
|||
|
|
"def Z(lambda1, lambda2, lambda3, G1, G2, U):\n",
|
|||
|
|
" return np.einsum(\"ab,cbd,de,feg,gh,ijcf->ahij\", lambda1, G1, lambda2, G2, lambda3, U, optimize=True)\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 95,
|
|||
|
|
"id": "a4ba9f23",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"615 µs ± 144 µs per loop (mean ± std. dev. of 2 runs, 10 loops each)\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# %timeit -n 2 -r 2 Z_naive(lambda1, lambda2, lambda3, G1, G2, U)\n",
|
|||
|
|
"%timeit -n 10 -r 2 Z(lambda1, lambda2, lambda3, G1, G2, U)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 96,
|
|||
|
|
"id": "19913a61",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"(3, 3, 3, 3)\n",
|
|||
|
|
"(3, 3, 3, 3)\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"print(Z_naive(lambda1, lambda2, lambda3, G1, G2, U).shape)\n",
|
|||
|
|
"print(Z(lambda1, lambda2, lambda3, G1, G2, U).shape)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 97,
|
|||
|
|
"id": "860ebc25",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"4.99716543697806e-26\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"print((np.sum(Z(lambda1, lambda2, lambda3, G1, G2, U) - Z_naive(lambda1, lambda2, lambda3, G1, G2, U))**2))"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 98,
|
|||
|
|
"id": "15b20588",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"['einsum_path', (0, 1), (0, 1), (0, 3), (1, 2), (0, 1)] Complete contraction: ab,cbd,de,feg,gh,ijcf->ahij\n",
|
|||
|
|
" Naive scaling: 10\n",
|
|||
|
|
" Optimized scaling: 6\n",
|
|||
|
|
" Naive FLOP count: 3.543e+05\n",
|
|||
|
|
" Optimized FLOP count: 2.431e+03\n",
|
|||
|
|
" Theoretical speedup: 145.740\n",
|
|||
|
|
" Largest intermediate: 8.100e+01 elements\n",
|
|||
|
|
"--------------------------------------------------------------------------\n",
|
|||
|
|
"scaling current remaining\n",
|
|||
|
|
"--------------------------------------------------------------------------\n",
|
|||
|
|
" 4 cbd,ab->acd de,feg,gh,ijcf,acd->ahij\n",
|
|||
|
|
" 4 feg,de->dfg gh,ijcf,acd,dfg->ahij\n",
|
|||
|
|
" 4 dfg,gh->dfh ijcf,acd,dfh->ahij\n",
|
|||
|
|
" 5 dfh,acd->acfh ijcf,acfh->ahij\n",
|
|||
|
|
" 6 acfh,ijcf->ahij ahij->ahij\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"print(*np.einsum_path(\"ab,cbd,de,feg,gh,ijcf->ahij\", lambda1, G1, lambda2, G2, lambda3, U, optimize=True))"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "87ec392f",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__np.einsum__ вычисляет свёртку тензоров, причём разбивает суммирование на несколько этапов:\n",
|
|||
|
|
" \n",
|
|||
|
|
"- вычисляет части свёртки по отдельности, создаёт промежуточные массивы(тензоры)\n",
|
|||
|
|
"- вычисляет свёртку от новых тензоров"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "2359e3ae",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__Наивный способ__\n",
|
|||
|
|
"\n",
|
|||
|
|
"Для $Z_{ahij}=\\underset{bcdefg}{\\sum}\\lambda^{(1)}_{ab}\\Gamma^{(1)}_{cbd}\\lambda^{(2)}_{de}\\Gamma^{(2)}_{feg}\\lambda^{(3)}_{gh}U_{ijcf}$ с размерностью индексов $\\chi$ есть 6 индексов суммирования, т.е. $Z_{ahij}$ вычисляется за $\\chi ^ 6$, весь тензор за $\\chi ^ 6 \\cdot\\chi ^ 4$"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "10dd2839",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__Оптимизированный способ__\n",
|
|||
|
|
"\n",
|
|||
|
|
"Делаются свёртки \n",
|
|||
|
|
"\n",
|
|||
|
|
"- $\\lambda^{(1)}_{ab}\\Gamma^{(1)}_{cbd} = A_{acd}$ за $\\chi$, т.е. $A$ за $\\chi^4$\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"- $\\lambda^{(2)}_{de}\\Gamma^{(2)}_{feg} = B_{dfg}$, за $\\chi$, т.е. $B$ за $\\chi^4$\n",
|
|||
|
|
"\n",
|
|||
|
|
"На этом моменте у нас есть два новых тензора $A, B$, которые мы вычислили за $2 \\chi^4$\n",
|
|||
|
|
"\n",
|
|||
|
|
"- $B_{dfg}\\lambda^{(3)}_{gh}=C_{dfh}$ за $\\chi$, т.е. $C$ за $\\chi^4$\n",
|
|||
|
|
"\n",
|
|||
|
|
"Добавляется ещё один тензор $C$, который мы вычислили за $\\chi^4$, так как тензоры $B, \\lambda^{(3)}$ были уже известны. На данный момент $3 \\chi^4$ операций\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"- $C_{dfh}A_{acd}=D_{acfh}$ за $\\chi$, т.е. $D$ за $\\chi^3$\n",
|
|||
|
|
"\n",
|
|||
|
|
"Новый тензор $D$, суммарно $3 \\chi^4 + \\chi^5$ операций\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"- $D_{acfh}U_{ijcf} = Z_{ahij}$ за $\\chi^{2}$, т.е. ещё $\\chi^{6}$ операций\n",
|
|||
|
|
"\n",
|
|||
|
|
"В итоге вычислили тензор $Z$ за $O(\\chi^{6})$ операций, наивный способ совершает $O(\\chi^{10})$ операций. Оптимизация сложности происходит за счёт дополнительного расхода памяти, упрощения выражения, поэтому для небольшого числа операций лучше использовать обычное суммирование."
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "2aafecb8",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"Для $\\chi = 1$ получается (x0.02)\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"naive: 3.59 µs ± 220 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)\n",
|
|||
|
|
"\n",
|
|||
|
|
"optimised: 225 µs ± 3.51 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "92128c85",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"Для $\\chi = 4$ получается (x4000)\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"naive: 1.59 s ± 25.5 ms per loop (mean ± std. dev. of 2 runs, 2 loops each)\n",
|
|||
|
|
"\n",
|
|||
|
|
"optimised: 401 µs ± 17.5 µs per loop (mean ± std. dev. of 2 runs, 100 loops each)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "db7a9dcc",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"Для $\\chi=50$: \n",
|
|||
|
|
"\n",
|
|||
|
|
"optimised: 438 ms ± 17 ms per loop (mean ± std. dev. of 2 runs, 10 loops each)\n",
|
|||
|
|
"\n",
|
|||
|
|
"Для $\\chi=5-50$ naive заexplorить, к сожалению, не получается "
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 326,
|
|||
|
|
"id": "d332b860",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAEQCAYAAACutU7EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuy0lEQVR4nO3dd5gV9dn/8fdNx4CogICPItifiIqKgp1g7NHoZSMaaxRje0wz+SnmsST2xPJENGIK2KJGo1GKiCVCgFWKDUFQsFOVIh2WvX9/fGfleLLl7O6cmVM+r+uaa/bMzM7ce+R8vc+3mrsjIiIiUkyapR2AiIiISEMpgREREZGiowRGREREio4SGBERESk6SmBERESk6CiBERERkaKjBEYSZ2bDzOyjtOMQkeJgZh+Z2bBCv2cOz+xvZm5m/ZN8bqlSAiO1MrNzow/bOjPrXsP5EUpERMqLmXU3sz9GCcA6M1tkZk+b2YFNvO+xZnZdTGFKGWiRdgBSFFoBVwM/jul+F6LkWaTomNlBwKjo5QPATKAbcC7wbzO7wt3/0MjbHwtcClxXw7ldgapG3rc2+binJEgJjOTiTeA8M7vJ3T9p6s3cfUPTQxKRJJnZlsCTwBrgIHefk3Hu98AY4C4zm+ruE+N8truvi/N++bqnJEvfgiUXN0f7wXVdFDU5vWhmC6Kq5ffN7Coza5Z13Tf6wJjZdDMbX8s93zezf2W8NjO73MzeMbO1UfX1n82sU6P/OhHJxUVAV+DKzOQFwN3XAOcADvxv9fGMPh9nmtn1ZjbPzFab2Rgz2znjumGE2hei66u3HtGxb/RXybjvGWZ2rZl9bmYrzewfZralmbUyszvMbGF0fLiZtc2MuaY+MGZ2SVS2rDKz5Wb2ppldlHVNNzP7U0Y5N9PMLs5+s8xsWzN7JrrXIjO7E2jdgPdb6qEaGMnFJ8BfgB+Z2Y111MJcSqhSHgWsBQ4HbgI6AP+vjvs/BtxgZtu6+2fVB81sX2An4PcZ194H/AgYDtwDbAdcDuxvZvu5+9pG/H0iUr/jCZ/rJ2o66e4fmtm/gQFm1jZKaqr9CmgO/A7YErgCeMXM9nT3JcD9wDbAEcBZGb+3uJ6YfhnFdAuhrLgc2AisBroA1wN9gbOBj8lIrrKZ2Y+AIYRapnuAlsDuwIFRfJjZ1kBF9LfcCywilHP3mllHd/9tdF1b4CWgO/B/wDzgTGBAPX+PNIS7a9NW40Zo13agHyFRWAfcn3F+BPBRxuvNarjHUGAl0Drj2LCs39spes5Ps373NmAD0Cl6fWB03dlZ1x0cHR+U9numTVupbsBS4M16rrk7+izuEb3uH71eCGyRcd2A6PhvM47dE/6XVON9PwKGZbyuvu8MoFXG8UcJ/VrGAJZxfCLwWT33fBqYXs/fNxRYAHTOOv4AIWnaInp9RRTfaRnXtAXei473T/u/ZylsRdWEZGaHmtmzUXWhm9m5jbjH3WY2JWp++KiG89dlVWFmblvH8XcUI3f/lFALc56ZbV/LNasBzKx5VI3bCXgV+BawWx33/gCYCpyedeo04CV3/yLj9UrgeTPrVL0RCoWFwHca/QdK2UuifImu2cPMXjWzNdGz/tfMrKnxJ6A9sKKea6rPb551/EF3X1b9wt1fBt4FvtfEmB5y9/UZr18DDPirR1lDxvFtzKyuJpzlwLZmtl9NJ6P/RqcAIwHPKoNeICQofaPLjyWUSU9W/76HGqk/NeivkzoVVQIDtAOmE7LbNfVcW5tmhOaHB2s5/ztCr/rM7VXgX+6+qJHPLBU3Eb491NgXxswONrNxhG8iSwjVvw9HpzvUc+/HgL4Zbd4HANtHx6vtQvg3sDC6d+bWBSjbBFNikffyxcw2B8YS/g3vFz3rSuBnjXxeklYQkpi6VJ/PTnTer+Ha2UCPJsaU3Zy9PNp/WsNxA7ao4163EuJ+3czmWBgqnvmlqDOh+et8/rP8qW5Wqy6DtgfmuHv2KKfZdf410iBF1QfG3UcRDeHL7nwVHWsF/IbQ1rgVIcO/xt3HZNzj8ujaXwBH1vCMlYRv+dX33A44hG+2y5Yld//UzP4MXGBmN2aeM7MdgBcJH9CfEgqWtcA+hIKhvmT5cUKT0enR9acTmqyezrimGfAlMLCWeyxtyN8jkimJ8iX63c2Ac6Jv5NPNbDfgZ2Z2R1atQaGZCextZq299hE8exKafWtKWPJhYwOP11rT5e4zzWxXQu3JUYTaoYvM7F53v5RNZdjfCLXRNXm3/pAlLkWVwOTgr8COwBnAZ4R/iM9FnTvfauQ9f0T4H+NT8YRY9G4ivCfXZB0/gdDD/nh3/7j6oJn1zOWmUXI0ETjdzG4HTgWed/flGZfNIXTyq4gSTZEkxVG+HACM9292cB1DSIx6AB/GF27sRhDiP5VNNatfi2pPDwFezPr7AHbOvp5Qo/pRxuvUk7eoGfxJ4Ekza0Hor3eJmd1E6PuyAmjh7i/Wc6uPgb3MrFlWLcwueQi7bBVbE1KtzGxH4AeETlPj3H2uu99D+EZ1Ud2/Xes9mxOqCx+q4xtHWfEwSuhPhCGTmX1hqr/xfP0NJ2pvvqwBt38M2Jvwnm/DN5uPINTSNKOGkQTV/W4a8CyRnMVYvnQlNB9lWphxrpDdTxh1c1v2F5No1M1fCZ//G2r43bPNbIuM6wcQRviMzLhmVXQulc+xmXXMfO3ulcA70cst3H0jIbk50cz2quH3O2e8HEVo1j4l43xb4IK44y5npVQDsw/hwzMjqz9ca+DlRt7zaMLomweaFlrJuZlQC9OL8E0DwrfI9cAIM7uf8L6fRcNmuvw7cBdwB6EfzXOZJ919nJkNAa40sz2jZ64jjGI6hZDYDGvUXyRSt3yUL0XF3b80s+pOrG+YWfVMvF0JIxZ3Aq7wmiexWwhMiJqgtwB+AswnfNarTYn295jZaKASeM7dV8X/19ToBTNbBPybUNtSPSz7bcLfCWE6iP7ApOjvf5fQL6Y3cBLQJrruAcKXt+HRdBCfAz8klFcSk1JKYJoRqiD3I7TBZmpsh7xBwER3n9GUwEqNu39mZn8imngqOjbbzE4kNDHdBnxB6Mj4L0IP/Vzuu9DCpHWHA4/XVHC5+2VmNo2wrMGNhELuE0InurL4H4mkIq7yZQHhm3mmLhnnCpq7j4++PFxNGBXYjdBBdgJwnrtPqOVXbyVM3X8lIYEZD1zu7l9mXPMPwheYH0SbAT2JamYScD+hefAnhFFU8wh9XX5b3Qzk7ovMrC/wa+BE4GLCgIWZwM+rb+Tuq83scOAPhERmNfAIMBp4Ppk/p/RZYfcZq52ZrQQuc/dh0etdgFnAAHd/JYff/0X0+z1qOb8N4X+MF1Q/Q0TKQ77KFwsztt4KbO3RpItmdjXhy8C2Bd6Jt8EsrLr8CvADd89uEhZpkqKqgTGzdoRqPQjfiLqbWW9gSVQD8AgwzMx+DkwjjBToD8x1939E99iJMFxyG6BV9PsAM7LmEzifkPnXOOukiJSWhMqXR4Fro/v8ltCp8/8B15da8iKSb0VVA5ORzWcb7u7nmllLwhwlZwPbEqr2XicUDlOje/wLOKyGe/R094+iawyYC4x290vi/StEpBAlWL7sQZiyfn/CCMc/AjeUYgKjGhjJp6JKYEREpHgogZF8UgIjIiIiRadk5oERERGR8pFoJ14z60ZY9vxYwpoZc4GL3f3Vun6vU6dO3qNHj/wHKCINNnXq1C/cvXP9VxYulTEihamu8iWxBCaahXECYZKg4wgLYO1AmNmxTj169GDKlCn1XSYiKTCzj+u/qrCpjBEpTHWVL0nWwPwSmO/uZ2ccK+R1P0RERKRAJdkH5kTgNTN73MwWmdmbZnaZZc3LLSIiIlKfJBOYHYBLCP1ejgLuJvSHubSuXxIRERHJlmQTUjNgirtfFb1+w8x2JiQw92RfbGaDCGsR0b1798SCFBERkcKXZA3MfCB7UcSZQI3ZibsPdfc+7t6nc+eiHuAgIiIiMUsygZlAWI000y5A0Y9gEBERkWQlmcDcCfQzs8FmtpOZnQr8D2FNEBEpEFdeGTYRkdhNnQpHHgkzZzb5Vok
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 576x288 with 2 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAEQCAYAAACutU7EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqBElEQVR4nO3deXgUVb7G8e8vkS0E2WUJFxFw44IgE5HFERwVFBGHwW0QNRcVUPSK48JV1OvKqDACCurgdZAZiAvqMA4RNySiAmqC7CqoCApIQARkEbOc+0d1QhICdEinqrvzfp6nnnTX6VS9SQn8PHXqHHPOISIiIhJLEoIOICIiIlJeKmBEREQk5qiAERERkZijAkZERERijgoYERERiTkqYERERCTmqIARERGRmBNTBYyZnWlmr5vZBjNzZpZ2BMeYaGZZZvaLmX1bRvt9oWOXtR0TiZ9DREREKiamChggGVgB3AzsPcJjJADTgL8fpH0c0KzU9j6Q6ZzLOcJzioiISATFVAHjnHvDOXeXc+4VoKB0u5lVN7NHzex7M9tjZp+aWZ9Sx7jJOfcksPog59jlnPuhcAOqAb8Fnq2EH0lERESOwFFBB4iwqUAbYBDwPdAX+LeZneacW3qEx7wG+Al4NTIRRUREpKLipoAxszbAH4FWzrn1od2TzOwcYBhwwxEcMxEYAvzDObcvYmFFRESkQuKmgAE6AwasMrPi+2sA7x3hMc8D/gPdPhIREYkq8VTAJAAOOA3ILdV2pAN+hwILnHOrKhJMREREIiueCpjP8Hpgmjrn5lX0YGbWHLgAuLaixxIREZHIiqkCxsySgbahtwlASzPrBGxzzq02sxnA82Z2K7AYaAD0Ar5xzr0WOkZbvMexmwPVQ98PsMo592ux0w0BdgMvV+oPJSIiIuVmzrmgM4TNzHoBZfWuTHPOpZlZNWA0cBXQAtgGfALc75zLDh0jE+hZxjGOc859G/qMAd8Ac5xz5R78KyIiIpUrpgoYEREREYixiexEREREQAWMiIiIxKCYGMTbqFEj16pVq6BjiIiIiI+ys7O3Oucal9XmawFjZs2AR/Cm+K+DN1D2eufc+4f6vlatWpGVleVDQhEREYkWZrbuYG2+FTBmVg/4CPgQb36VLUBrQCs8i4iISLn42QNzB7DJOXdVsX1rfTy/iIiIxAk/B/H+HvjYzF4ysxwzW2JmN1qphYtEREREDsfPAqY13orQ3wB9gIl442FGlPVhMxtqZllmlrVlyxb/UoqIiEjU820iOzP7FchyznUvtm8MMMA5d/Khvjc1NdVpEK+IiEjVYmbZzrnUstr8HAOzCSi9qvPnwM0VPfDOnTvJyckhN7f0ItQSD2rXrk2LFi1ISNC0RSIi4vGzgPkIOLHUvhOAgz4iFY6dO3eyefNmUlJSqFWrFhpSE18KCgrYsGEDW7du5Zhjjgk6joiIhDgHf/4z9OsHp5zi//n9/F/a8UBXMxttZm3N7BLgv4HJFTloTk4OKSkpJCUlqXiJQwkJCTRp0oQdO3YEHUVERIpZswZGj4YPPwzm/L4VMM65T/GeRLoUWAE8DNwDPFWR4+bm5lKrVq0K55PoVa1aNfLy8oKOISIixWRkeF8vuCCY8/s6E69zLgPIiPRx1fMS33R9RUSiT0YG/Od/wrHHBnN+jYoUERGRcvn5Z5g/P7jeF1ABE/Puu+8+2rdvH3QMERGpQt55B3JzoW/f4DKogAlIWloaZsaDDz5YYn9mZiZmxtatW8M6zm233cb77x9yLUwREZGIysiAunWhe/fDf7ayqIAJUM2aNRk7diwVmWk4OTmZhg0bRjCViIjIwRUUwBtvQJ8+UK1acDlUwATorLPOolWrVgf0whTKz8/nmmuu4bjjjqNWrVocf/zxPPbYYxQUFBR9pvgtpLfffpvq1avz448/ljjOXXfdxSnFHtJfsGABPXv2JCkpiZSUFK6//np27txZCT+hiIjEmyVL4Icfgh3/AipgApWQkMAjjzzCM888w9dff31Ae0FBASkpKbz88st8/vnnPPzww4wZM4apU6eWebyzzz6bRo0aMXPmzKJ9zjnS09MZPHgwAMuXL6d3797079+fpUuX8tprr7FkyRKGDBlSOT+kiIjElYwMMIPzzw82h6+PUftl5EivQvRTp04wYUL5v69v37706NGD0aNH8+KLL5Zoq1atGg888EDR+1atWrF48WJeeOEFrrnmmgOOlZiYyOWXX86MGTMYPnw4AB999BHfffcdgwYNAmDs2LFcdtll3HrrrUXf9/TTT3PqqaeSk5Oj2W5FROSQMjKgSxdo3DjYHOqBiQKPPvooM2fOJDs7+4C2Z555htTUVBo3bkxycjLjx49n/fr1Bz3W4MGD+eijj1i3zluhYcaMGfTs2ZMWLVoAkJ2dzfTp00lOTi7aevToAVBmL5CIiEihLVvgk0+Cv30EcdoDcyQ9IUHq0qULAwcO5I477uCee+4p2v/SSy8xcuRIxo0bR/fu3Tn66KOZPHky//znPw96rM6dO3PSSSeRnp7ObbfdxsyZM3nssceK2gsKCrj22mu55ZZbDvjelJSUyP5gIiISV+bM8dZACvLx6UJxWcDEojFjxtCuXTvefPPNon0ffvghp59+OjfeeGPRvnB6SQYPHsyMGTNo3749u3fv5uKLLy5q69y5MytXrqRt27aR/QFERCTuZWRA06Zw6qlBJ9EtpKjRtm1bhg4dysSJE4v2nXDCCSxevJg5c+awZs0aHnzwwbDmfLniiitYtWoV99xzDxdeeCFHH310UduoUaP45JNPGD58OJ999hlfffUVs2fPZtiwYZXyc4mISHzIy4O33vJ6XxKioHqIgghS6N577+Woo/Z3ig0bNoxLL72UQYMGcdppp/Htt9+WGHx7MMceeyxnnHEGS5cuLXr6qNApp5zC/Pnz+fbbb+nZsycdO3bkzjvvpEmTJhH/eUREJH4sWAA7dkTH+BcAc84FneGwUlNTXVZWVpltn3/+OSeffLLPicRvus4iIsEaNQrGj4cff4Q6dfw5p5llO+dSy2pTD4yIiIgcVkYGnHmmf8XL4aiAERERkUNatw5Wroye20egAkZEREQOIyPD+xoNj08XUgEjIiIih5SRAW3awAknBJ1kPxUwIiIiclB798J773m3j8yCTrOfChgRERE5qHnz4Jdfomv8C6iAERERkUPIyIDataFnz6CTlORbAWNm95mZK7X94Nf5RUREpHyc8wqYc86BGjWCTlOS3z0wXwLNim0dfD6/iIiIhGnVKu8R6mi7fQT+FzB5zrkfim1bfD5/lZKZmYmZsXXr1godJy0tjX79+kUoVdn69etHWlpapZ5DRETKJxofny7kdwHT2sw2mtlaM3vRzFr7fP6os2HDBoYOHUqLFi2oXr06KSkpXHfddXz//fflOk6rVq0YN25ciX3du3dn06ZNNGzYsEIZJ06cyPTp0yt0DBERiT0ZGdCxI6SkBJ3kQH4WMB8DacB5wHVAU2CBmZX5r6uZDTWzLDPL2rIlPjtq1q5dS2pqKitWrGDatGl89dVXTJ8+nZUrVxYt3lgR1atXp2nTplgFn3urW7cu9erVq9AxREQktmzfDh99FJ23j8DHAsY5N8c597Jzbplz7l2gX+j8Vx/k81Occ6nOudTGjRv7FdNXI0aMICEhgXfffZezzz6bli1bctZZZ/Huu++SkJDAiBEjAOjVqxfDhw/n5ptvpn79+tSvX5/bb7+dgoKCovZ169Zx++23Y2ZFBUvpW0jPP/88ycnJzJkzh5NOOomkpCT69+/Pjh07eOWVVzj++OOpW7cuV155JXv37i3KWfoW0vz58+natSvJycnUrVuXLl26sGLFiqL2BQsW0LNnT5KSkkhJSeH6669n586dRe179uwhLS2N5ORkmjRpwpgxYyrvlywiIkfk7bchP18FzAGcc7uAlcDxQWUI0rZt23jzzTcZMWIESUlJJdqSkpK44YYbmDNnDj/99BMAM2bMoKCggIULF/LXv/6VKVOmMGHCBABee+01WrRowb333sumTZvYtGnTQc+7b98+/vKXvzB
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 576x288 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"N, O = [], []\n",
|
|||
|
|
"C = np.arange(3, 51)\n",
|
|||
|
|
"\n",
|
|||
|
|
"for c in C:\n",
|
|||
|
|
" np.random.seed(42)\n",
|
|||
|
|
" lambda1 = np.random.normal(size=(c, c))\n",
|
|||
|
|
" lambda2 = np.random.normal(size=(c, c))\n",
|
|||
|
|
" lambda3 = np.random.normal(size=(c, c))\n",
|
|||
|
|
" G1 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
" G2 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
" U = np.random.normal(size=(c, c, c, c))\n",
|
|||
|
|
" \n",
|
|||
|
|
" res = np.einsum_path(\"ab,cbd,de,feg,gh,ijcf->ahij\", lambda1, G1, lambda2, G2, lambda3, U, optimize=True)\n",
|
|||
|
|
" n, o = res[1].split('\\n')[3:5]\n",
|
|||
|
|
" n = np.float64(n.split(' ')[-1])\n",
|
|||
|
|
" o = np.float64(o.split(' ')[-1])\n",
|
|||
|
|
" N.append(n)\n",
|
|||
|
|
" O.append(o)\n",
|
|||
|
|
"\n",
|
|||
|
|
" \n",
|
|||
|
|
"plt.rcParams.update({'font.size': 14})\n",
|
|||
|
|
"plt.figure(figsize=(8,4))\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.subplot(121)\n",
|
|||
|
|
"plt.plot(C, N, color=\"blue\")\n",
|
|||
|
|
"plt.title(\"Naive\")\n",
|
|||
|
|
"plt.xlabel(\"$\\chi$\")\n",
|
|||
|
|
"plt.ylabel(\"FLOPs\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.subplot(122)\n",
|
|||
|
|
"plt.plot(C, O, color=\"red\")\n",
|
|||
|
|
"plt.title(\"Optimised\")\n",
|
|||
|
|
"plt.xlabel(\"$\\chi$\")\n",
|
|||
|
|
"plt.ylabel(\"FLOPs\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.tight_layout()\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.figure(figsize=(8,4))\n",
|
|||
|
|
"plt.plot(C, N, color=\"blue\", label=\"Naive\")\n",
|
|||
|
|
"plt.plot(C, O, color=\"red\", label=\"Optimised\")\n",
|
|||
|
|
"plt.xlabel(\"$\\chi$\")\n",
|
|||
|
|
"plt.ylabel(\"FLOPs\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.legend()\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.tight_layout()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 179,
|
|||
|
|
"id": "688ca9b4",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"c=2\n",
|
|||
|
|
"np.random.seed(42)\n",
|
|||
|
|
"lambda1 = np.random.normal(size=(c, c))\n",
|
|||
|
|
"lambda2 = np.random.normal(size=(c, c))\n",
|
|||
|
|
"lambda3 = np.random.normal(size=(c, c))\n",
|
|||
|
|
"G1 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
"G2 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
"U = np.random.normal(size=(c, c, c, c))"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 180,
|
|||
|
|
"id": "6eeac349",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"text/plain": [
|
|||
|
|
"array([[[ 0.24196227, -1.91328024],\n",
|
|||
|
|
" [-1.72491783, -0.56228753]],\n",
|
|||
|
|
"\n",
|
|||
|
|
" [[-1.01283112, 0.31424733],\n",
|
|||
|
|
" [-0.90802408, -1.4123037 ]]])"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"execution_count": 180,
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "execute_result"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"lambda1\n",
|
|||
|
|
"G1"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 319,
|
|||
|
|
"id": "9e977aa4",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def mY_Z(lambda1, lambda2, lambda3, G1, G2, U):\n",
|
|||
|
|
" c=2\n",
|
|||
|
|
" A = np.tensordot(lambda1,G1, axes=([1], [1]))\n",
|
|||
|
|
" B = np.tensordot(lambda2, G2, axes=([1],[1]))\n",
|
|||
|
|
" C = np.tensordot(B, lambda3, axes=([2], [0]))\n",
|
|||
|
|
" D = np.tensordot(A, C, axes=([2], [0]))\n",
|
|||
|
|
" Z = np.tensordot(D, U, axes=([1, 2], [2,3]))\n",
|
|||
|
|
" \n",
|
|||
|
|
" return Z"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 365,
|
|||
|
|
"id": "324465ef",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"assert np.allclose(mY_Z(lambda1, lambda2, lambda3, G1, G2, U), Z(lambda1, lambda2, lambda3, G1, G2, U))"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 366,
|
|||
|
|
"id": "207b4224",
|
|||
|
|
"metadata": {
|
|||
|
|
"scrolled": false
|
|||
|
|
},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"O, my, Oerr, myerr = [], [], [], []\n",
|
|||
|
|
"C = np.arange(3, 50)\n",
|
|||
|
|
"\n",
|
|||
|
|
"for c in C:\n",
|
|||
|
|
" np.random.seed(42)\n",
|
|||
|
|
" lambda1 = np.random.normal(size=(c, c))\n",
|
|||
|
|
" lambda2 = np.random.normal(size=(c, c))\n",
|
|||
|
|
" lambda3 = np.random.normal(size=(c, c))\n",
|
|||
|
|
" G1 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
" G2 = np.random.normal(size=(c, c, c))\n",
|
|||
|
|
" U = np.random.normal(size=(c, c, c, c))\n",
|
|||
|
|
" \n",
|
|||
|
|
"\n",
|
|||
|
|
" time1 = %timeit -o -q -n 2 -r 2 mY_Z(lambda1, lambda2, lambda3, G1, G2, U)\n",
|
|||
|
|
" time2 = %timeit -o -q -n 2 -r 2 Z(lambda1, lambda2, lambda3, G1, G2, U)\n",
|
|||
|
|
" \n",
|
|||
|
|
" my.append(time1.average)\n",
|
|||
|
|
" myerr.append(3*time1.stdev)\n",
|
|||
|
|
" O.append(3*time2.average)\n",
|
|||
|
|
" Oerr.append(time2.stdev)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 367,
|
|||
|
|
"id": "f7cb2b04",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAEQCAYAAACutU7EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA2DklEQVR4nO3dd5hU5dnH8e9NB1FAQMAYxK5gwbgWrFhIwViixhYLMUosQYkliYmvLWoSotgJYGJAjIlYghogooKKgmUBRbArYBRpUkRZXBbu94/njAzjttmdmTPl97muc83OOc+cuWcv9uGep5q7IyIiIlJImsQdgIiIiEi6lMCIiIhIwVECIyIiIgVHCYyIiIgUHCUwIiIiUnCUwIiIiEjBUQIjIiIFz8zmm9mofL9nPd6zr5m5mfXN5fsWIiUwkjVmNiD6Q3QzO6SGMu9H15/NcXgikiVm1t3MhkcJwFdmtsTM/m1mBzbyvv3N7NoMhSkFrlncAUhJWAucDkxNPmlmBwA7RNdFpAiY2UHAhOjpPcBbQDdgAPCCmV3i7nc28Pb9gYuAa6u5tguwoYH3rUk27ikZogRGcmEC8GMzu9jd1yWdPx14G1gfT1gikklm1gF4GKgADnL3D5Ku3QI8CdxmZjPcfVom39vdv8rk/bJ1T8kcdSFJLvwT2BL4XuKEmTUFTgEeSDrXxMw+MrPHU29gZs3MbLGZPZiLgEWkQX4OdAWuSE5eANy9AjgbcODqxPmkMR8/MbPrzGyhma0xsyfNbKekcqMIrS8kdU27mfWIzm0yXiXpvqeb2TVm9omZfWFmj5pZBzNrYWZDo3rlCzMbbWatk2OubgyMmV1oZm+Y2ZdmtsrMXjOzn6eU6WZmfzWzRVEX2ltmdkHqL8vMtjGzcdG9lpjZrUDLNH7fJU0tMJILHxO6j04H/hOdOwrYipDAnALg7hvM7H7gcjPr6O6fJd3ju1H5+3IWtYik6xhCl/DY6i66+zwzewE4wsxaR0lNwq+BpsDNQAfgEmCKme3p7suBEcDWQD/gzKTXLa0jpl9FMf0R2BEYRGj1XQN0Aa4D9gfOAhaQlFylMrOfAXcTWpnuApoDvYADo/gws62Al6LPMgxYAhwJDIvqtRuicq2BZ4DuwB3AQuAnwBF1fB6JKIGRXHkAGGpmm7n7l4Q/1Jfd/QMzSy53H3AlIakZlnT+DEJF9WSO4hWR9PUE3qmj6+V14DBCMvFG0vkuwC7uvhLAzKYQ/oO/FLjK3aeb2btAP3e/P42YWgD7uXtldN/OwKnAU8D3PexoPCxq7TmHWhIY4IfAXHf/cS1lbiC0ouzh7onkariZ3QP81szuij7jQGBn4BR3HxvFNhKYlcZnK2nqQpJceYjwbeX46JvH8cA/Ugu5+9vAKyR9wzKztsBxwD/dvSon0YpIQ2wOrK6jTOL6Finn70skLwDuPhmYS0gaGmNMInmJvAwY8PcoeUk+v7WZ1daFswrYxsz2re6ihW9jJwHjATezTokDmAS0JrT2QBiQvJjQmgN83c3217Q+XQlTAiM5ETUBP0loeTkWaAPUNJ5lNHCAme0QPT8hKj8m23GKSKOsJiQxtUlcT0103qum7LtAj0bG9FHK81XR4/+qOW9A+1ru9SdC3K+Y2QfRVPHDk653JnR/nUNoMU4+Et1qW0WP2wIfuHvqLKd3a/008jV1IUkuPUDoItoCeNrdl9RQ7l/ArYRuo+uix7fdvTwnUYpIQ70F7G1mLWvpRtoTWEf1CUs21DTLsabzVsN53P0tM9uF0HryPULr0M/NbJi7X8TGRoF/AvfWcJu5dYcs9aEERnLpMeAr4CDCbIRquftyM/sPcEbUJ3wEtfdLi0h++A/QB/gx8I1xKtGMoUMIX2AqUi7vlFqeMEZkftJzr6ZMTrn7GkK3z8Nm1gwYBVxoZjcBiwgtNM3c/ek6brUA2MvMmqS0wuychbCLkrqQJGeiP/wLCK0q/66j+GjCIL/bCP9O0xm0JyLxGEGYdTPEzLZLvhCNffs7oYXj+mpee5aZtU8qfwRhhs/4pDJfRtc6ZDbs+jGzjsnPozF5iYHI7d19PSG5Od7M9qrm9Z2Tnk4gDFw+Kel6a+DcTMddrNQCIznl7vUdxzKR0G98MvCsu6f2Y4tInnH3z8wsMYh1VjTz5i3C2jADCF9KLqlhEbvFwItm9jfCOJTBwKfA0KQyiW7ku8xsIlAFPBHNbMyFSWa2BHiB0NqSmJY9m/A5AX4D9AWmR59/LmFcTG/gR0CrqNw9wC+A0Wa2D/AJobtci+fVkxIYyUvuvs7M/glcjAbvihQMd59qZnsCvyV8AelGGCD7IvBTd3+xhpf+ibB0/xWEBGYqMChlPahHCa2yp0WHAdsRtczkwAjCelaDCWP5FhLGutyQ6AZy9yVmtj/wf4TZlhcAywkJzmWJG7n7GjM7EriTkMisIczMnAj8Nzcfp7DZprPIRPKHmf2Z8Ifdxd0/jzseEck8C7suTwFOc/d/xRuNFBKNgZG8FK3FcCbwbyUvIiKSSl1IkleiZbiPIvQVb0WYTi0iIrIJJTCSb3oS+oGXAr9091djjkdERPKQxsCIiIhIwdEYGBERESk4RduF1KlTJ+/Ro0fcYYgUlRkzZixz9851lyw9qnNEMq+2OqdoE5gePXpQXq6tc0QyycwWxB1DgpkdClwO7ANsTVhjZFQt5XsA86q59AN3/29SucMIi6f1IqzzMcTdh9cVj+ockcyrrc5RF5KIFKq2wBzgEiB1X53afJ+wuFrimJy4EC1/PwGYBuwN/AG408xOzFDMIpIhRdsCIyLFzd0nEJINzGxUGi/9zN0X1XDtfGChuw+Knr8Vrap6OfBIQ2MVkcxTC4yIlJpHzWyJmb0Y7duTrA8wKeXck0CZmTXPTXgiUh9KYESkVHxBaEk5GegPPAM8aGZnJJXpSthUMNliQmt1p1wEKSL1oy4kESkJ7r4MuCXpVLmZdQJ+BdzfkHua2UBgIED37t0bHaOI1J9aYESklL0M7JT0fBHQJaVMF6AKWJb6Yncf6e5l7l7WubNml4vkkhIYkRL32GNwzDGw7Bv/PZeE3sCnSc+nA/1SyvQDyt19Xa6CEilqY8bA0UdDVVWjbqMuJJES99xz8PTT0L593JGkx8zaAjtGT5sA3c2sN7Dc3T8ysz8A+7n7kVH5s4F1wCxgA3AMcBHw66TbDgd+YWa3ASOAg4ABwGnZ/jwiJWPOHHjmGWjWuBRECYxIiZs5E/baq9F1SRzKgClJz6+LjtGEpKMbsEPKa64CtgXWA+8C57j71+Nf3H2emfUn7IJ+AWEhu4vdXVOoRTJl5cqMfGMqvCpLRDJmwwaYNQt+8pO4I0mfuz8LWC3XB6Q8H01Ibuq673PAdxoZnojUJEMJjMbAiJSwDz+Ezz+HffaJOxIRKRlKYESksWbMCI/fUXuDiOSKEhgRaayZM6FFC+jVK+5IRKRkKIERkcaaORP22CMkMSIiOaEERkQawz10Ian7SERyxl0JjIg0zoIFsGKFEhgRyaG1a6GyEjp0aPStlMCIlKiZM8OjEhgRyZmVK8OjWmBEpKFmzICmTWHPPeOORERKhhIYEWmsmTPD7KNWreKORERKhhIYEWkMDeAVkVgogRGRxli4EJYu1Qq8IpJjSmBEpDG0Aq+IxEIJjIg0xsyZYBZ2oRYRyZlEAtOuXaNvpQRGpATNnAm77gqbbRZ3JCJSUlauDDMHMjB7QAmMSAmaMUPjX0QkBitWZKT7CJTAiJScRYvCIF6NfxGRnMvQNgKgBEak5MyaFR6VwIhIzimBEZGGSmwhsPfe8cYhIiVICYyINFR5Oey0E2yxRdyRiEjJUQIjIg3hDtOnwwEHxB2JiJQkJTAi0hDz58PixdCnT9yRNJ6ZHWpmj5vZJ2bmZjagjvJ9zewxM/vUzNaY2WwzO6eaMl7NsWtWP4xIKXDPaALTLCN3EZGCMH16eCyGBAZoC8wB7ouOuhwIvAEMAT4FvgeMNLO17v5AStlewPKk50sbH65IiauogHXrlMCISPqmTQuL1+2+e9yRNJ67TwAmAJjZqHqUvynl1F/
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 576x288 with 2 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAEQCAYAAACutU7EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3wElEQVR4nO3dd5hU5dnH8e+9DXbpUpYmgkAEK5qNXVxFJGLFEo0VxYLY0KiYaDQaX31jwd4VBdEEg8bYG0pQUXgXUQQioosoHaQKS9nd5/3jmWGHYXbZgdk5M7O/z3Wd68ypc88cYG6eas45RERERNJJVtABiIiIiMRLCYyIiIikHSUwIiIiknaUwIiIiEjaUQIjIiIiaUcJjIiIiKSdnKADCFKrVq1c586dgw5DREREYpgyZcoy51zrWMfqdQLTuXNnSkpKgg5DREREYjCzudUdUxWSiIiIpB0lMCIiIpJ2lMCIiIhI2klqAmNmvc3sNTObb2bOzAZu4/zOofOil99GnXe4mU0xs/VmVmpmg+v0g4iIiEigkl0C0xiYDlwFlMVx3W+BdhHLh+EDZtYFeAuYCOwL3Ak8ZGanJChmERERSTFJ7YXknHsLn2xgZs/FcenPzrlF1RwbDCxwzl0R2v6vmR0AXAu8vL2xAqxevZolS5awadOmHbmNpKjc3FzatGlD06ZNgw5FRETilC7dqF8xs4bAbOA+59zYiGMHAe9Fnf8ucJ6Z5Trntiv7WL16NYsXL6ZDhw7k5+djZtsXuaQk5xxlZWXMnz8fQEmMiEiaSfVGvL/gS1J+B/QHxgFjzOzsiHPaAoujrluMT85aRd/QzC42sxIzK1m6dGm1b7xkyRI6dOhAQUGBkpcMZGYUFBTQoUMHlixZEnQ4IiISp5QugXHOLQPujdhVYmatgOuB0dt5zyeBJwGKiopcdedt2rSJ/Pz87XkLSSP5+fmqIhQR2RHFxX49fnxS3zbVS2BimQR0j9heBBRGnVMIlAPLduSNVPKS+fSMRUR2kKu2LKBOpWMC0wtYGLH9GdA36py+QMn2tn9JtOLiqgRVREQkY2zYAJ98AvPmJf2tk1qFZGaNgW6hzSygk5n1ApY75340szuB/Z1zfULnnwdsAqYClcDxwGXAsIjbPg5cbmb3A08AhwADgd/X9ecRERGp1+bOhcpKyEl+i5Rkl8AU4ZORqUA+cGvo9W2h4+2ArlHX3ASUAP8HnAFc4Jy7L3zQOTcH38C3N/AlcCNwpXNuh7pQp6uBAwdiZgwaNGirY8OGDcPMOO644wKITEREMk5pqV8H0GY02ePAjAeqbXTgnBsYtT0SGFmL+/4H2G8Hw8sYO++8My+99BIPPvggjRo1AqC8vJxRo0bRqVOngKMTEZGMEU5gGjZM+lunYxsY2Ya9996b7t2789JLL23e9+abb9KwYUOKQ41xJkyYQG5uLosWbTk+4I033sjee++dzHBFRCRdlZb65OXTT5P+1kpgMtSgQYMYMWLE5u0RI0Zw/vnnb+5107t3b7p27cqoUaM2n1NZWcmoUaNiVj+JiIhspbQUunSBrOSnEyk9DkyqGToUvvwy/uvC12xPT6ReveD+++O/7swzz+Taa69l9uzZNGnShHfeeYeHHnqIm2++efM5F154Ic888wzXX389AO+++y5Llizh7LPPru62IiIiVUpLYdddA3lrlcBkqBYtWjBgwABGjBjByJEjKS4u3qr9y3nnnUdpaSkTJ04EfCnNSSedRMuWLYMIWURE0olzgSYwKoGJw/aUhEBggxRywQUXcN5559G4cWNuu+22rY63bt2aE044gREjRrDbbrvx2muv8frrryc3SBERSU8//wxr1iiBkcTr06cPeXl5LFu2jJNOOinmORdddBGnnnoqu+66K23btuWoo45KbpAiIpKewj2QlMBIopkZ06ZNwzlHgwYNYp7Tt29fWrZsya233soNN9xAVgANsUREJA0FnMDo1yrDNWnShKZNm1Z73Mw4//zz2bRpE+eff34SIxMRkbQWTmC6dAnk7VUCk2Gee+65uI8vXLiQPn360Llz5zqJSUREMlBpKRQWQmjA1GRTAlOPrVq1ipkzZzJq1KgtBr0TERHZpgB7IIESmKRIdu+j2jrxxBOZPHkygwYN4thjjw06HBERSSelpXDooYG9vRKYemx8qmZWIiKS2jZuhJ9+CrQERo14RUREJD4//giVlUpgREREJI0E3IUalMCIiIhIvJTAiIiISNopLYUGDaB9+8BCUAIjIiIi8Skt9QPYBTh6uxKYZCgurprRUUREJN0FPAYMKIGR7TB+/HjMjGXLlu3QfQYOHMhxxx2XoKhiO+644xg4cGCdvoeISL3iHHz/vRIYqRvz58/n4osvpmPHjuTl5dGhQwcuuugi5s2bF9d9OnfuzD333LPFvoMPPpiFCxfSsmXLHYrxgQceYPTo0Tt0DxERSbIVK2D1aiUwknhz5syhqKiI6dOnM3LkSL777jtGjx7NjBkz+M1vfsMPP/ywQ/fPy8ujbdu2mNkO3adZs2Y0b958h+4hIiJJlgI9kEAJTEa67LLLyMrK4oMPPqBPnz506tSJI444gg8++ICsrCwuu+wyAIqLixk8eDBXXXUVLVq0oEWLFlx33XVUVlZuPj537lyuu+46zGxzwhJdhfTcc8/RuHFj3n77bXr06EFBQQEnnHACq1atYuzYsXTv3p1mzZpxzjnnUFZWtjnO6CqkCRMmcOCBB9K4cWOaNWvG/vvvz/Tp0zcfnzhxIocffjgFBQV06NCBSy+9lNWrV28+vm7dOgYOHEjjxo0pLCzkjjvuqLsvWUSkvqqPCYyZ9Taz18xsvpk5Mxu4jfOLzezfZrbQzNaZ2TQzuyDGOS7G0qNOP0yKWr58Oe+88w6XXXYZBQUFWxwrKChgyJAhvP3226xYsQKAF154gcrKSj777DOeeOIJnnzySe6//34AXnnlFTp27MjNN9/MwoULWbhwYbXvu2HDBu69915eeOEFxo0bR0lJCaeccgojR47k5Zdf5tVXX+WNN97g0UcfjXl9eXk5J554IoceeihfffUVkyZNYujQoWRnZwPw9ddfc/TRR3PCCSfw1Vdf8corr/Dll19ywQVVfxyuvfZa3n//fV5++WXGjRvH1KlTmTBhwo58nSIiEi2cwHTpEmgYyZ4LqTEwHRgVWrblYOBr4C5gIdAPeNLM1jvnXow6dw9gecT20h0PN8rQofDll/FfF75me3oi9eoFoYSiNmbPno1zjp49e8Y8vvvuu+OcY/bs2QC0a9eOBx98EDOjR48efPvttwwfPpxrrrmGnXbaiezsbJo0aULbtm1rfN/y8nIeeeQRdtttNwDOPPNM7rvvPhYvXkyrVq0AP3nkRx99xB/+8Ietrl+9ejUrV67k+OOPp2vXrgD06FGVg959992cfvrpW1z72GOPse+++7JkyRIKCgp45plnGDFiBP369QPg2WefpWPHjrX96kREpDZKS6FNG2jcONAwkprAOOfeAt4CMLPnanF+dB3AY2Z2BHAKEJ3ALHHO7Vi3mHrowAMP3KIty0EHHcSf//xnVq9eTdOmTWt9nwYNGmxOXgAKCwtp27bt5uQlvG/mzJkxr99pp50YOHAg/fr1o0+fPvTp04dTTz2VTp06ATBlyhS+++47xowZs/ka5xwA33//PQUFBWzcuJGDDjpo8/HGjRuz11571foziIhILaRADyRIz9momwKxutKUmFkDYCZwu3Puo4S/cxwlIVsIl7wkYfbnbt26YWbMnDmTAQMGbHV85syZmBndunVL6Pvm5Gz5R8nMyM3N3WpfuH1NLM8++yxDhw7lnXfe4bXXXuPGG2/k1VdfpV+/flRWVnLhhRdy9dVXb3Vdhw4d+PbbbxPzQUREpGalpXDwwUFHkV6NeM3sOKAP8GTE7oXApfhSmZOBWcA4MzusmntcbGYlZlaydGnia5mC1rJlS/r168ejjz7KunXrtji2bt06HnnkEY455hh22mknACZNmrS5JAPg888/p3379ptLX/Ly8qioqEha/Pvssw/Dhg1j/PjxFBcXM3LkSAD
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 576x288 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"plt.rcParams.update({'font.size': 14})\n",
|
|||
|
|
"plt.figure(figsize=(8,4))\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.subplot(121)\n",
|
|||
|
|
"plt.plot(C, my, color=\"blue\")\n",
|
|||
|
|
"plt.title(\"My\")\n",
|
|||
|
|
"plt.xlabel(\"$\\chi$\")\n",
|
|||
|
|
"plt.ylabel(\"time,sec\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.subplot(122)\n",
|
|||
|
|
"plt.plot(C, O, color=\"red\")\n",
|
|||
|
|
"plt.title(\"Optimised\")\n",
|
|||
|
|
"plt.xlabel(\"$\\chi$\")\n",
|
|||
|
|
"plt.ylabel(\"time,sec\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.tight_layout()\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.figure(figsize=(8,4))\n",
|
|||
|
|
"plt.errorbar(C, my, myerr, color=\"blue\", label=\"My\")\n",
|
|||
|
|
"plt.errorbar(C, O, Oerr, color=\"red\", label=\"Optimised\")\n",
|
|||
|
|
"plt.xlabel(\"$\\chi$\")\n",
|
|||
|
|
"plt.ylabel(\"time,sec\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.legend()\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.tight_layout()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "1b180a50",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"### Задача 7"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 64,
|
|||
|
|
"id": "1a00a71e",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"N = 1000\n",
|
|||
|
|
"K = 3\n",
|
|||
|
|
"d = 3\n",
|
|||
|
|
"L = 10\n",
|
|||
|
|
"\n",
|
|||
|
|
"np.random.seed(42)\n",
|
|||
|
|
"mu_true = np.random.uniform(-L, L, size = (K, d))\n",
|
|||
|
|
"data = np.random.normal(mu_true, size = (N, K, d))\n",
|
|||
|
|
"data = np.vstack(data)\n",
|
|||
|
|
"np.random.shuffle(data)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 65,
|
|||
|
|
"id": "381a2bc2",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def scatter_plot(data, col=None):\n",
|
|||
|
|
" from mpl_toolkits.mplot3d import Axes3D \n",
|
|||
|
|
" %matplotlib inline\n",
|
|||
|
|
" fig = plt.figure()\n",
|
|||
|
|
" ax = fig.add_subplot(111, projection='3d')\n",
|
|||
|
|
" ax.scatter(data[:,0], data[:,1], data[:,2], s = 0.5, color=col)\n",
|
|||
|
|
" plt.show()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 66,
|
|||
|
|
"id": "1cc70b18",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAADwCAYAAAAzS5nVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABmq0lEQVR4nO29eXwbd53//xydlmXL933Edpz7sHMnpbSlJz2TcrUUKNBCYVm6/NhdCizLUhZYoMC3sCwUlpulpfSkpS2FFiilLW2OJj5jx47t+LYsWbbua2Z+fzgzlRXZliz5SKLX45FHEtkz85E0r3m/P+/j9RZkWSaNNNI4t6BZ7gWkkUYaqUea2GmkcQ4iTew00jgHkSZ2Gmmcg0gTO400zkHo5vl5OmSeRhqLDyHVJ0xb7DTSOAeRJnYaaZyDSBM7jTTOQaSJnUYa5yDSxE4jjXMQaWKnkcY5iDSx00jjHESa2GmkcQ4iTew00jgHkSZ2Gmmcg0gTO400zkGkiZ1GGucg0sROI41zEGlip5HGOYg0sdNI4xxEmtjLAFmWCYVCiKJIWiU2jcXAfEILaaQYkiQRCoXw+XwACIKATqdT/2g0GgQh5X33aZxnEOaxGGlzkiLIsowoioRCIQDC4bD6uiRJ6u8pRNfr9Wi12jTRzw+k/AtOE3sJEOl6KyRVCB7rdyVJQpZlBEFAEAT0er1q0ZXX0jinkCb22QZJkggGgzOIKssywWAwLoJGEh1Ao9Gg1+tVi54m+jmBNLHPFsiyTDgcJhwOn0G+RIgdfU5ghkXXaDSIokhWVpbquqdx1iEtZng2QCFuLFInA+VcWq1WDbQBvP7663g8HpxOJ06nE6/XSygUmrF3T+P8QjoqnmIoVjrS9Z4Nyu8sFMqxCtkVix4MBgkGg8C06x4djEvj3Eea2CmCLMu4XC5cLhcFBQVzEihyr51KKETXarXqmhTvIZLoSjAuTfRzF2lipwBKbtrlcmG1WiksLFzuJQHE3NvLskwgECAQCAAzia5E3dM4+5EmdhKIzk0rlnKlIh6ia7XaGRY9TfSzE2liLxDRuWklQn02lYjGIrokSfj9fvW1NNHPTqSJvQDEyk0Di7JvXkqkiX7uIE3sBBDpeisWOhKCIJxTKabZiO7z+RgaGqKwsBCz2Zwm+gpEmthxQokuS5I0axrrbLfY8yFyy+F2uykoKFCJHhmRTxN9+ZEmdhyYzfWOxrlO7GgoJFc8l0iLrnxGkZ1raaIvHdLEngPRZaHz5XzPJ2LHKq6JtOjK74iiqHayAWqxTLpFdXGRJvYsUHLTc7ne0TjbouLJIJ6quVh79EiiKxVzaaKnHmliR0HJ605NTWGxWBKq9T7fLHaiiIfoadGJ1CBN7Agorrfb7aavr4/GxsaEjo83Ku5yuWhtbUWv15Ofn09eXh6ZmZln3U2c7HpjET0cDqsFP2miLxxpYp9GpOsd2VCRCOaz2LIsMzQ0xMDAABs2bECWZRwOBz09PXi9XrKyslSiZ2RkJPN2Fh3JNrDEQiyih0KhGUSXZZmMjAz0en2a6HPgvCd2rNy0RqNZUD56LmKHw2Ha29sRBIHdu3cjSRKSJGE2m6msrESWZdxuNxMTE3R0dBAMBrFYLOTn55Obm4vBYEj2raYUi0HsaCh78MhrtrW1UVNTo3o40Q0taaJP47wm9my56YUSe7bgmcvloqWlhVWrVlFRUQFwxvkFQSA7O5vs7GxWrVqFJEk4nU4mJiYYHBxEkiRycnLIz88nJycHnW55v7qlIHY0lAen0oKqfH/RDS1pdZnzmNhz5aYXGt2OZbEHBwfp7+9ny5YtZGdnx30ujUZDbm4uubm5wLTFn5qaYmJigt7eXvXn+fn5yxKwW64gYeQDJZZFV4iueF+RvejnE9HPO2LHk5tOpjRUueFFUaS9vR1Zltm9e3fSFlan01FQUEBBQQEwLaYwOTnJ6OgoXq+XY8eOqfvzrKysJbmBl4MkkiTNWk8QSfTzXXTivCJ2PGWhsHBXXDmf2+2mpaWFqqoqKioqFoUABoOB4uJiiouLcblcrF+/HofDQX9/P263m8zMTPLy8sjPz8dkMqV8DcvhisPcxI5ELNEJOH+Ift4QOzqNMtdNmcwNGwwGaW5uTtj1ThYZGRmUlZVRVlaGLMt4vV4cDgfd3d34fD6ys7NVohuNxqSvt1yuuPJQThTnm7rMOU/sudRCUwlRFDl+/DiiKKbE9U4GgiBgNptnRNxdLhcTExO0t7cTDodnRNz1ev2Cr7PUSJWnkKi6zNlG9HOa2AspC10IFNe7srKSqampZY9YR0MQBCwWCxaLhZqaGkRRVCPu/f39yLJMbm4ueXl55ObmxqUEs1yu+GJ9j/EQ3el0UlxcrKrEruRA3Mq6A1MEJTfd3d2t3rCLheHhYfr6+lTXe2BgYNGulSpotVry8vLUzyUcDjM5OcnExAQ9PT3qz/Pz88nOzp4zWHWuIhbRu7q6sFgs6msruUX1nCN2pGRRKBSa0VmUSiiudzgcXpDrvdAA3WJAp9NRWFioijAGg0EcDgcDg0O8PjDFmkITFcX5WHLzsAc01BRkLvOKlx4K0SP36CtZXebs2TTEAUmSCAQCqg7ZYpHH4/Fw8OBBLBYLDQ0NC3K9V3KziMFgoKSkhKKq1bT6ctHmVaDT6TjY1sN3f3+UFw41EwwG8Xq9K/p9pBLR71O5v7Rarbr/Vojudru588476e7uXqbVniMWe7bctFarTTmxR0ZG6O3tZfPmzTPcsuXGYux5SyxGPntVPbmm6f2kJiufqtoQlWZobm6mu7sbv9+vRtzz8vJSEnFfiRBFcc7YQ7Tr3tvbu6yxlrOe2HPlpgVBQBTFlFxHFEU6OjoIhULLHvVeLHRbPWg0UFdoVl/Ly3wjYv5/rw1RnW9iw45y9Ho9W7duRZIkXC4XJwbH+cbzvVxZraG2NE8l+rnyOSnNQfFCaepZLpzVn7qyj55NsihZi62c1+v10tzcTHl5OVVVVSmzjMsVWY4FWZb5z2e7MGgEfvzehjN+Lsky++ryWFOUSeeYm4A47ZraPSEMhkyKSssZDbl42WkiqzgTrdNJf38/oiRzxGHgTWuK2FRdnJT2+nK6/aIoJpTu8ng8aWIninglizQazaxzqOeDUvc9NjZGT08PmzZtIicnJ5llLwiHTk1SmZtBWc7itnEKgsBdl9eh1Qj02LxoNQKr8k3qz0edAf7QPk5wbQEPHR0hXw7z5n3wo5f7Meg0fPLSOm6/oIrX+hzk5uSwumy6OMfl9fP0c130DtvwWgfQ6nQUFUyXvs4VcY+F5XwQzueKRyMYDC7rtuSsI3YiuelkgmeCINDe3k4oFGLXrl0JFXGk6gYUJZnmISeyLC86sQFyTHqaBp282GVHqxX4rxvW8/v2cbQagRNjbiY8AZ7vHOeydYUMnXICsHNVDi92O3AHwoREGV9IYm3xtCvvD4n83+FRbt5Tw4bSbB45OoIvEKQqw8jw8DAul4uMjAzVbTebzXN+bstJ7ERd8eXGWUPs6HE68TzpF0psr9eL2+2mqKiITZs2JXQzKZY+FTegViPwgb1VaDXzn2uh15RlmXF3kKIsA2POAIdOTeELibyldrrZxO4Ocrh/Ck8gTGG2AYcnxKs9E+zMnr7Ji7MNFJn13PNcD55giNv2VmH3BGkbdvGm+gJqCzLJNU0/FLdVWgiJEqWl2RSXlNA/4aPAKOOcmqSvr091XxWim0ymGWuNt058MZCIK74StlhnBbFjjdOJBwsh9ujoKCdPnlTLMRP9ghLVPeuxeXAFRBoqYkfY5yO11RUguACnRJRkXh+YIhCW+EO7lVt2VdBYlUOPzc19f7Nz/dZSBEGgoSKL/335FKvyTDi9IYqyjYxM+fndeIjCjnEeOjJMMCwx6vRhc4exOQNsKM/mtb4pcjP1XL2pWL1mYZaBbz3fw9ZKC1vKs/nFq4O8f28l6yoqqKioUMUmHA4HJ06cIBAIYLFYVKLLsrysxE7EYi83uVc8sePV9I4FrVYbd1RckiQ6Ozvx+Xzs3r2blpaWlKuoRGLM6WfC7adr3EfzsAtBlui
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 432x288 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"if d == 3:\n",
|
|||
|
|
" scatter_plot(data, None)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 67,
|
|||
|
|
"id": "e32355c0",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def naive_dist_i(x, mu):\n",
|
|||
|
|
" # x: N datapoints, mu: N cluster centers\n",
|
|||
|
|
" # returns: D_{i}, squared distances from x[i] to mu[i]\n",
|
|||
|
|
" dist = np.zeros(x.shape[0])\n",
|
|||
|
|
" for i in range(x.shape[0]):\n",
|
|||
|
|
" dist[i] = np.sum((x[i] - mu[i])**2)\n",
|
|||
|
|
" return dist\n",
|
|||
|
|
"\n",
|
|||
|
|
"def naive_dist_ij(x, mu):\n",
|
|||
|
|
" # x: N datapoints, mu: K cluster centers\n",
|
|||
|
|
" # returns: D_{ij}, squared distances from x[i] to mu[j]\n",
|
|||
|
|
" dist = np.zeros((x.shape[0], mu.shape[0]))\n",
|
|||
|
|
" for i in range(x.shape[0]):\n",
|
|||
|
|
" for j in range(mu.shape[0]):\n",
|
|||
|
|
" dist[i, j] += np.sum((x[i] - mu[j])**2)\n",
|
|||
|
|
" return dist\n",
|
|||
|
|
"\n",
|
|||
|
|
"def naive_kmeans():\n",
|
|||
|
|
" ss_list=[]\n",
|
|||
|
|
" mu = data[np.random.choice(range(data.shape[0]), K, replace=False)]\n",
|
|||
|
|
" c = np.random.randint(low=0, high=K-1, size=data.shape[0])\n",
|
|||
|
|
" for n in range(10):\n",
|
|||
|
|
" c = np.argmin(naive_dist_ij(data, mu), axis = 1) \n",
|
|||
|
|
" ss = np.mean(naive_dist_i(data, mu[c]))\n",
|
|||
|
|
" ss_list.append(ss) \n",
|
|||
|
|
" for i in range(K):\n",
|
|||
|
|
" cluster_members = data[c == i]\n",
|
|||
|
|
" cluster_members = cluster_members.mean(axis = 0)\n",
|
|||
|
|
" mu[i] = cluster_members\n",
|
|||
|
|
" return ss_list, c"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 79,
|
|||
|
|
"id": "6b0e811d",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def dist_i(x, mu):\n",
|
|||
|
|
" # x: N datapoints, mu: N cluster centers\n",
|
|||
|
|
" # returns: D_{i}, squared distances from x[i] to mu[i]\n",
|
|||
|
|
" return np.sum((x - mu)**2)\n",
|
|||
|
|
"\n",
|
|||
|
|
"def dist_ij(x, mu):\n",
|
|||
|
|
" # x: N datapoints, mu: K cluster centers\n",
|
|||
|
|
" # returns: D_{ij}, squared distances from x[i] to mu[j]\n",
|
|||
|
|
" n, k = x.shape\n",
|
|||
|
|
" x_mod = np.transpose(x[None,...].repeat(k, axis=0), axes=(1, 0, 2))\n",
|
|||
|
|
" mu_mod = mu[None, ...]\n",
|
|||
|
|
" return np.sum((x_mod-mu_mod)*(x_mod-mu_mod), axis=2)\n",
|
|||
|
|
" \n",
|
|||
|
|
"def kmeans():\n",
|
|||
|
|
" ss_list=[]\n",
|
|||
|
|
" mu = data[np.random.choice(range(data.shape[0]), K, replace=False)]\n",
|
|||
|
|
" c = np.random.randint(low=0, high=K-1, size=data.shape[0])\n",
|
|||
|
|
" \n",
|
|||
|
|
" for _ in range(10):\n",
|
|||
|
|
" c = np.argmin(dist_ij(data, mu), axis = 1) \n",
|
|||
|
|
" ss = np.mean(dist_i(data, mu[c]))\n",
|
|||
|
|
" ss_list.append(ss) \n",
|
|||
|
|
" for i in range(K):\n",
|
|||
|
|
" mu[i] = data[c==i].mean(axis=0)\n",
|
|||
|
|
" return ss_list, c"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 82,
|
|||
|
|
"id": "84416419",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"885 ms ± 51.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"name": "stderr",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"/var/folders/h9/ll8t8frd1r575llt96l1glhm0000gn/T/ipykernel_21445/3230926319.py:24: RuntimeWarning: Mean of empty slice.\n",
|
|||
|
|
" mu[i] = data[c==i].mean(axis=0)\n",
|
|||
|
|
"/Users/goloshch/.conda/envs/NPM/lib/python3.10/site-packages/numpy/core/_methods.py:181: RuntimeWarning: invalid value encountered in true_divide\n",
|
|||
|
|
" ret = um.true_divide(\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"7.57 ms ± 419 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"%timeit naive_kmeans()\n",
|
|||
|
|
"%timeit kmeans()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 83,
|
|||
|
|
"id": "84fb1f46",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPYAAADwCAYAAAAzS5nVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABjO0lEQVR4nO29eXgcZ5nu/Xurqru1tvZ9seV93+0sBBJCFpKY2GENOSTDECAznMnkY2ACDMMAhwFCgBMYBhiGgYE5wIQkdggkmUAChC2QxMa2ZMu2LFuyFmvf1VIvVfV+f7Sq3Gq3pF4l2e77unxZW3W9XV13Pc/7LPcjpJSkkUYalxaUhV5AGmmkkXykiZ1GGpcg0sROI41LEGlip5HGJYg0sdNI4xKENsfv0yHzNNJIPUSyXzBtsdNI4xJEmthppHEJIk3sNNK4BJEmdhppXIJIEzuNNC5BpImdRhqXINLETiONSxBpYqeRxiWINLHTSOMSRJrYaaRxCSJN7DTSuASRJnYaaVyCSBM7jTQuQaSJnUYalyDSxE4jjUsQaWIvAKSUBAIBDMMgrRKbRiowl9BCGkmGaZoEAgEmJycBEEKgaZr9T1EUhEh6330alxnEHBYjbU6SBCklhmEQCAQA0HXd/rlpmvbfWUR3OByoqpom+uWBpH/AaWLPA0Jdb4ukFsEj/a1pmkgpEUIghMDhcNgW3fpZGpcU0sS+2GCaJn6/fxpRpZT4/f6oCBpKdABFUXA4HLZFTxP9kkCa2BcLpJTouo6u6xeQLxZih78mMM2iK4qCYRjk5OTYrnsaFx3SYoYXAyziRiJ1IrBeS1VVO9AG8Oc//xmPx8Po6Cijo6NMTEwQCASm7d3TuLyQjoonGZaVDnW9Z4L1N/HCOtYiu2XR/X4/fr8fCLru4cG4NC59pImdJEgpGRsbY2xsjKKiolkJFLrXTiYsoquqaq/J8h5CiW4F49JEv3SRJnYSYOWmx8bG6O3tpbi4eKGXBBBxby+lxOfz4fP5gOlEt6LuaVz8SBM7AYTnpi1LuVgRDdFVVZ1m0dNEvziRJnacCM9NWxHqi6lENBLRTdPE6/XaP0sT/eJEmthxIFJuGkjJvnk+kSb6pYM0sWNAqOttWehQCCEuqRTTTESfnJyks7OT4uJisrOz00RfhEgTO0pY0WXTNGdMY13sFnsuhG45xsfHKSoqsokeGpFPE33hkSZ2FJjJ9Q7HpU7scFgktzyXUItuXaPQzrU00ecPaWLPgvCy0LlyvpcTsSMV14RadOtvDMOwO9kAu1gm3aKaWqSJPQOs3PRsrnc4LraoeCKIpmou0h49lOhWxVya6MlHmthhsPK6IyMjuN3umGq9LzeLHSuiIXpadCI5SBM7BJbrPT4+TmtrK1u2bInp+Gij4mNjYxw9ehSHw0FhYSEFBQVkZWVddDdxouuNRHRd1+2CnzTR40ea2FMIdb1DGypiwVwWW0pJZ2cn7e3trF27FiklQ0NDnDlzhomJCXJycmyiZ2RkJPJ2Uo5EG1giIRLRA4HANKJLKcnIyMDhcKSJPgsue2JHyk0rihJXPno2Yuu6TmNjI0IIdu3ahWmamKZJdnY21dXVSCkZHx9ncHCQEydO4Pf7cbvdFBYWkp+fj9PpTPStJhWpIHY4rD146DmPHTvG0qVLbQ8nvKElTfQgLmtiz5SbjpfYMwXPxsbGaGhoYMmSJVRVVQFc8PpCCHJzc8nNzWXJkiWYpsno6CiDg4N0dHRgmiZ5eXkUFhaSl5eHpi3sRzcfxA6H9eC0WlCtzy+8oSWtLnMZE3u23HS80e1IFrujo4O2tjY2btxIbm5u1K+lKAr5+fnk5+cDQYs/MjLC4OAgLS0t9u8LCwsXJGC3UEHC0AdKJItuEd3yvkJ70S8nol92xI4mN51Iaah1wxuGQWNjI1JKdu3albCF1TSNoqIiioqKgKCYwvDwMN3d3UxMTHD48GF7f56TkzMvN/BCkMQ0zRnrCUKJfrmLTlxWxI6mLBTid8Wt1xsfH6ehoYGamhqqqqpSQgCn00lpaSmlpaWMjY2xZs0ahoaGaGtrY3x8nKysLAoKCigsLCQzMzPpa1gIVxxmJ3YoIolOwOVD9MuG2OFplNluykRuWL/fT319fcyud6LIyMigoqKCiooKpJRMTEwwNDREc3Mzk5OT5Obm2kR3uVwJn2+hXHHroRwrLjd1mUue2LOphSYThmFw/PhxDMNIiuudCIQQZGdnT4u4j42NMTg4SGNjI7quT4u4OxyOuM8z30iWpxCruszFRvRLmtjxlIXGA8v1rq6uZmRkZMEj1uEQQuB2u3G73SxduhTDMOyIe1tbG1JK8vPzKSgoID8/PyolmIVyxVP1OUZD9NHRUUpLS22V2MUciFtcd2CSYOWmm5ub7Rs2VTh37hytra22693e3p6ycyULqqpSUFBgXxdd1xkeHmZwcJAzZ87Yvy8sLCQ3N3fWYNWlikhEP3XqFG632/7ZYm5RveSIHSpZFAgEpnUWJROW663relyud7wBulRA0zSKi4ttEUa/38/Q0BAdHefo7NQoKvJRUuImP78Qvz+XwsLLox4+FBbRQ/foi1ld5uLZNEQB0zTx+Xy2DlmqyOPxeHjllVdwu91s3rw5Ltd7MTeLOJ1OysrKqK1dRyCwkcLC1WiaxvHj3Tz/iwEOHDiF3+9nYmJiUb+PZCL8fVr3l6qq9v7bIvr4+Dj3338/zc3NC7TaS8Riz5SbVlU16cTu6uqipaWFDRs2THPLFhqp2PPm5kquf72PjEwHQlSSmytYulSSleWkvr6X5uZmvF6vHXEvKChISsR9McIwjFljD+Gue0tLy4LGWi56Ys+WmxZCYBhGUs5jGAYnTpwgEAgseNQ7VejvFwgBRUXnrVNm1vmv/3zQQV6+yZYtOTgcDjZt2oRpmoyNjdHZOc5vfuOlsvIEZWWZNtEvletkNQdFC6upZ6FwUV91ax89k2RRohbbet2JiQnq6+uprKykpqYmaZZxoSLLkSAlvPBCBqoCb3v7ZMTf1y4xKC426etVMIygVzQxoeJy5VFens+Rwxl4PMUghxgd7QlG3E3ByGgdy5er1NRkJ6S9vpBuv2EYMaW7PB5PmtixIlrJIkVRZpxDPResuu+enh7OnDnD+vXrycvLS2TZceF3qspS06QmxTe1EHDddV4UAYMDCkKRFBScP+fYmKDppIZhBKg/4sQ0lwDw8ssONBVe+zo/O3b6aW9Tcee5KS8P3tSTkwF+/SuVrq4uenubUFWNoqJgpmK2iHskLOSDcC5XPBx+v39BtyUXHbFjyU0nEjwTQtDY2EggEGDnzp0xFXEk6wY0gFcVBQnUJGlLMRsyMqCrS+X0aRVFkdxyi58TJzQURdLXp+CZgKYmjRUrdM62Ba16TbVBS4uG3w+mKdB1QUlJ8JoHAnDoUDbbd+iUlVVRf2Qpfn+AjIxuzp07x9jYGBkZGbbbnp2dPet1W0hix+qKLzQuGmKHj9OJ5kkfL7EnJiYYHx+npKSE9evXx3QzWZY+GTegCjwQCBDN7RTvOaUEj0eQnS0ZG1Nob1MJBGD58uCDZGJC0NGu4vcLsrJMJicV2toEWVlBYufkSLKyTV78tQufH3buDDDhEXT3KNTVGRQWmmRmBi1/ZZWBYSiUlZVTVlbO0JDA6fQwMjJEa2ur7b5aRM/MzJy21mjrxFOBWFzxxbDFuiiIHWmcTjSIh9jd3d2cPn3aLseM9QOKVfesSQhGhGDnDOuci9RdQuCLw5KYJnR2Kug6nDzpYOvWAFVVBgMD8Kc/uli3VkcIqKgwePlPDvLyDVSfICfbZHRUobdvCU0nFY4ccWGYkrFRwfi4gmccysokbe0amZk+1qw5X0eQnW3y29+4qKg0qajQOXDAyY7tgqqqLKqqqmyxiaGhIZqamvD5fLjdbpvoUsoFJXYsFnuhyb3oiR2tpnckqKoadVTcNE1OnjzJ5OQku3btoqGhIekqKqE4B/QrCsdUlQOKArpOk6Kw2jTZMXVeL3BKUVhjmoRuBAaAH2kae3Sdf3M6KS8o4OoY1zk2JvjD7114vYIlS3QO/dl
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 432x288 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"ss_list, c = kmeans()\n",
|
|||
|
|
"plt.plot(ss_list)\n",
|
|||
|
|
"\n",
|
|||
|
|
"colors = np.array([plt.cm.cool(i/(K-1)) for i in range(K)])\n",
|
|||
|
|
"if d == 3:\n",
|
|||
|
|
" scatter_plot(data, colors[c])"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "717e0a58",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__Вывод:__\n",
|
|||
|
|
"\n",
|
|||
|
|
"- При N=1000 naive_kmeans работает за 874 ms ± 77.6 ms ms, kmeans за 7.7 ms ± 265 µs, т.е. быстрее примерно в 100 раз\n",
|
|||
|
|
"\n",
|
|||
|
|
"- naive_kmeans для N=10000 мой компухтер уже не вывозит"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "3f82aca7",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"### Задача 8"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 69,
|
|||
|
|
"id": "46c3a97e",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def HWseq1(n):\n",
|
|||
|
|
" if n > 2:\n",
|
|||
|
|
" a = np.zeros(n+1, dtype=int)\n",
|
|||
|
|
" a[1] = a[2] = 1\n",
|
|||
|
|
" for i in range(3, n+1):\n",
|
|||
|
|
" a[i] = a[a[i-1]] + a[i-a[i-1]] \n",
|
|||
|
|
" else:\n",
|
|||
|
|
" a = np.array([0, 1, 1])[1:n+1]\n",
|
|||
|
|
" \n",
|
|||
|
|
" return a[1:]\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"def HWseq2(n):\n",
|
|||
|
|
" if n > 2:\n",
|
|||
|
|
" a = [0, 1, 1]\n",
|
|||
|
|
" a[1] = a[2] = 1\n",
|
|||
|
|
" for i in range(3, n+1):\n",
|
|||
|
|
" a.append(a[a[i-1]] + a[i-a[i-1]])\n",
|
|||
|
|
" else:\n",
|
|||
|
|
" a = [0, 1, 1][1:n+1]\n",
|
|||
|
|
" \n",
|
|||
|
|
" return np.array(a[1:])\n",
|
|||
|
|
"\n",
|
|||
|
|
"\n",
|
|||
|
|
"@jit(nopython=True)\n",
|
|||
|
|
"def HWseq3(n):\n",
|
|||
|
|
" if n > 2:\n",
|
|||
|
|
" a = np.zeros(n+1, dtype=np.int32)\n",
|
|||
|
|
" a[1] = a[2] = 1\n",
|
|||
|
|
" for i in range(3, n+1):\n",
|
|||
|
|
" a[i] = a[a[i-1]] + a[i-a[i-1]] \n",
|
|||
|
|
" else:\n",
|
|||
|
|
" a = np.array([0, 1, 1], dtype=np.int32)[1:n+1]\n",
|
|||
|
|
" \n",
|
|||
|
|
" return a[1:]"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 70,
|
|||
|
|
"id": "f89b5663",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"79 ms ± 2.52 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
|
|||
|
|
"32.8 ms ± 258 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
|
|||
|
|
"520 µs ± 10.7 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"n = 10**5\n",
|
|||
|
|
"\n",
|
|||
|
|
"%timeit HWseq1(n)\n",
|
|||
|
|
"%timeit HWseq2(n)\n",
|
|||
|
|
"%timeit HWseq3(n)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 73,
|
|||
|
|
"id": "8465b982",
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"686 ms ± 21.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
|
|||
|
|
"52736107\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAFSCAYAAACzGKivAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuKklEQVR4nO3debhdZX3//fc3CSFMMiUoQkIYggQqMkRwQMA6oVRwRFCrIAV9+qC1qIitRR6qVfn5q0MVFRVRKyBwtZgqlYJCQRnDIIUgIQSEBIEkJIEMZDrf54977Z6dnTPsk5yzz0rO+3Vd+zp7DXut+6y1196ffd/3WisyE0mSJA2vUcNdAEmSJBnKJEmSasFQJkmSVAOGMkmSpBowlEmSJNWAoUySJKkGDGXSZiQiRkXEdyNiYURkRBzdwXVPq9Y5uVPrlKTNiaFMGmYRcXFE/KKH8RsSct4CnAK8FdgVuLmfdZ8cEUsHUt52DeWym9YxNiI+FRF3R8TyiHgmIm6NiA9HxJZDuW5JGmxjhrsAkgbVPsCfMrPPMLapiYgtMnN1y7ixwDXAwcA5wE3AYuDlwJnAg8ANHS2oJG0Ea8qkTUhEHBkRt0XE8xHxVER8tQonRMTFwFeBSVUN26NNr7k1IpZGxJKIuD0i/qxq2vwhsE01f0bEudVr3h8Rd0TEcxHxdERcERG7tZTlmIj4Q1WWm4B9m6b1teyxEfHliJhb1W7dERFvan5tNf9bqrKuAt7E+j4OHAW8PjO/kZl3Z+YjmXk58Crgrmp5W0bE16rt9Xy1LY7oYX2vq7bt8oiYERGHNM3zp4g4sWn4t9W2GVMN71MtY/f+tl8UsyPiky3bc0q1jEPoQURMjIifV7WBy6tt31ym3SLisohYVD1+GRFTWpZxVkQ8Wb0XfhwR5zbeJ9X09Wptq3nuaxl3SkTMrLbnrIj424gY1TQ9I+L06v9eFhFzIuL9Lct4cUT8NEpT+/KIuCciXts0/a0RcWe1jkci4guN97q0uTKUSZuI6kv9P4G7KbVDpwInAV+sZvkb4DxgLqXp8uVVaPg58FvgZcDhwNeAtZSmzY8Dy6v5dwW+Ui1rLPC56jV/AYwHLm0qy0TgKuBa4CDgX4Dzm4rb17J/SAlT7wX+DPgR8B8R8bKWf/nLwGeB/YDbetgk7wOuy8wZrRMysyszn60GzwfeA3yIst3+B/hVROza8rIvAmcDhwALgZ9GRFTT/hs4uvrft6bUxq0EplXTjwYezsy51XCv2y/Lve1+QGlmbvYh4J7MvKuH/xXgAmBr4LXAAZTtu7ipTNcDz1O27SuBPwHXVdOIiBOAz1flOoRSk3hmL+vqVUScBvwTpXZyKvAJ4NPAX7fMeg7lvfcy4GfARRExqVrGNpRtOhl4G/BSynu3sY43AT8Fvln9rx8C3lWtV9p8ZaYPHz6G8QFcDKwBlrY8lgMJTK7m+wLwEDCq6bUnU8LB1tXwJ4FHm6bvVC3jqF7WfTKwtI0y7lctZ/dq+J+AWUA0zfPZlvKut2xgb6ALmNQy/irggur50dVy3tlPmZYDX+9nnm2AVcAHmsaNBh4GPt+yvjc1zfPqlv/3I8CD1fPXAw9U++0z1bh/Bb4/gO33ImA18IqmMs0DzuhjGfcCn+tl2oeq90a0/J8LgROq4ZuB77W87rqW98vFwC9a5jkXuK9p+DHgL1vm+Tgws2k4gS82DY+p9tf7q+HTgOeA8b38PzcC/9Ay7m2U4yJ6eo0PH5vDw5oyqR5upNQ4NT/e2zLPVODWzOxqGvdbSq3MPj0tNDOfoXzRXlM1Z53ZqK3oS0QcUjWV/TEingMatVGN1zbKkk0vu6W/5VJqaAKYWTWhLY1yMsCxlMDW7H9rwJrnjYjvNEa3sb69gS2A3zVGZObaqqz7t8x7b9PzJ6q/u1R/bwD2rWrXjqbUSt1QPYdSO3VDU3n73H6Z+STwC0qYAjiGEqB/2sf/8nXgsxFxS0R8PiIObZp2KLAn8FzTNl0C7Ej3dp3K+vuonX32vyJiAjAR+G7L/vsS6++//92embkGmE/39jwYuDczF/SyqkOBv29ZxyWUkP2igZRZ2pTY0V+qh+WZObt5RETsMIDXZ68TMk+JiK9RvviPA74QEW/LzGt6mr9qWrqGUovyl8DTlOa3mygBcGOMqsr6ckpNUbMVLcPLmp4f1PS80Sw5ixI0NlTrNlvdw7RRAJn5h4h4ktJ0eDQlIN0BfDMipgK7U4WyAWy/7wOXRMTHKeHs3zNzUa+FzfxBRFxDOcP29cDNEfHFzDy3Kuc9wIk9vPSZ3pbZgy7WD7tbND1v/JD/CP2c2cv6+zdpv8vMKOD/A67oYdr8NpchbXKsKZM2HQ8Ar2juUA0cQWmee7ivF2bm7zPzy5l5NCU8fLCatIrSzNVsP0qI+LvMvDEz/0B3DUdzWQ5v6nMF8IqWeXpa9t2UL/0XZebslse8PsrfPN/T1ehLgNdHxLTW+aNcr+0FlO2yitIc2Zg2mtLnamZv6+vFf1Nq9KYBN2Tmo8AC4CzW7U/WzvYD+BUlYH6EcgmTi/orQGbOzcwLM/MESp+t06tJd1FqSxf0sF0boewB1t9HrcPzKf3/mh3UtP6nKLWIe/ewntm0727gwIgY38v0u4D9elpHVesmbZYMZdKm4wLgxcAFETE1Io6lNBt9MzOX9/SCiNgzIr4UEa+KiD2qs9sOpDuQPAqMi4g3RMT4qlP4Y5R+amdExF7Vev6xZdHfoXTS/lpEvCQi3kUJF83WW3ZmzqI00V0cEe+qlj8tIj4ZEe8Y4Pb4GqX59tqI+FhEHFT9v++oxh+SmcuAbwNfjnI259Rq+IXV9hyIG4ATgNmZOb9p3PtZ99Ib7Wy/RjPqRZQTDOYBv+5r5RHx9ShnvO4VEQdRaj4b+/GnwFPAzyPiqGo7HBkR/ze6z8D8OvDBiDgtypmen6Gc+NHsN8DBEfGhKGeUnkVToK18DjgryhmXL4lyJu8HquW16xJKDeLPI+I11f90XHSffXke8N6IOK9a/n7V++X83hcpbQaGu1ObDx8j/UEPnaur8dNo6jhfjTuScibiSsqX8FeBLZumt3b0fyHwb5Qv/ZWUwHA+sEXTPN+m1PgkcG417j2UWqbngdspl6RI4Oim1x1LOYPveUqfrff1UN6elr0FpfP4HEot1pPAdODQavrR1fw9dgJv2UZbUs78+z2l+XMRcCvwYWBs0zxfq7bXymr6EU3LWG99lMCZwLSmcY3O+t9sGndyNe79LeXqd/tV8+1RjT+njf/1Xyid+Z+n1GhdBuzWsq9/SAk7K4FHKKGv+f/6TDW90Ufr3Ob3SzXPuZQzN5dQgus/0dTRv5rnJEpt1vPVNv8tcGLT9ATe1fKaR4FPNg3vTjkrczHlJIC7W95fb6Q0+S6n1CjOoI8TIXz42BwekdlrVxRJ0hCKiMMpgXavzHxsGNb/SUrQmdzpdUtanx39JanDotwCagKlWfPfhyOQSaof+5RJUuedBPyRckLAgC/gKmnzZPOlJElSDVhTJkmSVAOGMkmSpBrY5Dv6jx8/PidPnjzcxZAkSerXnXfeuSAzJ/Q0bZMPZZMnT2bGjBn9zyhJkjTMIuKPvU2z+VKSJKkGDGWSJEk1YCiTJEmqAUOZJElSDRjKJEmSasBQJkmSVAOGMkmSpBroWCiLiIsi4umIuK+X6RER34iI2RFxb0Qc0qmySZIkDbdO1pRdDBzTx/Q3A1Oqx+nAtztQJkmSpFroWCjLzBuBZ/qY5Xjgx1ncCuwQEbt2pnR9ePRR+N73YOHC4S6JJEnajNWpT9luwONNw3OrceuJiNMjYkZEzJg/f/7Qluqee+D00+Hxx/udVZIkaUPVKZS1LTMvzMxpmTltwoQe7+kpSZK0SalTKJsHTGwa3r0aJ0mStNmrUyibDnygOgvzFcCSzPzTcBdKkiSpE8Z0akURcSlwNDA+IuYCnwO2AMjM7wBXA28BZgPLgVM6VTZJkqTh1rFQlpkn9TM9gf+3Q8WRJEmqlTo1X0qSJI1YhjJJkqQaMJRJkiTVgKFMkiSpBgxl7coc7hJIkqTNmKGsPxH
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 720x360 with 1 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {
|
|||
|
|
"needs_background": "light"
|
|||
|
|
},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"n=10**8\n",
|
|||
|
|
"\n",
|
|||
|
|
"x = np.arange(1, n+1)\n",
|
|||
|
|
"%timeit y = HWseq3(n)\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.figure(figsize=(10,5))\n",
|
|||
|
|
"plt.plot(x, y/x, color=\"red\")\n",
|
|||
|
|
"plt.ylabel(\"a(n)/n\", fontsize=fontsize)\n",
|
|||
|
|
"plt.xlabel(\"n\", fontsize=fontsize)\n",
|
|||
|
|
"plt.title(\"Hofstadter-Conway sequence\", fontsize=fontsize)\n",
|
|||
|
|
"\n",
|
|||
|
|
"print(y[-1])"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"id": "af6b8830",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"__Вывод:__\n",
|
|||
|
|
"\n",
|
|||
|
|
"- $a(10^8) = 52736107$\n",
|
|||
|
|
"- $jit$-версия работает лучше всего"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"metadata": {
|
|||
|
|
"kernelspec": {
|
|||
|
|
"display_name": "Python 3 (ipykernel)",
|
|||
|
|
"language": "python",
|
|||
|
|
"name": "python3"
|
|||
|
|
},
|
|||
|
|
"language_info": {
|
|||
|
|
"codemirror_mode": {
|
|||
|
|
"name": "ipython",
|
|||
|
|
"version": 3
|
|||
|
|
},
|
|||
|
|
"file_extension": ".py",
|
|||
|
|
"mimetype": "text/x-python",
|
|||
|
|
"name": "python",
|
|||
|
|
"nbconvert_exporter": "python",
|
|||
|
|
"pygments_lexer": "ipython3",
|
|||
|
|
"version": "3.10.0"
|
|||
|
|
}
|
|||
|
|
},
|
|||
|
|
"nbformat": 4,
|
|||
|
|
"nbformat_minor": 5
|
|||
|
|
}
|