example_tikhonov.ipynb

{"cells":[{"cell_type":"code","execution_count":1,"source":["import matplotlib.pyplot as plt\n","import numpy as np\n","import torch\n","import random\n","from skimage.metrics import structural_similarity as ssim\n","from methods import RMSELoss, TikhonovTomography, GCVloss\n","\n","plt.rcParams['font.family'] = 'sans-serif'\n","plt.rcParams['xtick.direction'] = 'in'\n","plt.rcParams['ytick.direction'] = 'in'\n","plt.rcParams['xtick.major.width'] = 1.0 \n","plt.rcParams['ytick.major.width'] = 1.0 \n","plt.rcParams['font.size'] = 9  \n","plt.rcParams['axes.linewidth'] = 1.0 \n","\n","torch.manual_seed(5)\n","np.random.seed(5)\n","random.seed(5)\n","D_DIM = 12  # detector dimensition\n","F_DIM = 64  # reconstructed image dimenstion\n","NOISE_RATE = 0.125 \n","device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")"],"outputs":[],"metadata":{}},{"cell_type":"code","execution_count":2,"source":["# d = Rf\n","R = np.load(f\"R_{D_DIM}_{F_DIM}.npy\")\n","f = np.load(f\"f_{F_DIM}.npy\")[:, ::-1]\n","d = R.dot(f.flatten())\n","ls_image = np.load(f\"ls_image_{F_DIM}.npy\")\n","model = TikhonovTomography(R=R, d=d, ls_img=ls_image, device=device, noise_rate=NOISE_RATE)\n","train_y_img = torch.tensor(f.flatten())\n","train_y = model.d.detach().to(device)\n","\n","model = model.to(device, dtype=torch.float32)\n","train_y_img = train_y_img.to(device, dtype=torch.float32)\n","\n","loss_function2 = RMSELoss()\n","loss_function = GCVloss()"],"outputs":[],"metadata":{}},{"cell_type":"code","execution_count":3,"source":["optimizer = torch.optim.Adam([\n","    {'params': model.parameters()},\n","], lr=0.1)\n","\n","model.train()\n","training_iter = 100\n","prev_loss = torch.tensor(20000.0).to(device)\n","\n","for i in range(training_iter):\n","    optimizer.zero_grad()\n","    f_out, d_out, omega_sum = model()\n","    loss = loss_function(d_out, train_y, omega_sum)\n","    loss_img = loss_function2(f_out, train_y_img)\n","    if prev_loss > loss:\n","        prev_loss = loss\n","        prev_loss_img = loss_img\n","        prev_f = f_out\n","        for n, p in model.named_parameters():\n","            print(n, p)\n","    loss.backward(retain_graph=True)\n","    optimizer.step()\n","\n","for n, p in model.named_parameters():\n","    print(n, p)"],"outputs":[{"output_type":"stream","name":"stdout","text":["alpha Parameter containing:\n","tensor(1., device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.9000, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.8001, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.7005, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.6013, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.5027, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.4051, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.3091, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.2153, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.1249, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.0394, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(-0.0250, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.0113, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.0033, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(-0.0027, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(0.0012, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(-0.0004, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(-0.0001, device='cuda:0', requires_grad=True)\n","alpha Parameter containing:\n","tensor(9.6650e-05, device='cuda:0', requires_grad=True)\n"]}],"metadata":{}},{"cell_type":"code","execution_count":4,"source":["out_f = prev_f.detach().cpu().numpy().reshape(*f.shape)\n","rmse_loss = prev_loss_img.detach().cpu().numpy()\n","ssim_loss = ssim(out_f, f,\n","                 data_range=f.max() - f.min())\n","plt.axes().set_aspect('equal')\n","plt.title(f\"RMSE:{rmse_loss:.2f}, SSIM: {ssim_loss:.2f}\")\n","y = np.linspace(-0.45, 0.45, F_DIM)\n","x = np.linspace(1.35, 2.25, F_DIM)\n","X, Y = np.meshgrid(x, y)\n","plt.ylabel(\"Height (m)\")\n","plt.xlabel(\"Major radius (m)\")\n","plt.pcolor(X, Y, out_f, cmap=\"bwr\", vmax=1., vmin=-1.)\n","plt.colorbar(orientation=\"horizontal\", shrink=0.45)\n","# plt.savefig(f\"L2_{D_DIM}_{NOISE_RATE:.2f}.png\")"],"outputs":[{"output_type":"execute_result","data":{"text/plain":["<matplotlib.colorbar.Colorbar at 0x7f4e6c6f0110>"]},"metadata":{},"execution_count":4},{"output_type":"display_data","data":{"text/plain":["<Figure size 432x288 with 2 Axes>"],"image/png":"iVBORw0KGgoAAAANSUhEUgAAANMAAADpCAYAAAC+5gISAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deZylRXX3v+d2T09Pz0zP0DMMw7AIBIkaY0jkjUkUiWbRkGiMmpi4RGJcyCZZzWLEQY2SYN74GjWLS1BcE+OuL0ZUkIDBKCKaqMAnICEIiDAM4zDM0id/PHXuPbfu6ec+3X37dt+e+n0+/Xmerqqn6jxPd9U5deosoqoUFBQsHq3lJqCgYLWgTKaCggGhTKaCggGhTKaCggGhTKaCggGhTKaCggHhsJ1MIqIico2IfElErhWRJ7q6nan+0a7sOBE5JCJvT7+LiLxKRP4zPf8VEXl6qjtLRO5K/dvPC+eg40Ei8jkRuU5ELhWRHXO0e1tqc42IfEJETpxvH1l/kyLyJhH5j0T/l0Tkx1LdFhF5b3qnL4vI50XkoanuQhF5btPv1ICOJ4nI10XkhkTPeNDmOBG5JLX7soi8WUTWproxEXldKv+6iPxxk3GXBKp6WP4ACoyn+zOBe93vO4GrgTe59n+Syt6efv8F4BJgIv2+AXhguj/L2jWg4zLgKen+HOAdc7R7AjCW7n8N+Nh8+8j6exFwIdBKv28Bjk/3bwBe5toeA2xL9xcCz236nfrQMA3cCpwACPA+4HlBu6OBR6T7FvAu4EXp97OBDwBjwDrg88Ajl+N/6rDlTBk+TTUZZlzZx4FHichk+v2XqP6IhmOAO1V1P4Cq7lHV6+czqIgcBTyY6p8I4M3Az4qI5G1V9cOqeij9+jng+Pn2keEY4HZVnU39f1tVb3Z133Rj/4+q3jFHP/2+Ux0eD1ypqjdpNTP+Hvj5vJGqflNVr0r3s1QT5vhU/b3AJ1X1kKreB3wm0TB0lMlU4eeAy7J/mEPAR4EnicgjgBuBb7n6fwROS6LFG0Xkydk/8JmZmPdEABF5ooi8KbU5FvhvNVapugfYC2zrQ+9vAh9ZZB9vAZ6daPtrEflxV/c64FUi8lkReXV6/7nQ7zuRxohEz+OAb7jfb05lc0JE1gHPofP+XwSemMTWI4DH0ZloQ0WPfHqY4fMisgnYAfxYUP9W4FVUf/C3UXEvoFqtReTBwKOARwOvAX6SSuyASgx7Zt6hqn4I+NBCCRaR36XiRL++0D4SHV9K+67HAKcD7xGRC1T1fFX9hIgcT/VNfhT4pIg8X1XfOUd3c36nNNapi6HVICJjwDuBS1T14lT8D8DJwL8BdwFXAFsHMd58cbhzptNU9UTgJcA7nKgCgKpeCxwF/BTw4fxhVT2gqp9W1fOoxJOnz3P8W4DjjKOJyAZgimxlN6SN/7OBM5NIM+8+MvrvU9WPqeofU03Op7u63ar6flU9B3g58Is1/dR+pxr8N/AA9/vxqawH6f3+gWpv+9tu7EOq+keqeqqqPhY4AHxtHjQMDIf7ZDJcANwOPD+o+33g91X1fl8oIg9Pq7fh+6lEnMZQ1dup/vBPTkW/CnzI9jHZeE8Hfg/4SVW9q2kfIvJJEfnBoL9HiciRruhUo19EflxENqb7MeBhDd4t/E59cDHwSBE5If3+fOC9c7R9A7Ae+BUTaRN960RkOt0/BHgq8LfzoGFgONzFPABUVUXk94F3isjfZXWXzvHYVuBvEic4CHwbeJarP1NErnG/f1RVX5z2Tk9U1eem8rOBt4rI+VSarWfYA+n5M1X1VipR6lbg44kJ3auqp9f1ISIt4PuouFeOE4HXOVX0jcBvpPtTgdeIiFJpyb4IvHSO7wDUfqf8Pfwzu0XkN4F/SZP2Mirug4icRqVRPFNEHpne8T+BL6T3/7Sq/g5wBJUYqlR/h2c7RcpQIW6SF6wyiMj3Ab/lJm7BEqJMpoKCAaHsmQoKBoQymQoKBoQymQoKBoQymQoKBoSRVI0nNWhBwZJCVfvZN3ZhJCcTgM72nGuueCjV32b//k6Z3fsye7VWIDfs3Vtdd+/ulNn9Xekod8+e3r62b++UnXJKdT1m5r5O4c03d3dy8GCnbvPm6rrVWelMJmORiYnqOu7+lew+eoERgSyA9tF924KCFYaR5UwrHfftq7jQvn2dMluwjbv4es8IDPPlTMaRfP9RH8ZU2lwI4I47ujszjgMd9hYRZ3X9JIW6PlYJVtfbFBQsI5aFM4nIWVRGjUpl7nJ10OY84BmqevKQyZs3Ds129qnGae68s7ru2tVpZ5zJc6towTamEC3ctrfyz1k74zh++7IhOUP4PdOWqbRXusGZ7BlnMhY5Pd07gOdW+b7I1WlrrIduMYIjVp1fRxRDpz45cL0Q+CEqj86LqHyCfJujgFOGTVtBwWKwHGLeDwKXq+p+Vb0R2GjBMRxeQuVsVlAwMlgOvroFuNv9vosq9sI3AUTkgcAGVb22LoyBV12+9Nxz2blz51LQ2gsvX6X71viadpFJSaYMuMVJUibh+C5MstngfFPtfmqqd/iovbWzfr00ZlrtY491nVx7XXW92XkqmDwaDZAPDm2Z8sBsJdLtD0RX33zd5Hh3JXRkVvtokXp9hES/5aD0LmCz+31TKjPsBM7t18konjMVrG4sx2S6CniFiKyhCuG0J/POPAl4feJKR4vIa1U1jDk3VNjkDU5cxa2eU1PrgW7uYDBG4NXgts/3a8Nkl/N8N4eye19m7e0q+9xhrHGcq51f3g03VFevGrf3ytmch5MGjCMZB44Onbv1FZWUsca/nLHqnENVD/R2En3UFYShTyZVvVtE3kDlVanAOSJyKvATqnqBqv6wtRWRG1bERCooaIBlEUhV9S1UoaY8rgnarXi1eEGBYXR2d8sNE0H8OYmZGjiRaM101W7z5k1A91mRSVWRSOSPdQwmEc240Jhbpg9UN3aQBXBrZpTnTSDMosEfeFm7OrEqOlNysEftGr1TdBY2Pd1R3I5NzM7dSYQVLuYVC4iCggGhcCYHb8lgsBW1XeNXc1v1PSdIq/6RD6pYjW28/aPRKu6Ra6fb3Ag6yoP/+q9O2W23ddPRzw7OBvDs0BQPVrZ5c2979+6TG7rfxX8WY96eiVu9J21mZh0AEnGkRM+h8Q4n27unmxzPMNeML79XTuFMBQUDwmHPmZRezhHJ+mujw0Nbev1+xLhDsnWbnj6qXWWcxov+VubPSI05tN2H/OHqdenA9WsuaOmtt3YT7rmK3fcry0+Ko/d1HNiOA6amNgIxZ+rnt2XfYaON7T78AaqD8H2O6edc0GvZTQIQlo9DFc5UUDAglMlUUDAgHLZinikbvHgym2lqoSPmjY9XJ/5jkfjj5UJ7OD3oLRV27Ojt30S5HTt6263dnWLvezHvppuqqzf6M5276dC9zBgNsC1lm/FinslMNZYekdehWTSMO/vECLkvIXT0N7PJmqLlXDciV5McgZkkrVZHbB+2yFc4U0HBgHDYciZDxFSiFc/qujiT7aD9Tti4QtIijLuz1eNTzgyvpTYmYXUAG/cljmTKBuNGEJ/8Gkcy03Df2Qkn9JYlbqUTHbVzexWP/OKtzI9pL5GOAjY5NbtxGq+XiSzmc02+/y51QWUiLhcZnq8PrO6XEoUzFRQMCGUyFRQMCIetmGeiQiRGRAoIK1vrD4lMvPOWBCnYwj37KhHKS0Ymhfl9/1HbknjlRTmzbrAyb4dnso2PYWdi3kkndV+hLeZ9Z3xTu2jPru6uAKamqo37pg3pnbyywRp60S+3WnCdTW/e0tMkir6UG1ZEOo8IUUgKu/cKH1NGrJscjiKicKaCggHh8OJMbvU0BaqpvKFHq91VZpicdFbPxh22dRKbf6dVWQTsSszEW3xbhKA1u7/dKbw2qbi9rZ2pwnNXct+hV38byzOO5JQNd++r7N/udL6BtupH77lpOhX2C+6XW384VjKW2m3d2rH+iDTuNn5k3xeV5fD0G2eKOJpxqLUTS8uhCmcqKBgQDg/OFOm8E8YmOuuJ2XdFq6GteD6Od6tV7UP8Im4LtS3cppkGWHPH/1Q3/hDW7n2Zqb9tT+b3R8YFHTdsH8gmjvTtXR1ua13tDWzccvd46Bxmj3UbvlVX//28Sbjv1N2PuQ95dOKoB8bXtcuMpsi9yhCpv5vsv3xZJ8Rf50B3KazMC2cqKBgQymQqKBgQDg8xL0LgW90KwvrmG2cvLpn44DXGJmaY+nvT+Hc6lZFIZ2GJfcemXDDxzsuKpmxwYp5uPRLoiHSensjyIFf3++FN5PIWDaH/hA2Wy7V+AHNahLbiZI0TWTel+4mZjT001qXaiaT2yAKiLvjRmiX4zy+cqaBgQFh1nMk7+7XtzSJjrobpTGw1y3z+uuqiaKztQPk3OetuW6m90ZrR4SPrWyeZYsGX3UdnI78rdRtxITvE7BMHpl1m58OtbZ3+N5o63mtf7B2MQ/kPY515NmGs2luvJ1X+unSdnt7YrjJ6PcPOD3ejyLie++QczLdvK1pag1NEFM5UUDAglMlUUDAgrDoxr/soqWLlpliQ6G2duHcwiQ9eJLLNvEkxtzpLglzZAPCQh1TXtXtTbgLvxBdl5bMNuTeVyHPIOhHwnr1rurqKuo38F6OyyHYusnub2lGFfB6LMgmY6OeVDfaRvIxm7+SVKTZo6nfqhI09Vf49ozOzNo1T3fRDR1qus6IYJApnKigYEFZc5kAReRHwFOAgcDXwQlXtu0uM3NANues5dJQTPlZepP7O99m+zlZvb4zQXnhvzh70hHhWZhty30niUodS5B+/OhsjyA0QIHaos/u69v5+f8Cd7Z03es6U7/i9ciJKm2hlURb3pGAZd8Gwo8BPJh3syeLn+a48ifn3iL7LIDF0zuQyB/4o8EzgtVmT96vqI1T1kcBRwGOHS2FBwcKw4jIHqur1ru39VByqoGDFY8VlDjSIyBlU+Zs+E3WSZw58ybnnAfVGkN3HTNJTlp8pQa+7gtcTmIT2Az/QKRu7+t+rG3Op8Aajdu9FusCi4d69KffRHd00eHqjI7O6wElRdGTfPj+TidwhahNF+feMzvVMbvPioMmvqU5mDzl6xnqam8hnZVHojSC5Ybuuy0nQXso/sEjZbyVmDkREHgacDzxhrv1SnjnwUEkkWLDMWHGZA0XkZKrcTU9R1Tvn6KOCm1CmWvBKhnxV7jrttny0rd72HnmMe+8N8dCHJpr/7bOdwmtSmilbdT0rO/nk3k7Sztm4EfQqOyJXiTr1d7TY1kX5gV7rCc8NjRPMzKxvl4m9l72Lfye7jzQcUbKAjo9EL5E1iAJFeeZpQ7UtVLw2IzLmi/L4zgND3zOp6t2AZQ58F/DbInKqiPxBavIaKs71VhG5VER+etg0FhQsBCsuc6Cq/sxi+o7Un21HsCB862xgKe4PLG11sxXPxyqRq79Q3Vx7bafQDi+j1TZKQ57uZ4NgInXcJNoD5bltPSKnOY/cstqrxu3eq+jN4ru954uM6PwRgJV5Tm339sEdJ5ucrP4unpFZd9F7GlOJfBrb36ofe15k0vFyaFtQMCCUyVRQMCCsHtu8tjq2U9QW76KQNUm02B9IJ5F9l0k1csftnco7A/1ILot4OSUPju9om5joBL43kTKiJ+/Kt7er7L+/U5nG0umOS4XBi3J5dB//yaKIyZuOTe9pKv0ozp4X6TI7PKDzUYPvsmHDWFcTX210e3HcPnuU9tae6wptHYl0RcwrKFgZGF3ONDsbql67EpvsyQzS3HJ+//7q0DY6oPXoscT2LME2316lmps0+/bGpYKAdZMBA4vC1UUq4PWT6bBzd2C4l8YS962mpioVtz8QbaKd9p/7vn0pWqqxBO/cGNkg2seN2Im9sHvR6enqL+k5k9For+dpzbmzrw+TLvQLhbQAFM5UUDAglMlUUDAgzCnmicg64HnAUwEzHLsd+GfgTaoauGgNEbOzscO/l11MHjAxwokYkatB3Ybfut0/0QmAP3FSdb9uomNTltubhTJaENc3EjGjKDz5mVIXcVEepcAsYjyJeVEOgsj2L/h8nXqTq7yywd4vinDkkdv1dSkgervNlSSexihVVs/fMTqEjP6HFog6zvR2qsnzBFV9kKo+CHgicAdw0aJGLShYhZiTM6nqU4Kye4B3p5+VhYgz5azGrTwTKWue1x1ElgfWnZl1RTHypqY6VhQzM0cAsbPa2vFDPXRY9r7dTstuY0YW3IbIzrDWO9CVmXW2p9sejVTLkdKjzejsxn/IXFPgXyIydw9MNizI/vR0x3nT6IjCJEeWIZEg0EbEghepgOirzROR7cDTqNwkBEBVz13UqAUFqxBNVOMfphL5ru3XcKjIDcxy5xVwG53EchxbWduq6rqSl22r2Mk9uzuroXUR5R2zbqNzyCAWClu3VpxgYqLDEfYG2my7jw4nbfE8cNAFoa/b8NREWekKmpJtX6IVvjYli+/MxvIfpiZ5QuiQldqZjZ6vbtKVJ2lNK0kE+4IIMhF7WyCaTKa7VPX/LWqUgoLDAE0m04Ui8ldUVt0KoKpvW1KqCgpGEE0m068DXwGOXWJa5oX79gkHD3bsHVqt6n699+E1Fm7edl5Gi2zFkkw2MXNMu8jEvChunol5XrLMN/JetZvblnkyIpeKSHsbOe+tmco28NEmPxBhIhfv/DGAMZKYtL8m9FPwAj5UNRNrydEOXx3AnvVUR86POboUM21/kposAANEk8l0r6r+2sBHLihYZWgymW4Vkd8EvkhHzLtySalqgL1748iekzs61tFjuVebd1u2Os+Z0vI8sb3DmQxR5sA6BUH0e+QpHe3V8xU4YjShnVxg49bzYH6fjdnmFlF+l0iNHCkP2qy1N19wdxeOc82BOo7tA7B02H6gZIhyykScKfKqnAeaTKYx4OHpB6oJteyTqaBgpaHvZFLVXxkGIQUFo44627yXAq/PIwSJyJHAb6jqziWmrRb79/fPhjeWH4FHcdz8A0kcHJs90C7yTnsQb9ojR73oeMfuvVgYkRaFdOvQ0/2cv7dNu/QR6QyhAiAyu6gzxYjOtqx/l3zbRLooWUBNF2EEIvveXdE7Iq/GutSB0QBLaAFxMfD3IrIVuCmVnQh8C/jzRY1aULAKUWebdxXwZBHZADwoFX9NVffM9cwwMTHRwL8rNyrzu/woAr5pNFzZ5OSaObuIfOBMTW4G01GYODnY4XzmZu8d1yzabOSwVxdgp5bT9Dvdr2MPde0jwoI+IpV+buQe/T0jK43O7z4RA72d5KYkkbWDHyAyTJwHmuyZ9gCfX9QoBQWHAYpzYEHBgDCyMSDmOrHvKs9ZuZcTrMwfVuUbVnqd8SJ7Tm8BYSKfWT54EVD2fqe6ibz9HEx4GZsIDpqaxC6I5MKmybGbRr6sO1BrKOblsfwjAwX/vSNpzbA++nvmipN+zn9L6BwIgIi8J/v9bxY1YkHBKkWdavxY4AHAySLyI6l4DfCQxQ7aJ3PgJPBm4HjgZuBXVbUn712rFatNuzb3ue41UjZEO2LXLl+Uo/C7UbB440xje+7pVEZhdaJTd6uvO7mvO83v5yG3wED5IeoUF47G8fE1PUPmBgp1fp2+naHb7b7qfyxSKNRZcHgsoWr84cCTgB1UsSAA9gOvXsyALnPgDwHHULnAP8o1OYtKa/gMETk3/f63ixmzoGAYqFONfxD4oIg8WFW/OsAx25kDgRtFZKOIrHVpZc4A/iLdfxh4EcFkOnhwjoUm2i8EzoE9gU/8vZO7bT9k10g1HuUAay+Qe4NNQsQlPPJ9TrRnijhTHXeIkr3WIfq4dc9F7R2NrYAz5a8ZabWjISKJoJ3QrN8p+Vyd5cQtAE2enhGRfxSRT4rIp0TkU4sace7MgVF9XtfG9LQwM1P9nH/+zkWSVFCweDTR5v0N8AKq/csg0C9zoK/vySpo2L1bF6t8KSgYKJpMputV9bP9mzVGbeZAqiRoZ1J59p6Zfu9BHjavbZ/mzP6lTswzdwwvT5gywLlqzBxfXU2xcJeb2pE5WLTX7SI6R5MkvE0VEHViWF2qwaaoE4P6iFezoQtGfJ2rLP/eoQ4mEtsiX5a8Tf7sAlCnzfsTG05E3gd8gY4/0ysXOqCq3i0iljlQgXNE5FTgJ1T1AuBC4C0icjlwC1Cs1gtGAnVT8X+y68DQJ3PgfcAvNemnbtUCWJtzpih1tz/ks5XUtcv3rtGhY5SZ3OrW+l1yFNV0oZhvxrt+QeZyrhOt2E1Vx6n9IWfXXcdN6vQs0fl2Y6123TtFZUsVnUhV37qongsKDjM0CUJ5PZWFi1IZ594PXA+8WFW/vLTkFRSMDprsuN5Ldd5zFdUZ0ZOB9wN/B/xIzXNLilarAbuvswyIziDssKJmI9pP+jFxxqTH6e3r23USKQ+aiFz9HNjqzqUi27zIvCC3lJjnZtxHIors8JoGjmxSFr1S7Tlj1FlT15R5oElPP6SqV6rqoaTVOy0FVBngBqCgYPTRNDrRG4B/B/4P8E0RGQeChK7DQ6sVW3Cvm3QOcndmkfi9ebfFLfYr044d1dWZetsCZqpxn8kuz1jjkQdGqvrYCMQ5Z7tQx2miFbVuJ5/36e9rln/Paeq4SoQm5oN++DpDBY/ccj/UwYRRkmo40wAD9zfhTM8CPgEcCVwCPFNVD6rqUxc1ckHBKkPdOdMWVf02sJ1qv3RVqtoO3DrXc8PC+HjMmbqWQys0tuKj6BuX8quWmXo7DmYRQmdmpKf7KG5fvo+KDnS7kgU02Uz0808yNF1t6/pN135q6vki2orlAVL6uW3l7aPzWc9RJR+03z5wCa3GXwz8LvCOrFyBxy5q1IKCVYi6c6bfTdfHDI+cgoLRRZNzppOoQntNAz8NvFBV/+9SE9YPuQIi1HCaPGAKBd/INASe9Zt45zUKqd3MzLo5u/CiXO4pH4qikX92XbD9QJ7pCoqfEIY2bqrEyPqdrxFxPxO3unzBUaCgOpPFKF5hKMHl4l2/f5hF2uY1UUD8PfByYI2qHqSaUAUFBRmaTEVV1WtF+gdZHyZarSztiaUS8StaHvGkn+V0lDYv3U/WpIiJTO6sq7Ut50a/O4jVFzmp5Uuwo8fs3eK9ssXbCyzn+0QzNWv7Ohu6uoXb160Z743fNz6eUsTMkzNFOpI6rt8VOzB/IPKZn+slFoAmnOkuEfkFYK2I/CxVBvaCgoIMTSbT84Dvp/J+/WHg7CWlqKBgRCGqNUl/VyhERGdntZulN4mX0NQ5LBC5LCFznZ4AXNB/c/HwOaHMHCLSWEQhjuzqzr0OpQyJdSmYPPmhUsL6qnGR6CcF5cYFXuRqR4gK/CeiMevCI9f9eSLTwq7/ibq/f5//BdmwAVWd196m7tD2eugJXi1Ue6hT5jNIQcHhgLpzpgfavYh8emTOm+riwzXkTN71PXcA9F3YarhuwmWws7y5t91WXS0ZLnS4lF+Cc/W9R2D014TJdtfZ4tp5p7roqnVu991x6rqvXVn8aqLOtpxypO4EwBBZTITKhjq7xKjjRVo7RGhqfz56smBBwZAxsrHGBe1/ELlAnxyvuc4zz/hF17Y5k9ucKtoM9owj3XRT5wGLxuLptn1RVJbe5QCdhGtGj6cjVxXXhYnz93VxWqJtXWSl3S5rGGdPgtQ5deHQa92N+hnz1bmtD9CPyVC3Z/oEHe/aU0XkX+jsmX5y4JQUFIw46pbu5w6NioKCVYA6BcQ3hknIohGJFkneOTTbq9aOYJKCd60wnUGk1Y5CJm+0G3vAKyBMORGJGDNB4NpEfxRUyb+LjV9jyteFOuWFiYj99vHWrn0U4GXjps543Vlp+9KadzHWzzWliWv6kN3WCwoKGmBkFRA9CHasOUeK4phEOgrPfYwrRAFg7d57w2+sswOMTkQjzUa2oh4MHAyjbqOF2OzkIitzj1xpESUoWD/llLrGee8KMpXltoXRAAEiBUR9u847tXMC++9Y92EaaTjmh8KZCgoGhDKZCgoGhKGLeSIyA7yNKsPFNVTOhurqj6VylW+ln3NUtSfbuyId9wLgwGxyIQj2wfPNSB8dX+Tx8Hz/XbmEtib5KLCra8tOkRlFREAaNDLbixQQdu1ygUh9eCHPYlDYOY9/l0jqaVs37HEvb/KvfZBIhvZl9h2cONia6Ha4bBq9OGpvYqxEDzSMtewtXxaC5eBMLwLeo6qnA+uBx2X19wI/n+qfB/zVkOkrKFgQlkMBcQaVGzxUkWLPAC62SlV1SWC5n253vzYE7dpUN1EoeA6Sq5Pn6iPvzy9utij7NDP37KiU45tOOqm7ke/ME5KH3PGEJMLXb3Z2b5t7nQNzE0RTvAC0JtaSIzIWaDtXRh8rYstWH6m86+LVubL8e/f7W7Q5ZHvM3u7DsRraZsoireaWgzPNUGUEhJrMgCIyBrwW+LOwvtWi1RJaLWHnzp1LQWdBwbywJJwpTYQrgqqPUjkZbqKaSHNmBqSKZf7/VfWSqFJnZ8P41gUFy4UlmUyqeogqm3oPRGQDVUbAd6br+4I2rwa+qaqvm3OM7NykLaaM94p+Jh5MTjqDVOY+f4liC5ho5MMdm6TjJSJr/7CHnQjAOj/LTZS7M4gs7RUVuVLC9WHxHaLglrnk5emO0J2zoPoOayIzkNzaN+o4+mhBBmd1Ymfu7uG7GCOJdHuDMevMM5oaNwdRnkZRzPsL4BkpM+AB4F8AROQ1InKkiJwGnAOcLiKXisg/LQONBQXzxtAVECnkck+4MFX97XT7LXA+B/OE38TmLtsS2G2JW2FNZTw7uSZvFiYfNC4VhUk2Ok4++bvaZUeYA6B3y7BOvAIij//rXmp/wA1zDulpjGgzDuCZoZkGHrV1spsuiJUSdSGLjH7vbp/U4HsdbXkEou5oU7O9Y9rLRGYg9lKeG0ZRiQxBWT8rkX4oh7YFBQPCyNrm5fJt+9AuknujrGSR/J3KJiY6nMlSyERW4/aoX8RvuaW7nVebn3DCFgBOPnVLu2zsztt7O85W9vv2d/Z6xnX8mHZv9oN+S2bjewYSpcexzzExUY11hDXKacthK7znrIk7HBhf16FxV29XNmbE3Camq79B1yGsPRyZzhsX8ns5e4fAa6KvPn0AAA8eSURBVDI6oB3FPVNBwapEmUwFBQPCyIp50H/DmNf33V4mkcJLFraHPvbY6uqlmVwEhF51r1cAmJ+g73/r1qMAmHDKAINt1r3YZhKOFx+t3vr3/ogmFvp9udEdGTTYddu2Ti7e7Skv79rZ+3of6MiH7aoDU5t6aIxiV9RZrVjdRq8lsU6sMnoBjzxgfx+1eVFAFBSsEIw0Z/KoO4Rtbyzr0tFBe+WKzlltNff7clsYvSp6T+Arl/cVxY+P7AEjH7soO05uHF2X4gZixYkt7MblLL2vf6ft2zsKhc2bq3s7LL9/f+e733Vnd19Qr6L3XDOnu9XqKIPW28tHDpVNDnSXICKRR+FMBQUDQplMBQUDwkiLef5coO6cqV1Xt9MljmJkj6yd6A2Av2Wm1fUcdMSZ/DpXWWRSZmLPuslqzM2bnf1YFBQ/yUv37K7aeQXErSmVd2QV4enIz6o86pz3Jielpy9TPPi+okSNdcGDIm+V8Q2Ve8varYFsbLJtt8EhADrZEU8Ndf8nC0XhTAUFA8LIciZFGp9Yt5UTwal3FBTIY00rs14OTCB8rxvTarlxa7Uq3rOnN2B+ZJAdcaax/VWleE1BpBtPHW5KnWxySQAeeErSJDhzh/sOVpt6ryAwbmbdexojuu0+8ruLrOnrrB0iRBGl2txt+gjARSSCkDNp4mR1/yejbjVeULAqMbKcqckqslAZuGv1zPXNkc47imifnvM+VKbijlbsWudGvzwbS7CUNdC7MfK65u3bu6/AurSXOM4diB63I3Gzh1bm4/fs6/gd1e2jIh8qGz6ywwtcnNrfJYqbVxfcZnJDO35uZy9ZYw0e/c9E++6FonCmgoIBoUymgoIBYWTFPIg3j/Nl1ZFatu0CD70ebJF7dhQp3xwN+wQiily2293tbSIDOhpNHvN68MirMbdZg478lUw8Nm3b1q5q33vzj/QS9+2Tnu5z8c0P1TSxY90rRwGRTBnRVKFQaymzQBTOVFAwIIwsZ8pV4004UtTecyEf4bSNfPmsycQObtUMzMeixOpR7JF2MJG8EXSWfc8l7D6yWYt00dbOq9yNOOvf15ka3hOZ1O/rkr/75NaOw2OksDAS2yloHG3f2Ss9ZBsipUT0SgtVHiyWG3kUzlRQMCCUyVRQMCCMrJiXI1JA1LHwOoVFVx9ZDDtvh7d/bqOIULERuRyE4d4OZpF/PA22q3dWDj2Gb15eCiIFhSkS83M0L+ZFWhKrT5oHcX0dYaGO/ItGfh+pzFwr7p3tDUoVRY3uSkyQ0O/v2ATlnKmgYIVgZDnTQjaO/bhWxK3y1SpS3/bbOBsiRtNWgNTttP3ybIN5BYS1j/LNWB+eM0Vu37kpeWTpEbHgyCPROKXP05vHAoSeUEsbZjpKjMjKvOvIogGaHJcslht5FM5UUDAgDH0yiciMiHxERC4Xkb8WkXBpSO3uEpFnDpvGgoKFYDnEPEt2dpGIvIUq2dnFQbs/Bq6s66if2NYbnWj+ho75M/4sKhLbcsPVSOIKLSycPNM2zIwCK1o7r4Cw+2jQ/DnoFhsNeUCLKMqlFwvrUio2sAwBOqJnO1R15wxqcrJXGVH392wiovfbGoyiAuIM4CPp3pKddUFEjgeOBnrSbxYUrFQsB2dqkuzspVRJzp42VydzWUA0VUxE7ZvYcnVZTExUZVG0oYg5tBlCED7Ir4pWbeljxiJO4p3gUpoWe87T03ZudNztQKtq38VMNlcOd+t3pPeLPAe9P7xpCOpsFiP4utxgr2EWP/tWi+Ek8/1/aYIl4UwiMiYi/xb8vIROsjMIkp2JyPcCqqpfrSW8VSWIllarZA4sWBFYicnOHg58t4hcDJwMfEdErlPVz/lGs7PxijIIS/K6vVPU/xr3Fe3eDne7UtzkOVmhvTpHjnSds9LO/mFsQzWAj1O3OzER4zR+8d+ccuCun+oQsi/YFtn4U1NVvzMzR7br1kacxgzwcrU8dPZCUQI390F0vHtf1G/fa/V1f4vFcKtRtBqvTXamqheq6umq+njg7cAF+UQqKFiJWInJznzZzmHQVFAwCIysBUQTzFc1GrXLN7tNlR7tfLp1+aLouG/URUnqqkuq+Sinbe7HGPU1F3JzPa9PWDuzuXsg6FUyeJHORD6vOAmc9wyNQlvPUV/Xfr4YRdX4ULBSlBKFjpVFAywdHaI6ONXgsCAiqn2WXGm16NdmGOhHhykqIs5UF0m1jjNF57NTU9JW2phWOzIaN3iDb7tfczBIKWPEOi6kU510NAZBe75Fztn7cYb5Khnm4lLSarW/xVycr9USVHVerGrVcqaCgmFjZDnTctNQsPoxX840kpOpoGAlooh5BQUDQplMBQUDwshPJhH5uIh8S0T+tKbNeSJyw3LSISJ/KCKXiMilIvLYYdMgIieJyGfS+J8WkWOXgoY01veLyBVpvE+JyElZfSOftiHQ8RpnN/pHix5QVUf6BzgWOAv40znqjwLeBdywXHQAPwW8cjm/BfBq4Nnp/izgz5eQju3AxnR/JnBRVn8+8Kx0/xbg8ctExwPTtQV8FviuxYw38pxJVW/p0+QlwKuWmY5fACZF5JMicpGIbKppu1Q0/AdgHoVHAHfUtF0sHbep6r3p1/uBPEpGX5+2YdChqten62yqy6J/zg8jP5nqICIPBDao6rXLTMoOYFZVfwy4isqLeNi4BHiBiFwLnA28aakHFJH1wCuAC7KqJj5tw6DD6p8B/Jeq3rSYcVb1ZAJ2Ai9fbiKofLbMNf9i4GHLQMOfU4l/D6P6Lq9cysFEZA3wHipx8j+z6lqftiHSgYj8OPArVAvMorDaJ9NJwOuTb9TRIvLaZaLjUuC0dH8asKTKkDkggLnP3sEScgMRaVG5z3xAVT8QNLmMag9Dul62HHSIyCOoFtunqup9ef28x0sbsJGFiLwR+BFgLfAVqlX3J1T1gqzdDap68nLQISJrgTcCx1H5cP2yqt42V19LRMP3AH9HtTdYA7xAVb8yaBoSHU8FLqQTw+PLwEeBI7UKpLMFeBswDVwL/FbatwybDnt/W2R+T1W/sODxRn0yFRSsFKx2Ma+gYGgok6mgYEAok6mgYEAok6mgYEAok6mgYEAok2nAEJETRERF5Fmu7M0icmOf596x9NS1x7ok0bldRP5ykX09TUSe07DthIh8KJ3/rDqsypdaAbgaeCpAOmM6jj52X6r6jKadz/XPKCJj86DRbNd+bz7PBDgbuKjhePupDEqfsMgxVyTKZFoa3A0cEJFtwM8AH7MKEXlMcoG4XEQ+KCKTqfyGdJ0SkX8SkctSu5NT+aUi8pci8nGqSLfW3wki8u8ichHwxpr+zxGRz4vIe0jWD+nZS9L9hSLyqHT/TBHZKRXemfr6tIg82r+kiDwYuF1VDzgaX5vcHT4oIi9Ixr2Xi4iFaPkY8PMD/t4rAmUyLR3+icpa/GnAu13551T1Map6OvC11Mbj+cCXVfUMqgQGf+HqPq+qj1PV67JnTgB+Q1WfE/WfJvVZVCGrnwcc3/AdZoAHAI9W1ccA/5rVP5Re06hLVPWxVFYYU8m49wtUqYMArgO+t+H4I4VVHYRymfEhKkvtu1X1Nuf/9j0i8gqqf7ajgN3Zc98N/HO6vxL4W1c3V76qr6iq9RP1f2JqcxDYLSJfC/rwpjACVfTdZKJ0kYjsBV4G9HN5+WK63gJc4+6X1DJ8JaBwpiVCMpx8P/CGrOrFwEsT5/kQ9ASB+zqVfR3p+nVXN9e+y5dH/d9INcnGRWQj8KCgj7uonAuhSp5gFtdvV9VnAp8Bfid75is4kTNB57i39zwlPbfqUDjTEkJVXx0Uvxt4s4h8HbiHXs70RuBtIvIZqn/G581z2J7+VfUOEXk7lS/VdVSTK8ebgHeJyNOpDD93AduAd4vIIWACeGH2fl8VkaNFZCIpF5rgTCoReNWhGLquEEiVNueU5aZjvhCRXwTWq+qbG7SdAN4LPGkprMSXG2UyrQAkrvFNVf2D5aalYOEok6mgYEAoCoiCggGhTKaCggGhTKaCggGhTKaCggGhTKaCggGhTKaCggFhZC0gRE5Q+MZyk1GwOvENVT1hvg+N7DmTiOjEhNJqVTlcx8eZ876ubrHPjEI7QauEtwcPVjloo/u6utXYbq5nZmcR5p81EIqYV1AwMJTJVFAwIJTJNGB86lM7l5uEEDtf9rLlJqEHOy9bkhDjy4YymQaMSy89b7lJCHHeK16x3CT04LzLL19uEgaKMpkKCgaEUdbm7aKT46egYJC4R1U392/WjZGdTAUFKw1FzCsoGBDKZCooGBBGajKJyM+JyFdFZF9Nm8eLyGfTz+PmajdgumZE5CMp2OJfi4vr5drcmII0XioiL15CWs4SkStF5AoR+YGsblJE3pHofIcFqBwG+tB1VvZ9jhkSTR8XkW+JyJ8GdfP/Vqo6Mj/AFmASuGGO+jGqtI6b08+XgLEh0HU+8Kx0/xbg8UGbkOYB03EEVWjmCapYef+a1Z8NvCTdnwucPaS/Wz+6zqJKXj3s/6dj5xp7Id9qpDiTqn5bVefkSlQx3G5U1V2qugu4id64bkuBM4CPpPsPp99zSAoxfLGInLpEdPwgcLmq7lfVG4GNKdb5fOhcDroAfllE/lVEXi5DCuyvqnUBNef9rUZqMjXAFqo434ZdDCeS6Ewaq27MR2gVYvgPqTKALwX6vb+vH9a3aULXB4EHU/3DPgBonMRgCTHvb7XiXDCkyuRwRVD1UVV9eZ/H76IS7wybUtmS0kX10TdRffRwTFW9M12/JCJ7ReQIVb07b7dI9Ht/Xz+wb7NYuvx3EJF3U8Ulb5RZYwkx72+14iaTqh6iCjC/EFwPnCgi0+n3E+kNLD9wukRkA1Wk0nem6/uy+rVUZ3r70uZ6Mx1ONkhcBbwihTU+Gtijqve7+ssSfdek67CM42rpEpHNSSwHeCzdIaGXC/P/VsPe9C1yw3g6VTD8ven65FT+GuDIdH8mVQ6gzwJnDomuLVQc6nLg9UDL00WVn+nqVH8l8JglpOU5aYwrgNOAU4E/SHXrgHclOt4FTA7xb1dH159RTbgrqBQ4a4ZE0xuB/6BacD+w2G9VLCAKCgaE1aaAKChYNpTJVFAwIJTJVFAwIJTJVFAwIJTJVFAwIJTJVFAwIJTJVFAwIJTJVFAwIPwv87mv91M5hTUAAAAASUVORK5CYII="},"metadata":{"needs_background":"light"}}],"metadata":{}},{"cell_type":"code","execution_count":null,"source":[],"outputs":[],"metadata":{}}],"nbformat":4,"nbformat_minor":2,"metadata":{"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.11"},"orig_nbformat":4,"kernelspec":{"name":"python3","display_name":"Python 3.7.11 64-bit"},"interpreter":{"hash":"31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"}}}