diff --git a/diffphys-code-ns.ipynb b/diffphys-code-ns.ipynb index b278b06..afe15d5 100644 --- a/diffphys-code-ns.ipynb +++ b/diffphys-code-ns.ipynb @@ -640,7 +640,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.6" + "version": "3.8.5" } }, "nbformat": 4, diff --git a/overview-ns-forw-v2.ipynb b/overview-ns-forw-v2.ipynb index ff5c46d..7ce9dcc 100644 --- a/overview-ns-forw-v2.ipynb +++ b/overview-ns-forw-v2.ipynb @@ -8,14 +8,16 @@ "source": [ "# Navier-Stokes Forward Simulation\n", "\n", - "Now let's target a somewhat more complex example: a fluid simulations based on the Navier-Stokes equations. This is still very simple with phiflow, as differentiable operators for all steps exist in ΦFlow. The Navier-Stokes equations (in their incompressible form) introduce an additional pressure field $p$, and a constraint for conservation of mass:\n", + "**TODO , double check here - sync with latest local-colab-notebook version, then update diffphys example**\n", + "\n", + "Now let's target a somewhat more complex example: a fluid simulations based on the Navier-Stokes equations. This is still very simple with ΦFlow, as differentiable operators for all steps exist. The Navier-Stokes equations (in their incompressible form) introduce an additional pressure field $p$, and a constraint for conservation of mass:\n", "\n", "$\\begin{align}\n", " \\frac{\\partial \\mathbf{u}}{\\partial{t}} + \\mathbf{u} \\nabla \\mathbf{u} &= - \\frac{dt}{\\rho} \\nabla p + \\nu \\nabla\\cdot \\nabla \\mathbf{u} + \\mathbf{g} \\\\\n", " \\nabla \\cdot \\mathbf{u} &= 0 \n", "\\end{align}$ \n", "\n", - "Here $\\mathbf{g}$ collects the forcing terms. Below we'll use a simple buoyancy model. We'll solve this PDE on a closed domain with Dirchlet boundary conditions $\\mathbf{u}=0$ for the velocity, and Neumann boundaries $\\frac{\\partial p}{\\partial x}=0$ for pressure, on a domain $\\Omega$ with $100 \\times 80$ units. \n", + "Here $\\mathbf{g}$ collects the forcing terms. Below we'll use a simple buoyancy model. We'll solve this PDE on a closed domain with Dirchlet boundary conditions $\\mathbf{u}=0$ for the velocity, and Neumann boundaries $\\frac{\\partial p}{\\partial x}=0$ for pressure, on a domain $\\Omega$ with a physical size of $100 \\times 80$ units. \n", "\n", "## Implementation\n", "\n", @@ -45,7 +47,8 @@ } ], "source": [ - "#!pip install --upgrade --quiet phiflow\n", + "#!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n", + "# [ !pip install --upgrade --quiet phiflow ]\n", "from phi.flow import * # The Dash GUI is not supported on Google Colab, ignore the warning\n", "import pylab" ] @@ -58,26 +61,25 @@ "source": [ "## Setting up the simulation\n", "\n", - "ΦFlow is object-oriented and centered around field data in the form of tensor object. I.e. you assemble your simulation by constructing a number of tensors, and updating them over the course of time steps.\n", + "ΦFlow is object-oriented and centered around field data in the form of grids (internally represented by a tensor object). I.e. you assemble your simulation by constructing a number of grids, and updating them over the course of time steps.\n", "\n", - "The following code sets up a simulation domain, an inflow object to emit a smoke density, and three field: a staggerend `velocity` grid, and two centered grids for the smoke density and a pressure field. We'll use $40\\times32$ cells to discretize our domain, introduce a slight viscosity via $\\nu$, and define the time step to be $\\Delta t=1.5$." + "The following code sets up a simulation domain with the physical size, an inflow object to emit a smoke density, and three field: a staggerend `velocity` grid, and two centered grids for the smoke density and a pressure field. We'll use $40\\times32$ cells to discretize our domain, introduce a slight viscosity via $\\nu$, and define the time step to be $\\Delta t=1.5$. \n", + "(As the following variables are constants, they have upper case names.)" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 10, "metadata": { "id": "WrA3IXDxv31P" }, "outputs": [], "source": [ - "domain = Domain(x=40, y=32, boundaries=CLOSED, bounds=Box[0:100, 0:80])\n", - "inflow = domain.grid(Sphere(center=(30, 15), radius=10)) * 0.2\n", + "DOMAIN = Domain(x=40, y=32, boundaries=CLOSED, bounds=Box[0:100, 0:80])\n", + "INFLOW = domain.grid(Sphere(center=(30, 15), radius=10)) * 0.2\n", "\n", - "dt = 1.5\n", - "nu = 0.01\n", - "velocity = domain.staggered_grid(0) # alternatively vector_grid(0)\n", - "smoke = pressure = divergence = domain.grid(0)" + "DT = 1.5\n", + "NU = 0.01" ] }, { @@ -86,16 +88,37 @@ "id": "ExA0Pi2sFVka" }, "source": [ - "The inflow will be used to inject smoke into the `smoke` tensor. Note that we've defined a `Box` of size $100x100$ above. This is the scale in terms of spatial units in our simulation, i.e., a velocity of magnitude $1$ will move the smoke density by 1 unit per 1 time unit. You could parametrize your simulation to directly resemble physical units, or keep conversion factors in mind for the simulation units. \n", + "The inflow will be used to inject smoke into the `smoke` grid. Note that we've defined a `Box` of size $100x80$ above. This is the physical scale in terms of spatial units in our simulation, i.e., a velocity of magnitude $1$ will move the smoke density by 1 unit per 1 time unit, which may be larger or smaller than a cell in the discretized grid, depending on the settings for `x,y`. You could parametrize your simulation grid to directly resemble real-world units, or keep conversion factors in mind for conversion. \n", "\n", - "The inflow sphere above is already using the \"world\" coordinates: it is located at $x=30$ along the first axis, and $y=15$ (within the $100x100$ domain box).\n", + "The inflow sphere above is already using the \"world\" coordinates: it is located at $x=30$ along the first axis, and $y=15$ (within the $100x80$ domain box).\n", "\n", + "Next, we create grids for the quantities we want to simulate. For this example, we require a velocity field and a smoke density field.\n", + "We sample the smoke field at the cell centers and the velocity in [staggered form](https://tum-pbs.github.io/PhiFlow/Staggered_Grids.html)." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# velocity = domain.staggered_grid(0) # alternatively vector_grid(0)\n", + "# smoke = pressure = divergence = domain.grid(0)\n", + "\n", + "smoke = DOMAIN.scalar_grid(0) # sampled at cell centers\n", + "velocity = DOMAIN.staggered_grid(0) # sampled in staggered form at face centers " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "Let's define the update step, and plot the marker density after one simulation frame." ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 12, "metadata": { "colab": { "base_uri": "https://localhost:8080/", @@ -115,10 +138,10 @@ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 3, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, @@ -140,11 +163,11 @@ " smoke = advect.semi_lagrangian(smoke, velocity, dt) + inflow\n", " buoyancy_force = smoke * (0, buoyancy_factor) >> velocity # resamples smoke to velocity sample points\n", " velocity = advect.semi_lagrangian(velocity, velocity, 1) + dt * buoyancy_force\n", - " velocity = field.diffuse(velocity, nu, dt)\n", + " velocity = diffuse.explicit(velocity, nu, dt)\n", " velocity, pressure, iterations, divergence = fluid.make_incompressible(velocity, domain, pressure_guess=pressure)\n", " return velocity, smoke, pressure\n", "\n", - "velocity, smoke, pressure = step(velocity, smoke, pressure)\n", + "velocity, smoke, pressure = step(velocity, smoke, None)\n", "\n", "print(\"Max. velocity and mean density: \" + format( [ math.max(velocity.values) , math.mean(smoke.values) ] ))\n", "\n", @@ -217,9 +240,9 @@ } ], "source": [ - "for frame in range(10):\n", + "for time_step in range(10):\n", " velocity, smoke, pressure = step(velocity, smoke, pressure)\n", - " print('Computed frame {}, max velocity {}'.format(frame , np.asarray(math.max(velocity.values)) ))\n" + " print('Computed frame {}, max velocity {}'.format(time_step , np.asarray(math.max(velocity.values)) ))\n" ] }, { @@ -246,7 +269,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 6, @@ -295,43 +318,17 @@ "name": "stdout", "output_type": "stream", "text": [ - "Computing frame 0\n", - "Computing frame 1\n", - "Computing frame 2\n", - "Computing frame 3\n", - "Computing frame 4\n", - "Computing frame 5\n", - "Computing frame 6\n", - "Computing frame 7\n", - "Computing frame 8\n", - "Computing frame 9\n", - "Computing frame 10\n", - "Computing frame 11\n", - "Computing frame 12\n", - "Computing frame 13\n", - "Computing frame 14\n", - "Computing frame 15\n", - "Computing frame 16\n", - "Computing frame 17\n", - "Computing frame 18\n", - "Computing frame 19\n" + "Computing time step 0\n", + "Computing time step 1\n", + "Computing time step 2\n", + "Computing time step 10\n" ] }, { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAACaCAYAAADB9pdsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAArNUlEQVR4nO3deZQcZ3nv8eftnn1GoxltthZbllcsO9gYx9gQg8FgDCTgBEIMCTgXgu+9IQlZzgmQnKznZj1JuEnIDfjGXBtCYkJY7IAN2I6JARsvwrIlS8iSbO3rjGbvnuntvX+o6erfK2kWSTVdXfp+zvFxPVPVXa+mn67umnqfepz33gAAAAAAiEum0QMAAAAAAKQbJ54AAAAAgFhx4gkAAAAAiBUnngAAAACAWHHiCQAAAACIFSeeAAAAAIBYtcznzpxz9G7BfBjw3i9txI7JccwH771rxH7Jb8wTjuFIuybNcf3oyWTaJW5zXRJ3mK5vDS53ZV30fJngUy38kAsHfUwc/KASrC9Xog0qwbbl4NnKwaPLVg6eu6jrvcYVX6oblz525n9Jahw3x+f1xPOo7PzvEmeY8s7G7p8cR5zCD7H5Rn4jbhzDkXbNmePOtUnc07FG4hXtV0i81l0g8fIuPe1Y2BadXgarLBuceZaD87NicGY5GXw05oN4rBA9YDx48FhZTxxHfE7i0cywbu8P6/rSPh1L4UhtuVDSx1YqUzqwhn+mx+X4OT7jVFvnXIdz7knn3LPOueedc39U/fka59wTzrltzrkvuDAbgSZBjiPNyG+kHTmOtCPHkRbOh9emww2cc2bW7b0fd861mtl3zewjZvabZvZl7/09zrlPmdmz3vt/nOG5PH9JRPzK67z3V892a3IczaU8p6m25DeaD8dwpF3z5Hg221tbvmLBu2XdR8/XmZQ/+eqXJG4/v0P33d+pcWfdVNy2YExuhmtjPrjkGc6fLZQk9FN1VzVzwdTYiSAe07g0pPvKD+jl2UMDCyR+6vCi2vK3D+rH9XenNkm8a/wxiaeKejW1ea+IHj/HZ7zi6Y8ar4at1f+8mb3BzP69+vO7zeyW0zNQYH6R40gz8htpR44j7chxpMWs7mrrnMs659ab2SEze9DMtpvZsPe16tk9ZrbyBI+93Tn3tHPu6dMwXiAW5DjSjPxG2pHjSDtyHGkwqxNP733Ze3+lma0ys2vM7GWz3YH3/g7v/dVzmVIAzDdyHGlGfiPtyHGkHTmONJjTXW2998POuUfM7Doz63POtVT/0rLKzPbGMUBgPpHjSDPyG2l35uR4WOZ9qm3Zw+YT9VLb7qEpzXeOd7cvry3ftmK5rPvpt78ocfZCrfl0bcFpRkvYT6WurjMTrAv7qxwTz7B9UCPq6tcH9aDZoF40W9K6ytagXrSzqHF/riDxRXt21ZZveEzv9/S5Fy6X+J+DU7FtYw9JXCodsTSZzV1tlzrn+qrLnWb2JjPbbGaPmNm7qpvdZmb3xjRGIFbkONKM/EbakeNIO3IcaTGbK57Lzexu51zWjp6o/pv3/mvOuU1mdo9z7n+Z2TNmdmeM4wTiRI4jzchvpB05jrQjx5EKM554eu+fM7NXHOfnL9rROeZAUyPHkWbkN9KOHEfakeNIiznVeAIAAOBkaA1aX89lEr9zwS0S/9qlQ7Xly985pc+0uGfaPfmxvMTl/RMSj2/RGrf1L55dW/76vm5Zd9/oRom3Dd+v+/KT044FzUbz9NLWN9SW33v5DlmXPUf7V1pQ++jLWjvp2lt1+/a6ir/udl0X1oe2B+tnUtF9+7bWE2xo5grat/OYHqFBzWfYI9R1ao1ntiuq61y1cFTWfaB1p8SH1+k9og6VXpD4yNhIMNpm7et51KlWowMAAAAAMC1OPAEAAAAAseLEEwAAAAAQK2o8AQAATjutlVvZd4PEd15ytcSvvynoibiss7bsy1qT6Q+PzWkkmX6tj+t9lT7fa68dri1fP3FY1v3uhk6J/2ndb0v81/sflPjgyONzGhuSxVlW4su7FteW+64cl3U+p7WRPq+1j1YKaiWDXpuZBXU9Ljs0p12n9r+07g6Nj6kBDbZfoLXKrr7mMx/UJYc9QbPZYH3Qy7Yj2Nc03FKtg116nfblvH6X1oc+OHapxEPuhxJ7T40nAAAAAAAnxIknAAAAACBWTLUFAAA4zVb1vUHih6+5UOLz33BQYtcSTN/z0fQ+1xm0gginGTqdwhhOabRKMFWwHEzXq5sS6YLH9l6rj/3NV+2R+Np7rpP4lzavlHjr0Fcl9hZMx0SiuIxOab28P8oH16N5WBnSNj/FQ/raTo3o9a2WNs2lzjVRnP2xRTqQfp2i6vsW6voOnT7uRrRtiR0Y1Lhuqq5fovtyBw8H2+r08mOm3oZTdcP2K/Va9b2aOUtbIf3Y4iGJz952tsQvBK9Hqdzc7Yu44gkAAAAAiBUnngAAAACAWHHiCQAAAACIFTWeAAAAp6itdZnEn7jgGokvuGm/xL4S1GG26bUA191Wtxy2kghqzsIatNAcajx9Qev0XF7bPfiCPva6n9T6uD/MXyzxf59aI/FYbuv0Y0VDtbX0Svzyhfnash/X9ilhTefBF7V+cf+o1mku6sxLfG7ncG25o1/bjIR558Zy4UA1DuoyKzs0L8eeivbdvUbfe5Wc1mi2vkzrSd0S/Z0cs+9Q+H6rf65Wfa8uXjQicV82qD916TpV44onAAAAACBWnHgCAAAAAGLFiScAAAAAIFbpmjgMADjjZTJdQaz1ceWS9nujryBOh2u6bpX4HT+1S2JfCGo4wzrNsBdnS936mWo6M8F1hFJQwxnWdE5TE+qCr4bea72aC+vXgt6O737fPom/9Zc/J/Hnpv6ublhB70U0XFuL1mV2t0R1nZW85lHukL72zx9eLPHjg9pr8+xOrQH9+SXjteWOHt3WVms/S3k/mFn54Q0SP3N/v8Sf3rpcxzK5vbZcCY75b+y5ROL3r9G8vOg87V3b+/qgp2jQq9N11fXknabe08ysWNB/11RFf8cVr7WuzY4rngAAAACAWHHiCQAAAACIFSeeAAAAAIBYUeMJAEiclmyfxDf2/JLE//uaqDbokl/WmiQ7S2t9jpENPvrKWu8z8qmNEv/Gt86vLX9h6C5Zl5/S2h+cWZyL6od/55Kg12ZGazZ9QXsFWhCGtWCura5n4rKwpkxr647RMkNfz7An4jTjmCtf0se/baXWrN07fl5teWhc6/TMTm3fOHVt2W6Ji5XoGlVlXJN2ZFRrNjePaZ3mN45ove9ir3n8MxdEedzfq7X5ftkSHVhQx/zk17Se9HXfv1viUmlI4q72c2vLK7teKeu+Of68xDs3XCTx+ydWS3yj3ylx7xVtEtuy6Pfigp6fPnh/FYv6Xj1oOu5KhRpPAAAAAABmjRNPAAAAAECsOPEEAAAAAMSKGk8AQMOFvTe/dvWHJH7j16/RB+Tz0fIT62VV5cmt0+7LdWp9nFusdUq9v3O9xHf+xmRt+e0f+gVZ976NX5B4YnJHsDdq1tKsu2NlbfnGN2u9r5/SbSuTQS4E/foyPUGfz5a6GtHBMV23OKhrDmo6/a7B4LmC6wxBvVxl33C0bfv0NWnHKGndn4zbzC7tG5F4aevFteUh09o6s6DfKOZdW1aPh+3Z6PUtBaXBEwU9lo4Wdf2g0/fEptz9Eg9P3FJbXlHUWns3cETi8n9prvz+Bq3xLJV0+7Af7dKuS2vLP9N3uawbL2qO783pP2TjiP47X34w6HU6pO/PlgVRXWb47vHD+kt8Yt9ZEh92zwSPCIvBmxtXPAEAAAAAseLEEwAAAAAQK048AQAAAACxosYTADDvwj6dj7/m/RJf9bnL9AFfeUjCysGoTqZUnGMNzHDQF+2g1ty4bQMSZxZ31pbfft9rZd1T79f6tWsf+0+JR3Nb5jY2NJUL219XW66Max75oFxxckj/1p8Peh52dmtd2YJLoyfIXLBS1vmuoGeoC+pDg7j0nW0ST7yklWcTw9FYuvv039GxWN9fLYuDHqLhJYygf+nyZaMSt/tOQ5Lo69XhtNdmqRKtrxR023JFX/wFQWpcm71K4vuc9kieKkW1yX5iUta5Yc2bTd/Q2tNHJz5r0/FBvfBocd8JtjR7xSJ9P6zu0b6cyzv0ucJ/d3lC3yPZfPRe1t+Y2eQmfX/9x56gj+fEcxJ7HxTONjmueAIAAAAAYjXjiadz7hzn3CPOuU3Oueedcx+p/nyRc+5B59zW6v/74x8ucPqR40gz8htpR44j7chxpMVsrniWzOy3vPdrzexaM/uwc26tmX3MzB723l9kZg9XY6AZkeNIM/IbaUeOI+3IcaTCjDWe3vv9Zra/ujzmnNtsZivN7B1mdkN1s7vN7Ntm9tFYRgnEiBxHmiUnv7XS5ZfO/lWJr3zwNRJX/vZfJS4f1rpMX1dSk9FyHHPtYc/CoMqmrPU8Pq/1O5WJIB6L+stlv6Q1nJd88mqJf+1tSyX+05f+Sp+rEjTCwymb3xzXXHrrolXRmjbtHzt1ULfdt69P4v0T2rt2cYc2/rzs3KjWODswPM0ozKyodWD+sPYV/Ob9qzQ+oG+aQl3Kv3Kxvj9uWH5Y4lUXa11zS39QXxq837zXOGtRIaAL/iV0vT2+Rh7Hy3WvX7mgr7UPXrD+Vq11bM8GtZAVPY5X6p7bLdLemH6JXrxdP6g7K5WGTzzoo88g0fhkVOO5/ogeh9+2XMcd3jWgr13H3d87/XHcT5aPu2xm9r1nztG4sEHiqaLeYyBt74o51Xg6584zs1eY2RNmdlb1jWBmdsDMzjrR44BmQY4jzchvpB05jrQjx9HMZn1XW+dcj5l9ycx+3Xs/6lz0VwrvvXfOHfeU3Dl3u5ndfqoDBeJGjiPNyG+kHTmOtCPH0exmdcXTOddqRxP98977L1d/fNA5t7y6frmZHTreY733d3jvr/beX3289UASkONIM/IbaUeOI+3IcaTBjFc83dE/p9xpZpu9939Tt+o+M7vNzP68+v97YxkhEDNyHGmWlPzubNc+hJ+4VfsKZj5zQOKpPVrvFv6ZNNtZVxsW1JS5tmDjoGCnEvb9DB8f1NRUctH2lRcnZF375hcl/oMPHZT4//7ZFRIfHHnccHrNZ4470557a3uj+t/yhObNyIDWcD53pE/iZ4f1K9glC7Sv53l7BmvLfW26rV+2WAe2/gUJN9yjNZz/c+vTEu8bfULitpbe2vLOiffKuu6s7qu3R/st9rXqezXTo++/zbuXSTzsot62PmX1a3Fp5HE8W3cRtVzU13aqrO+HQkWPpc9M7ZC4EtR4drRE7x9/0WpZ5/v6JD4wpbXFc1Uqj9eWN7gnZd13Bm6Q+OblQxIvXTgucXe//jsy7UF/06Go5vp7j+pn31/+UGs+9+afkjhtfTtDs5lq+xoze5+ZbXDOra/+7HfsaJL/m3Pug2a208zeHcsIgfiR40gz8htpR44j7chxpMJs7mr7XTvODdSqbjy9wwHmHzmONCO/kXbkONKOHEdazPrmQgCaTTQFJpPRKVwreq+TOBNMH9sz+pjEvqJTq7zpVJG03e4bp9/LO35K4tbff7PEU793j8TFsaAFQ7vmmKtL2ZbFmr+uTeNwmm4mmIpbGdGpTeV80G6lblaVD2bpTj26W+L2my+U+BLTlgAHjam2zcwFvXt256OvUeW85uzB0R6Jw6m1Dw/tl/hQXm9I+paJqO2Iv+BcWec7O4Nx6dT1v960UOK9I9+TuLVuaq2Z2bv7P1Bb7gvaERW9Hu/HJvTzpGdcpx22BHPbHzmkbTL08yX8LEGjtZi+vvXtVKYmNYdHi60SPzmgx86NI1+UuKNticTtrdFUW+vulnVudFTiXMlOSf0U1sEJnZr+lUOX6LiyZ0v86mBK8YWtgxLv363T6tcdjqan/+VLWkby3Jh+1lUq+XCklmZzaqcCAAAAAMBcceIJAAAAAIgVJ54AAAAAgFidwTWerm5J525nsz3hxqJcyUnsfVijQM0C5oPm7er+N0n8J6uvqi3f+gv7ZJ276Sqb1jMXSbjxbl39q+v00PHdsTtry5WK3nYcMDO7duFSid2/fF3i3KHpP458UFyZaYvi4gGt0WzV7g1WmdDioLDtRfgn2MnD+t6qbyGQyQaPDe7w33LjDRIvzDyoGwTvWz4vmotzmqe9rVE+FCb0tR2Y1Fq5oSnNne0lbWkyMqX1wYVCXX1wTuvsXRDbpL4HvjzyLxJ7r9u3ZpdL3F9X13nT2bptb6s+d1uL5my2U/9d+36oNZ3fPjQmccWfYrEe5lV9O5XRXIese2Fc4/+c0vtDOKcH1552zbsF3XWteEa0ptOKwXH9lC+VRXlbKA7Imh9O/afEd+9/tcSbh7XVS/ZFzfG9OW0ptL6uln9gYrOsq1SCVmEpr+kMccUTAAAAABArTjwBAAAAALHixBMAAAAAEKsU13hqP62udu2B9frOn6stf/hifeT1r9ilzxS07N34vPba+vRWnev9lfGvSjw8Xj+/m3oenCxNxHP7tGf05vf1Sdz6N2+IHvnidn2mLS/pU1e0fs5fqTWea299mcSPfOFrEt/ykQ/Xlr8+8vfBU2tNNM4cru4j5i3LtVdZZUjryMbHtB5uJK+1Q4t7NI96s9Hju68K+r/1a0+1TFGPu5n9WnO277vaiy7svzhZimr3lnbrOJb0T+hz33Ovrm/XfooueB+fWdU9za9S0Z6Vi1qj3CqV9G/5kxWNz+4KvkwEJW07xh6VeCz/7tqy239o2nHtf0DHVSxPnGDLo3JTeyV+aCz6jLi493xZ95qlmvMdnVp798IG7c34B+v1/fO93F0SVyrTjw3zTY9C415zrezPqS0PBsfl7x/S7w6D+S3Bc+t7oC2judHeHeWSOzKsDw2+l3SG9fVB/ehcDqZhL/Lx/A6JNxeHJP7huO6rNaufOVMl3b5S1/u8vn/onAeaQlzxBAAAAADEihNPAAAAAECsOPEEAAAAAMQqtTWe7a1ah/l3l7xP4g98KqoBqlys9Wy+76Zpn/vqcS3M+PGXdkr8C798s8Q/uzH6NR8ZezZ4tjN7rjdmb9GCKyT+wc3arLDtj39GH/DJz9cWh76vNQYjI53T7qu7+4Du+/KnJPYf0hz/6p9GtUlXffxnZd2zw58Nnp2cP1O4TFQPtLhde5cV9mhO7htZrHFQS3RgQus2X718d225Mqb1bTahz+1z2jdwcrfWDv3wsH5e/Mc+3Xdva1Sb99ql2qtx5bnayLO4QfvDfWtyo46FOv+mFtZrPX0k+ny/8UI9ti1s1bxz1iZxqaJ1z8WgTmyqrrbYRoK6yIVaYzYwqnHY9zbkvb5ntow9UFv++71vlHWbhi/RB2/tl/DZMX0PPJn7jMSF4vT1qUiWyfJoEGePu2xm9lJhWOJSWWv3wzwrBvd8yNQ/3ZDW3luX1v0v79Ccdk7fT+G+phd+D9HjcqkUNGg23Xf4Xg3X8z3nxLjiCQAAAACIFSeeAAAAAIBYceIJAAAAAIhVamo8W7J9Ev/emtsl/sCntT7IL4v6TmWeekbX7RnUJ6/oXG23TPuy+UvOk/iG/3eBxP/43qhH3H97XueF56a0PhSIaM+3mzreJnHfX2tvzdKffFHiDU8srS2PFbVPYbEy/d+cWoe1XqHjgNY/XLbzfol7/vDNteX3/u2LOo5RrSelr+eZw/uoxu2xAT1uXt6udV87gxrOxwf142mVrrYLd0X9k8/Oam1Q+8Wac8V9mr9bty2V+J6dWkv0QP4RiV/X+tra8ttXah1R15U6sGe+qLV2B8e/qQOn9idVvj98pLY8OqZ5t6JnXOINW/Q4nJs6KHHY47X+uO1fHtyLYon2zlzQoXnmwh6HMyjV1axtGf6qrNs+1qv7DurZjq2HQzObKurr+fRQdIx7/TJdlwlyNuP0uO29HnvzpSMS50ai7zELLtX+sdaux+VLevW7ektWj72Fir7fTs1MtfjU6p8srngCAAAAAGLFiScAAAAAIFaceAIAAAAAYpWaGs8VC66R+Lffu103yGgPxMInH6otD+3QeeRjE9rDLay16+7QPm19K3ZL3PVOrcV414ej2ok/+60bZd36Ke13BfxIJqP1Qr94vtaWuW89JvG6x8+W+OBklNeTZf0bU7GiOR1qzWgdWkdWa3oKW5ZL/JovRDVx71qjvbV+b3uPPpYazzOG91FPt/v2ah/CX7m8T+KF39Oeh4+NHpa4d0Tz6NpFUY3nea/U3nJulda/tS8aljj3lH70fWlU66PHJ/dKXOy7vrZc9vpeeuCzuq//sfVJicsl7YmHZqfHwg3FB2vLn9v2c7LumkXa03Cj36rPVAl6Hgb1v1kXxb5ba4dD57xa++RmHprr17u6fYW9F0sD4cZIsami3ovkvgPRfU+6WxbJurNb9LidyUyfdxOT+yWu1H038Qv0GG89CyRce7XmYfc6/c5Dv9jmwBVPAAAAAECsOPEEAAAAAMSKE08AAAAAQKyauMZTa9SuzrxS4uyt50o89hePSvzcC1GN2kRJfw0Vr88ddl0Lq+O6BnSO++XDL0nc/9FobD/WrnPUnw1eAm/6XDhztbX0SXzd2j0S7/+K1iK/MKa1F4W6Os62oGazp0XrlIJWtZYLakKHgj6gA1Oat2sejXp79S/RWr6O1n4dF3UYZ6TN7mmJH/uMHrOXdGiN2m6/UeKRvPY8zpV/PgoywVF5RHPQD2h/t4OTWhs0GdQ0hTVuD099tbb8wFPaM7Q1q7V3o7kXdCz07UwZfT3rX+//c+gBWffvgxdLvDe/Tp8p6HEY9vGsr/F0E5rTVghq/ldqfdyyrssl3lXYZ8BslIP7MKzLRTXwQ/tukHU5PyxxsRzew6Fi03GZ+hzXx/qgF237q/S4vfpOvQ/GkD037b6QDFzxBAAAAADEihNPAAAAAECsOPEEAAAAAMSqaWs8ndPem2v7tQat8tXHJX580yqJhwq6fb2Sn77HYX3dhZnZcFF/jWPbV0r85nuj2qYVXRfIOpfRnoe+Qo0njso4zdGD+7Sn1b4xrenZm9dehv1tUZ62BXVJI0XddiaFoExj76T+zWrXUG9tec/es2RdrvCNOe0L6bRv+DsSv/lprdm8quMdEg/ntFZ+qqh9PeuP4ZW9WsPpFgZ9CHdrv8TRIP+LpWGJwxrPkYn6uk19M+Sp4TzDRa//wNgGWTOU1X7i5bLmqZnWeJrT7wPt2Wi923dQt+3UGn8fFOrf0K69zT9r3zJgNrwvSlx/fNw29pCsyzj9/ls+psZz+uNjqRQdi90e7fFp/Qs1btfvRG9dpN/rnx0Jvk8Hx3EkA1c8AQAAAACxmvHE0zn3GefcIefcxrqfLXLOPeic21r9f/90zwEkGTmOtCPHkWbkN9KOHEdazGaq7V1m9kkz+2zdzz5mZg977//cOfexavzR0z+82RuY1Mv5A0/qdNktYzotxdWtDqcRZoOZtmGriXIwFbc1aFXx4oT+Wl/xg+jy/6NHBvW5K9pCAA1xlyUwx3NT2j7lvl06hXV1l07L3pvTPKxvCzTsdGphV0u4re57vKQ5PhnMCNs0pFNYLu3tqi1/bY8+Wbk0ami4u6zhOa55kS/o1MHHiv+kW1eCdhFBDu+djOLCAZ0a1t6h21Z0tVWCthWtQeuiY1v+BG8AJM1d1vD8NvNep3SXSkHizdRaIphqm6tv9XbgiG68YonGQUuhm5brvv75cJfElUo4JRIJd5fNW47rsbq+7U+pNCLrwhZAfoZjZSYoLxscj/Jy1W4tp3BtQUlc8EXlbSu0tdUn9i6TOB98h0IyzHjF03v/qJkFRzx7h5ndXV2+28xuOb3DAuYPOY60I8eRZuQ30o4cR1qc7M2FzvLe/6gK+ICZnXWiDZ1zt5vZ7Se5H6BRyHGk3axynPxGk+IYjrQjx9F0Tvmutt5775w74W2rvPd3mNkdZmbTbQckFTmOtJsux8lvNDuO4Ug7chzN4mRPPA8655Z77/c755abWVgQEzvvtTbyG7l1Er97/6USbxzWx6/oiualL23X92B3y/R1GGMlnaG8P69z3MP6t6w7p7a8pXJf8GzT7wsN0/AcD2snPrd/l8Qva9G2PTuKAxJf37K0trymR3N8UavWYZSDj6FMUE+3a1w3eLz8iMTrd0b3NNg98qism6nmAw3T0BwP6+H8TF+Fgtr6jrrDcMcVett993JtW9V5gdYOXf+CxhftfpPEzw99fobBoAk0/Bg+19rglqzWYS7rjuow/VVXyrrKhRdJ7C7SVhQ/ve7rEl+2450Sbxj+3JzGhkSapxyvz+OwpjM0/YG8JaP3W+luq/u+fKV+by9fdrnEblzvF3HN5q9I/LqN75H4G1N/NaexYX6cbDuV+8zsturybWZ27+kZDpAY5DjSjhxHmpHfSDtyHE1nNu1U/tXMHjezS5xze5xzHzSzPzezNznntprZG6sx0JTIcaQdOY40I7+RduQ40mLGqbbe+/ecYNWNp3ksQEOQ40g7chxpRn4j7chxpMUp31woKfaP/0Dij29aLXG7ac/DszqjmqCwyjJX1gvBYe1RMXjAhD61rfNPSfzgjm3Rc09pnR5zznEiPsjZZ4c/K/GWtuUSL+nS+gjvoxrPluBeAoWK1mmENZ5hX89C8IOR/E6JDxWfrn82A063BV0XSnxF33ht2fV368bl4KDcojXLZ6/UWqG3b9N6uZ2TWiM6nt8+l6ECJ2V1z09IfO5lUc9Ev6BHNy7qvSRCHZfq9r+4XD8v/mCKHMfJOLXvrCu7f1ziFWvq+oJ2dso6N6X3AbCSfrfIntcv8XtWa4/Qx4svk3hkYvNchoqYnGyNJwAAAAAAs8KJJwAAAAAgVpx4AgAAAABi1cQ1njrPfKp4QOKnxu6W+LyFWn99MH9lXaT1bp1aDnRM/dt4UD40kNd554XKuMS5qT11ETWdOFmaO5MF7dtW6jxf4uFCVIx8eEr/xtTbqs8V1i0PTul7IlfSDXo7z5H4cPD+A063ta1vlLi7NToQV/ZNyLrwL6p+JCdxuaBbrOnW/F7T/mqJN1D/hnnw4y1hnX6Ud5ntO3TdiNYpuzF9D5SHtNf5Zb1aL3dh++skXk+OYx5c4dZKXCnsqC27rTt040nNYSsW9bED+l37/B798n5xy2slfsqo8UwCrngCAAAAAGLFiScAAAAAIFaceAIAAAAAYtXENZ7T815rdsqmc8MP5aO4LdMq6wqtWt8W9jQcLegPDhe0dqIto/2znIvmnXtPj0OcHs61S9yR6ZV4uBDl2mgx/BuT5nhY4xnWMY+X9P3TmdH+Wc5F/bO8n76/HDAb2azm86qWPonzpei4W9qdl3UtQU1yJUjoyZx+9IV9blf4syTeVDeWcllr64CT1dKySOKlwQ0mJvZF303at+yRdW5oTJ9sSo/RxUNBXNHPgHPcUok3kOOIQUu2T+LFHXrsHTwY9WDu3qL3rcjkgxrPQHGvHvcny/p+WhF8T8mS44nAFU8AAAAAQKw48QQAAAAAxIoTTwAAAABArFJb42mmNT4lr3PFh3w0N7ynoHUVYS1E2Wv9T76k8bhpjacP9g3EI6hj81rHVl+XOVrQt3q5EtQxB888VtCf5IMa6YpRq4z5FR6Hj0xFdcVjOzW/u/J6TC7nNd8HR7R+dCiogS5xDEcDjAfF9nsPL6wtdz8/IOtaB7W+zQf9xId36z0ADk/pe6RILT4aYCLI8V3D0bF4yfNa49mZOyyxD2r3x17SnD4wqTk/RY4nElc8AQAAAACx4sQTAAAAABArTjwBAAAAALFKbY1n2EtwoqRzxUdboh5Yg4U2WVeqaF/PUlBbFPY0zLnxafdFX0PEIcyr8fIhiUczUQ3Q4KTmeLEtrGPW5x4uaL1oPqhjzpUHpx0LcKrCPmt7y0MS78wtqy0v2rdE1i0dmZC4UNKPuh1jCyQ+OKk1oHsz2jORnm+IQ6l0ROKtOe3Nuamu/q1lg9a3Ldmj3ztKwfeYnUf6NM7pvSx2OXIc8SuVhyXeXtCc3zK2uLa8YJsex88e1PeDr+hxfM+Q1uqT482BK54AAAAAgFhx4gkAAAAAiBUnngAAAACAWKW2xjM0PrlP4iM90dzvdq+9fwrFDokrpgVwo0FN56DpPPJwX8B8GM3vkvjQgtW15faS1i3nyloPFJR42ojPSXw4s1v3NaH7AuK2tfw9iZ8cuKW2XLEeWXdWTo/hxaBv7YsT+tH3/JDWMO+ZWneywwRO2sbKf0n87YM/VVvOlRfJuhVjmvNFrzm+dVyP8RuO6L0pdk0+edLjBE7WlvJ3JH788C215ZLvk3Xn5rokDrsrvzihOf7soN6bghxPJq54AgAAAABixYknAAAAACBWZ8xU2/C25fvHf1Bbznfpbfq7s0sl9sEF/rCVxFBu+7T7AuZDmHe7xx+vLY93aauVhW6lxGXTaVjjXrcfHNsy7b6AuA2Pb5T4ay76+Nqx53Wy7vxObZdSCPoF7Qxu6f9s4QGJx/N6TAfmw5Gx9RJ/ue7awPaJG2TdBd3dEhcrmuMv5bUVxbri1yUmx9EI4XH8vrrj+K7862XdBd06nbwU5PjOvJa9PV28X2JyPJm44gkAAAAAiBUnngAAAACAWHHiCQAAAACIlfM+bKQQ486c82bZedvfycpkeoKfaI1npZIzJFl5nff+6kbsuVlyPJvtldh7vQ05OZ5kZfM+6J0wT5olvzvbV0lcKmu7lGJpYD6HgznjGD6TrvbVEpcqeYkLRa3TR9KQ4zMhx5vd8XOcK54AAAAAgFid0omnc+5m59wW59w259zHTteggKQgx5F25DjSjPxG2pHjaCYnfeLpnMua2T+Y2VvMbK2Zvcc5t/Z0DQxoNHIcaUeOI83Ib6QdOY5mcyp9PK8xs23e+xfNzJxz95jZO8xs0+kYWCNVKuMzb4QzQWpzvFwebfQQkAypzPH81J5GDwHJkMr8NjPLTe1s9BCQDOQ4msqpTLVdaWa76+I91Z8BaUGOI+3IcaQZ+Y20I8fRVE7liuesOOduN7Pb494P0CjkONKM/EbakeNIO3IcSXEqJ557zeycunhV9WfCe3+Hmd1h9qNbOANNgxxH2s2Y4+Q3mhjHcKQdOY6mcionnk+Z2UXOuTV2NMlvNbP3zvCYAbPyTjNbcnQ5cZI6LrPkji2J41o98yazcrI5PmHJ+538SBJfL7PkjssseWM7XfltNvccT/ox3Cy5Y0vquMySN7ZGH8OTnONJHZdZcseWxHGR4yeW1HGZJXdsSRzXcXP8pE88vfcl59yvmNk37Wgn2s9475+f4TFLzcycc083qnHudJI6LrPkji2p4zodTjbHk/w7SerYkjous2SP7VTNNceTfgw3S+7Ykjous2SP7VTwPWV+JXVsSR3X6UCOz6+kji2p4zqeU6rx9N7fb2b3n6axAIlDjiPtyHGkGfmNtCPH0UxO5a62AAAAAADMqFEnnnc0aL8zSeq4zJI7tqSOq5GS/DtJ6tiSOi6zZI+tUZL8O0nq2JI6LrNkj61Rkvo7Seq4zJI7tqSOq9GS+ntJ6rjMkju2pI7rGM57bm4FAAAAAIgPU20BAAAAALGa1xNP59zNzrktzrltzrmPzee+jzOWzzjnDjnnNtb9bJFz7kHn3Nbq//sbMK5znHOPOOc2Oeeed859JEFj63DOPemce7Y6tj+q/nyNc+6J6uv6Bedc23yPLSnI8VmNK5E5Tn7PTlJyPKn5XR0HOd6kkpLf1bEkMseTmt/VMZDjMyDHZzUucjwm83bi6ZzLmtk/mNlbzGytmb3HObd2vvZ/HHeZ2c3Bzz5mZg977y8ys4er8Xwrmdlvee/Xmtm1Zvbh6u8pCWObMrM3eO+vMLMrzexm59y1ZvYXZvYJ7/2FZjZkZh9swNgajhyftaTmOPk9g4Tl+F2WzPw2I8ebUsLy2yy5OZ7U/DYjx6dFjs8aOR4X7/28/Gdm15nZN+vij5vZx+dr/ycY03lmtrEu3mJmy6vLy81sSyPHVx3HvWb2pqSNzcy6zOwHZvYqO9q0tuV4r/OZ9B85ftJjTFyOk98n/L0kKsebIb+rYyHHm+C/pOV3dQyJz/Ek5nd1DOT4sb8TcvzkxkiOn6b/5nOq7Uoz210X76n+LEnO8t7vry4fMLOzGjkY59x5ZvYKM3vCEjI251zWObfezA6Z2YNmtt3Mhr33peomSXxd5ws5PkdJy3Hye0ZJz/GG51CIHG8qSc9vswTkUL2k5Xd1TOT4iZHjc0SOn17cXOgE/NE/GTTslr/OuR4z+5KZ/br3frR+XSPH5r0ve++vNLNVZnaNmb2sEePAqSPHj0V+p0ej89uMHEe8Gp3jSczv6r7J8ZQgx4+vmXN8Pk8895rZOXXxqurPkuSgc265mVn1/4caMQjnXKsdTfTPe++/nKSx/Yj3ftjMHrGjl/P7nHMt1VVJfF3nCzk+S0nPcfL7hJKe44nJIXK8KSU9v80SkkNJz28zcvwEyPFZIsfjMZ8nnk+Z2UXVuy61mdmtZnbfPO5/Nu4zs9uqy7fZ0Tnd88o558zsTjPb7L3/m4SNbalzrq+63GlH57tvtqNJ/65Gji0hyPFZSGqOk9+zkvQcb3h+m5HjTSzp+W2WgBxPan5Xx0aOT48cnwVyPEbzWVBqZm81sxfs6Fzk321kcauZ/auZ7Tezoh2dC/1BM1tsR+9StdXMHjKzRQ0Y10/Y0Uv3z5nZ+up/b03I2F5uZs9Ux7bRzH6/+vPzzexJM9tmZl80s/ZGvrYNzityfOZxJTLHye9Z/54SkeNJze/q2MjxJv0vKfldHUsiczyp+V0dGzk+8++IHJ95XOR4TP+56mABAAAAAIgFNxcCAAAAAMSKE08AAAAAQKw48QQAAAAAxIoTTwAAAABArDjxBAAAAADEihNPAAAAAECsOPEEAAAAAMSKE08AAAAAQKz+P9lYtqYmy8aWAAAAAElFTkSuQmCC\n", "text/plain": [ - "" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAABUCAYAAACbU2yrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvDklEQVR4nO29aYxkWXbf9zv3vhd7brX3Vr1NL7MvnBkO6RmSsMRlKFEjWQZBmbRJi/DAsAWYsA17bAKGPvgD5UWGBBkUxiYhSqAt2hBlDmQZJC1zKMv27Oqe7pneu6u7a82lMjPWt9x7jz/cF5FZW3d1d1VmRM77NbIz6mVk5I0TEf933rlnEVWlpqampmbxMIe9gJqampqad0ct4DU1NTULSi3gNTU1NQtKLeA1NTU1C0ot4DU1NTULSnKQf0xE6pSXmpqaO4AABkEQk2CwGInfBcFgIP70ut9T9JpbesN/EOItDdf8HFWUQMzci/ebPtIBsKmqJ68/eKACHrEH/ydramqOCIJIijUduq176KYnOcnDdLTHMm1aJiExQtMYRMCKYPZpeNAoySHqMV4Vp4rf95WrQ1FyShyeUkpymeApmdCn1DGFH5K5XZyfUPohIRSoloC/S8/bv36zo4cg4DU1NTXvDsGS2GVW2md5zPwIZ5MVHltJ6CbQTZRUFCuQiCICBr1GwCGKOET/WVUAwamggFPBaxOvUAahDOAVMg9OYVQq/TKwUzgu2C12ZYPN4kXycpfC9QlhMn3kA7FHLeA1NTULg5gW7cZx7rcf433NVe7tWI43lYaBVHSfcDMTbrlOTG113MYfAtBAUYSgWgk7BASvU1GPQl40hbGzjJxldXKafnGClxodriZv0s/eJCu30JChuAOxRy3gNTULTxUPFlt9jx9rkRtzFGJcFyBUtwNoQGeX/vO8TSUY06CVrHKG45xqWc60la69Vrghivd+4d4Tc/Ydu/657v07qMyi5UHBq8yEPU+ELAipseykhuHgHoLx5Gmfwg/w6kBrAa+pqXkLhASbLLPafoiV5AE+wJMcb6bc3zWkZs/TDBov/8sQv4alMnbKVl5wVUdsmSusl8+TFVfJiiscZAjgnWHote7lHvMk719pc6oNLRPX6VQwxPj2fiG3EoV6KujXh1P2HnmPeIqr4uQIQcFqtW2p8TGbRknbhrWGkJoeS8NHsGnKG2HCpADncg7ChrWA19QsICINErvMWudR7jUf5Iwe52yvwUpDONNSbBULBmaxXReqOG4qZF5YSpus5CnHyi7NpE0/3WDTvEDudvF+iGpxqM/xegShYbr0dImlFBomCrKvdFIBjyAKGMUgMAuNQGIgVL2f9gv5rXKpjcT7GxH2xHh6W2hUJ4/VhrCTJqxma1yyXQoZ4DAcxImwFvCamoVDaKYnONF5P59NP8P9XcvJVvRGRWIMN6hQ6rXiZAVEYLURReVMG7wavDbZKR5mp3iI74f380r7OdbHz5Lllw8slvv2CEjCkj3DqnZYa8QT1DQ+HWPWEQOUKhh0loViRUlUZl64AImJJ7opgRgqmYZejOyJuJXojcf/yyz8kohyvGEYtw3Dcok35RRlOqJwV/eFq+4etYDX1CwQxnRI7TIf6vw57pPjnO1ZltO9rAuIArZ3+7oNPOLmHETRowo3LKWQGuFDssTy8GO80b6fl+yfMsk3cH6Xww+pxNj+kh6nlyRRUKvQUB72NhmFKLotlIYRWiaQVAKeGkVQEqE6ycVMk/3hlVSUUqfu+Z6Qe5UbV1Q9TsMonURYbliaeY/EtKm2Ubl7aYWRWsBrahaIxPbots7wodZpTrUNp1tReGabd0SxutmmHbP7VJtzxLCEAl2rtC2sNaBtGyyNT7PJw6gGhpMB8xAXFzH0tEfTxjKdkiozJOzF96O3DamJAp9UKYWJURomzOzjVfAIRZUqaEVo2hCvYIIQVBBzbbglaLyy2Y8h/q2mgW4itLIuiWnGDWS9Pv/lzlMLeE3NQiCINHmg9yM8po/z8JKhbZmJt1TCZVBSsyfSchMFn44ACBLT5rQKRViiIJ5sKd3EMtj6BG80z/Ld8vcPPSYuYrGmxXHTYSmVSrij+PZLyD1MXHzuiRGsCMHGq4pUIFXBBblmQxNg7A2ZF3ZKoWVjOmIZovB3E0Naee7XE/YdMihNq7QToaUd0pkHfvepBbzmCHLtJfBRIBaw9DgV7ueedpOW5QZhmRatGK4Np0yZRmSlCj8YFCRu8BmUQMzkSCUWxZxpp/jxcV5M18g04HzJ4dnUYCRWWVojOFVc5XWPnZJ7GJSBhhGaFpZSITWV0O47iXmNXrerwieXM8PIwdU8indihESgYeF4E1omblY2r4uX602ubW6V4XI3qQW85ggR+2AgCTHP2TMPl/53Amt7rHUe5bHWGk+uGqzENDdTbUxOC1cEbsiH3u8tatTsG0Rcqx8GiZ64CDyxAseaTV7Z/jQX7TPsjsZVufjB21MkwZomqYm5JVGIYeJhtwj0C896GNCmQUtSltMWjSqVsmWVlo3hk4k3XC0MlzNhfRJ4drzBVXOFLffq7G+dSd7PsXCSJzrLLDcMS6lwqqk0zI355vsJCioBryValePfbWoBr1lALEmyQjNZ4Uz7I6zoSU6xhkVIxMTLZ3QWKnAa2NQBu2aLi/nTTIotnO/PXZrcrREayQoPysdYSuWa64upOCsxV9leI+bxh/srD7XyuAHCPhEPN/MoifHdU+E0u8llBuYNvPfc7Y25myFiEDGoRln0VV577mFUBnZDxro5T5MOHe3xSDiDIiRVCCQVxYgydFG8v79d8AKvshlepihHlH6ENU2sabCjF8jNmMb4EY4XbU60LQ0jtC107DTGfq04B5VY5EOGD/mB2aUW8JoFw2JMk6XW/ZxI38dPtD/IvR3hySVHakL1wdJrPCSvhheHJzk/PsmfapsL9hl2xq/iDkmM3jmGVmONRxsn6KUxjjvdTDOVJ06VKRHzv2/Mdd67qTj2RHzqid/MVwxEsTqetLkcTrJhWgQ/PpTrGSEKeFDwIT5/p1AEZeRL+jJk172JNU1GZoXCnyFoPAGlEjcwIdro0jjwDN/jtZ0/RjUHDIldggSsaTAqN8hMn9JOmLjHCJNlOjalTKvwilFMlaYY7TTtowKl5JSzfih1GmFNzT4sS51H+HT6c3z6eI/Hlzy9JOYpl0EYOTsr6phe5qZV7PKRbsEDbcMTy2d5ducs3zN9vpX9b0yKK6hmh/eUbgNj2iwn9/L4iqVhogi5mTYIiQFRMCamwBVBKq9bZyGUIlzrYQdiMc80Vr5/tvm0ZHzqga82LEvjNRrJEqXrg3oOPIwiMQZuTVyzV/CBGAfHk8uEwo9wxRZDucyk8UOUwUZRVQFvuVpYzo0t38hfYiN/vhJvEElpJCusth9kTR6gJEcJlGRMmLDtGhwvE0SEbhL/vpnGoohXArsFbEw8Y93B+QwOIAccbkPAReQB4O8Dp4mv2pdV9W+JyDHg94CHgHPAz6vq9t1bas0PNpZW4zQnm0/ywaUej/U8D3YyymAogmHsklk+8Ow3qg96aoTl1NG2jl7imfgGPixzTp7gqkkYTs4xv7FywZoWLXq0q07MQa+tPtQQPexp6hzsFaFMZTtUXrfOfl/Iwl7eNOxls1xTFCPQskKTGF6IYYwDeNrXIVWvl+mVlTLN5VbCLOPG4EMUX191EVQFF2JK4K4zbGXKpn+VvJzmtsvsd1Pp0Ak9VDo4PCPZwWAIKGWIqYfT8E3QvZrMIkiMxZcluQ6qNRyMkW7HA3fAf6Sq3xGRJeDbIvLHwK8A/0xVf0NEvgR8CfhP795Sa35Qmfb8+LmlX+GJ1YQfPzlmp0w5N2oxCWb2Qb3p74qAh75rxHQvo7yvW/LRFcfx9c/wws4P8U/cb1G6/px64oZO8yQ9XZkJB8T2pl6jB5qYaVxWaCTKUhpiLnj1CAEwHnZLw8AZhi56jYWPv5saZht+TVuJ4b4V9FKhTYqR9CCf+DVI5YFPrw78vpNYiqUXljneeB8b+jx5uVv19q5aCCh4bzg/Es4NMzaH30PD9LVWUEfpRzjNMRhWdQkRIQvLNElmPcZDFXe3s8KeKOq7pbA+8bwi5xhMLlK4HQ4qNPe2Aq6ql4BL1e2BiDwH3Ad8AfiJ6m6/A3yVWsBr7gLNxilW2g/y4WMJZ1qB3TJh6AyTYCiuu1K9Pjtg5qlqFPMQhKFTIOH+juI15UR4ku3Jq0zyC8ybFy4IvfQMy9oDmHmAuYdSofBaibDQDVNRF/ZHTLwKeRCuFoarOWxkgSIoZQh0EkPLGlYaQsuCEZl57rP4+hwg+17ZaRXllFQsiVos6ex+058n1Z6AACOnjLRAtdjXfbEa8BAKijBkaPucYpmGMTSDnQ2HaJiYGx4fO26So0IeYLuA7aJkwAalH6EH1IkQ3mEMXEQeAj4OfB04XYk7wGViiOVmv/NF4IvvYY01P+CstB/kMT7Fp49NsKK8OGxThOmlbNzU2x8KuJ6glT80vexVy9gbHurkdGzK+65+nOeb+Zz1/qiQhBOcZdW2EJh5gWOvFF6ZOKVhhYaBtUYsbCmCVA2Y4sbaTMBzODcsecNdJZMxQQKr2RpL0iZok16618Uw5lDvbXYapGpPe7hjdFVjdlGsipxm2QipJkiVTxMqATXVhuM0Y2TkAn0ZVuml+89MAR8yMrfDwGyQmrN0kxivSquTY8vGvYYpvnpPZV7YzDzrbNMvzuPDNNXyYLhtAReRHvCPgF9T1b7sK/FSVb3VvEtV/TLw5eox5uR8XrMYxIyTD/IZPnlsiY3c4zTm/oYqHrm/uGJ/D5D9cVzYi1fG2G/c7OuXCYnAj51YwW9+ip30NYpyi/nJTImiuaKrtBMzKxvPPOzmgZHz7IaMrjTo2IS1ZkJqppkosVnTbmkZOOH8SHhqZ8jz8i/ZLd7AhxwjKUvN+1iVe0mzhylCQtNaWnba7CkK1Xau9Jng/DS74hAtUr3Ayl4oJagyIWc3XKB0IwCsiSe1XhJomBDzxr2Qy+TmD6yOvNylby6ybp/Ea4djzaQ6OQoNy+zkNo17X57Adh54Kn+DdX2ZrNhGw8Gmpt6WgItIShTv31XV368OXxGRe1T1kojcA6zfrUXWvBW2auTPvsb8R+M8KQjWdDjT6HB/Rxk6watQhmkK3fR+1xa0BCoh31/Asu+7rzzLiY/ZGve1AydMj9R2Kcurc2Y9Q4LBVMIRY9fK2HuGoWBH+jjtob5FEZJZv+q0auIEUWzWM88Fc56r45fIy6sAWNMhT/pMTI8sOApvZ2XkiYCpBHLkAmMZ4nx2IB32buTGS6uge6EUp4FccibFVXwloNNc+IYJs5O8U8VzM+9Y98Io5YDdZJtGsKyRxGycaYEUe/nnXomed5azzssM88vR+46nlrtkhxu5nSwUAX4LeE5V/+a+H30F+GXgN6rvf3BXVlhzEwTBgiR0W/fRbZxCsOS+z+741RjjOwJViMZ2WO08zBOrlo+sDHl6p8vYy6xpUVCQyuNMZc/7no7NthI3+abtVRX2pd/BZmFpWeV007HWTOiFMzGtMMxLGCXmPg8YM3YtvFoyr4ycsunHbJstNvzLdOxxVvQkD/uzQMy4meWCYxmU8J3iHBfzp8nyy4hp0EhWua/7KVKaNLWNIZanF15JTZXNU1U6XshHbMhrVcrlwQwquB4Rg5UU1RiXd9WJLPeBy7LBBq8xmLyGSEJie9jKW7ailcdscEHxtwyRBXwYk7vABf8M/eQMMvkAx5IWyw2DaRicxnj3oAj0S89T4ftsh9fZGb2IhuJQnKfb8cD/FeDfBJ4RkaeqY/85Ubj/FxH5VeB14OfvygpvQKgicrMjizEO6r0iGNOlkSxztvujnAxneKS1Us0EjBtOYyecHynrWcmmH/N8+BeMivUqtrt4nrmRBl17Eq8wKFPyfUNm48/3POpSFQuzNqPTbAWne57rNIVuioZ4WbySOtpJQiItDjvGezMSkiiwYeoBKiUlBRMytwOAsZYynK2e495nIxBLzjdDTJ1TpvsFKU06dMMy3Wqae1pNcp9SVH1GdqRP5vscXvhk70UT2SuaKYKS+cCubDApd6p8SoM1DazILNwRqnFocO1m6A1/QwMhFOTlLgbDuj1J6Y5RhDYjF92CifMMQs6ODNj2rzMpNg9NvOH2slD+BTfvSgnwZ+7sct6KSrjFIpJgpAEQNy1CVoUPptmhiyVUt4eh17qPk80n+dXTH+ATq2M+97lXSZ9cQe4/Ee8yGuNf3mTrG8r3L5zk77zw0zybvMQr7o+qCrp58Sxvj8S2OKH3MXZwMUuZ+CjEU5RYPg7gfTUhxeyJ9bTXM0Qx2i/gpgqxdIGTrYyltE0iVRtQrou/HCKCwWJnk2eKavMyl5xch2TlNj4UuDSnDB+NbVJn8xujJz12ytboBbwfEuP7BmMSGtpmlR7H0iYNKyQSNz+nIal+GfuMbPBazJvW+flsuaCMXSyh3/FvMMk3ATCmQWLb1V4AlCEKtldBRDF66xO04kE9WblF6Uf4pmM7OckypxBn8JSMdYdx2GRcbFG4XYIennjDAlRiGtMhsT0e6P0ID4SH+aHVZc52lZNNH1ODvOGNseHFXc9rxVWed19lUmxRlBvMy5vtvRE976X2A/zS2r/Oz9xT8FM/8zLJvV1k5SQAurEb72kEe3aNk2ccPz4Z8MhXPF+/9Bj/5att3ii/RX/00gKJuNBIlrhPjlctPoXcM8s+Sau0LgGaNvaznva8cFUIIBMYOWHsYiggVu8piZFZ7nPmhTdHXUal0qDzFh7a4aAEmiTYyvOMRThKW9tV21Io3ICgJUVbY3+QECO2XoWN3LCd+5jaVsWvNWRk5TbrzZdJ5Ala/hhNm85arxYBJjk8tzvisqxzdfIShetfk3p30PhQ4EKOquJCDKPsuJwNs0Ge7eCrvG7BIJhZxkhe7Ze4INWIBVv1QL/ZSVoBIfgxIRQMNZDZq+yYc7MB0aUb4UJGCMW+UOXh6cycC7illZ5kqXUvP5x8iCdWLT99ZpdH77nK0oNVCXVfOf/qKt/aOM7TOycZ73ySy/Z5Nt1OZdx5ySh4dwiWZrrK2fST/MRpx7/6w6+T/vlPoImFixvo5i66WxUlNBPk5BJy72l0ZYkHkudZ/r/e4H+/8DCZjBhll3B+wKLYxJomK42ExMSc5+lg3qLqhRFMFG9g1ikujg2LvSrKahNu4uMgXxdicUfLgq96RXuF7TLGluePcM2m4bSSEsAS857j8YLSh9kG23STVxFGDkbOo7Mr1OhphpAxLK7Qb55iOfRYCgnexLLzzCkTH3jTvMmWf4283KkKXw7PRjE9sJzZoAyQUTBiG+czghZQiayImW1iTifpeK26ML6FBz79S9ETV5wf4vyweswolarFXF3tz62AizRIk1U+3/s3+LP3WP6dXzmPeeIM+ujj6MoytDsAJEXB48MBT1ze4BcvbvKX/vYJvrr+ef6bkDPMLlK6LQ7byO8eodt+kB9p/EV+46MFH/lzF7E/+mHC//N93KWcnddSsknKpGgCkNpAu3WVpeOXaJ4R0s8+xOpf7vKbKy/zm3/yQX4Tz7n+n+D9gHm3iWBJpEnLxtjv2Amjapr6xFeN+0XoVq5WEYRS9j6sZYhZJlfzmOp1LhtQ4khJOGbbrDYSltN4/6uF4fIk53LxDD6MmR/bKKqOddmgVd7D/WJjtaU1GC+znOcQCkQcZQh4tRRhrxd25qEI18eulRByhtkFLoSC7cbrbJRP0ixbyNgwMNsMdYvN8XMUrk8IIw7dJhrwocDrtKwdhjJg5Ddi61aN7YOnNpFZDxhTbWQKnUToZL0o9PpWTowC/pqT57Rvyt7P54M5FXBLMz3BQ93P8uNnEn767CXMh+9Hz5yARopMJpBVBp1eFh5bhXaLj33sm9inA1+98pM80/0664PvEMKEeTL67SJYVpsP8Wi3w/s/egHT7cHFTSYvZAw3mqzv9MhdMovzCUora5LlKUv9jGNnriDLLVofXeaxbxY8Et7HefN1vB+xCF540JKJj6IUMwkg98qgLGkaS8PGSkzrYzgkhVmjI6cxdNIvAlt5waZZpyRnSY+xFJqz0IsVGDphoDlFOTiwJkS3T2AkO0zCCXxozXqe5BQUjKuqvwDs9SjZX7xiZjk5Nz6uakHudglast6qNnEFcrdL5nZj2OSQPe/ZatURgiMLHuuFhr3FtpzGqxY/s0WIaYABWtbSpMntz6s8/Of9dsylgFvb5Uz3Y/zSySf5tz/7It2/+gl0bQWyHHnpHGwP0VEeT7ENi6z24PQaurpC69/7HJ/6zvN88b9r8LuvfZY/TM/NZ4Xd7SAJ79eP8cMnAq1f/ATh//4+/X9yhedePcPIJYy9RfftsENMm7o8bpNuBx7/w21OnN2i/Usf4WOnX+FTW/fxTX+impg93wKuKIUfsZ7lnGh1EGDsAjtlyWWu0vEdur7BibJN0DhCq5cojWreYebjlJWLk4w3zQUuFE/FB268n1O6VjW4irHU14ewZTYp3M6hxnlvigb67gJX7UkyH9PjGirsmi1Gbp0QMlQ9IglldfJZSf2sfWrTpqRys7CBolrg3DbO7ZIVm1WjqgDq5s8OBHzI2dIB4pfopY1ZCInZ84vetw/5LF20l/jZDMzlNKEnzRgOUbcA8vz2zKGAC93mPXyYD/MLj1ym/cEOeAdfewbdGVO+OcENFJfFF81YJV3eJjl+GbPWRN5/P3Qa/NRHzvHi4DG+FZ7gYrmDhuEhP693Shyh9fFjXc52xugLbzB8tuD1N4+xlTfJgiHzZlbQMGU6xDU1gTe2l/FBOPvtV0gSz0dXS7pXT9KXcwswzCBQuAHnm5d4Qh+ll8T85ICyJecZSIeWdrmnfJCWFVZSpWND7NUsZlYGviW7XCyfYZxvkNgWoeFnQwHODYVhGfha9hoXyqdBHfPmdSnKqFhns/0650enaFcl3kPdIvd7VwyqIW7fyXSocdwnWGnAcsPGgQjcbMhujOWqlvsagh1+bPd6VB0+ZFw2r1Hq/TSLkwQCDdO7ZuNZtcD5DFc1SZ8ONLaqrDVhNW1gTQunxbw9xXfFHAq4odc4zYO9Bg/+2ARz5gSMc8rvbVFsQn+jRV4mFFVepjWBZurpXspprQ7prF2FdpPeTxzj8e+WnFh/kMvy7UMuAH7niFgS2+LhbuB4K8Of22HrSpcLow67pZ31fY4CXlUmEvs/x0wCy0bexPaVM69cxpgmD/fGtGR51tVtvlGcH7Otb1KGR2YNm1SVcdgik10m0iHzZ1GFjg20qkwURRl7xavSl22G2UW8HwPgyMlxZD6wkytbbsK5/P8lL7fn0OsECJRuyNBd4aLd5ni5TCqGPPRxfjLbTJtmSUybUE2J0+bN22TXxJjvXKOBoAUDdxlJDKt+CS/lnhc+vZv6a5pJGVGsBKzAUqIspTFP3IdkAZyYt2fuBFwk5X36cT6yqtif/Aj6jRfIvnqJ119cY1g0GJYJum+iCFQls9ue1kXPo4Mtuo8o6ec/xCO9SzxqT/NC0luo7AuIm7it9BgfXR2y2plw8etNXtha4/ykMUuNshKnjTQlzPo9F0HIqrDKwEUv/cRzyyQ2sNqe0NNVEtvB+Z1DfX63Q9CCfn6BV3WCIabMiQjjfBPnJxiTkJlP4TSezL1KVT0YhxWMXKAgpoShDueHrI+fZVNe4FtFrHL0vmCSX5jjQifF+QH98Tmeb+U07RJGUvqTN3FhzP7rr2kRThkMoapAbBilnRia6RoTdQt4JRpRPISC/vgck2STcXOL0k9wYVKlEN7okkz7nhNiK4LVNHCsaWk3TlSpgOODfyJ3mDkTcEEkZgn0kgCTDHdxzPb5NluTFmOXUAQze6m0Ss4HSMWQecvGRg+TDkgHY0RgKbVY07zF5eP8Mp0BaERx3tLPmlUb1b0+IAZQiRt2ML0Yjul2ToUCSJ1hd9IitYHSG0oZHFI/i3eOqqNwA95sXYDBfbRtwrYOcT5+aL2aqu+z4lQQjWGDUquKzaAEdbMBs6KO0g2rbnUBkaT62bzvj3hUC7JiG2cnGElxoToxVbnLU6aVqfsHMqQGmulK3LBcjJf+pigaxdrBSNZRjamFEG5aZBTTDQ2YUPWI0VgzkJxkYjbn/lW/HeZKwAWLNS3u76a0bIn/2stceqrN9zeOM/IWF/aEaprfKeisZBZgvL1KP2vy0f/vNUblKY43DYltx0nlB9in970iGAyGjayFV8PVvMF6btkuhKUUEjRmDExtUhlAqtzXMsRKutwb1rMWQxdHSl0N36H0i+GFqRYU5QbP7/4BL5kOa51HKcOE0g+ry19DETyFhyxE4U5EmfiYPldqoNBpe08fS7D9dBLLYqFaUrptnLfEjJNKuPZhKxd84g1JlYUSNGbonEoeJ092ce7qAa/8TjFN7QPvHeOQzcJGsZT9RvH2CmNvSKtZoUagbeEBfZJRY4OsuMIiXZXfjLkS8CmjUrlaJOSv5qwPVtnI01lq2PX9LKbex7RkeuDii/r464Fzow7nhmWVHrY44g3gQ8a42OLFYYPTpcUKDJzQL+NzTUQozLVjs6YFDkWI39cngXEqbBcJF7OE82OYuK0Fi/0pGgo80M/exFcVcNO+FxNKMt+iCLF9qGGvcVXsJBfbL+xl3SyeeEemm40gNwzMjWlxqnFOZBGubbXbtHAq3MN6usIom582Ae+OgGIQdTOnZa+ohlmPpNjoCjJf9aes8uJTo5y2S7whx9mRdO6zsd6OuRNwJTDxyk5p2L3c4sqkxdXSzLZgSr12jt80zjVt6j/xAqRsr3c4P7G87tfjlIwFO9Oq5hRul9dHQhESTjYDIxc7obWswVYlz9NRWBBDJ5nfE/GtvCTzCVdLy/kxvNzPYkOjBXvTKh4NOXmxvu/DahCqfOiglMGQihL2dWOKw2dTRBqVx7rIwgUzL/QmZeAicXbjdIiv1XhlGluqwgnboS1riKQLdgK/GeEmr+S+I2KqsFq0RVrZY2qL482EbnY8ZqP4w+mueKeYKwFXPN6P+bZ7jvTq+/nRzVVeHTV4bSgcb0bjd5M413AqWr4SraGLZcMXx4HzVgh6D9/ZdLwRvltV1y0eQQv+ef88p3ePcbrVYjPPGGtB0y6x3BBWG9AySsMooepWNxBhI4Pt3PMCr1C4MS9euJfLvMj25FXycpPDbsr/zpkK136iF5pgaRjhWCOwlARaNjBylpY17BQpZ4aP0G+dZzB+ZQFi3bfLdeJd7R21rGUphWMNR9d6mtaz5OJMx/5qypWtD9Jvn2c4fm2BbfE2YlvNzkxMHBE3tUUncRTB0rIJo+UGFyYPs919jM3BM3M6C/X2mCsBjy0dHbv+Ahezh3hh0OXiRNjJA73ExCbz1T3j8KQ9b2s6rWQrjw3bRRqcc1uMy80qrWjxzrKqjkv6HCNzD5PsAcYyxomnDL2qeEWrnsfRHr7qJDcNpYx0i7HbYiiXGeaXycvNI+KJQuxMmdAkpWmFVMKsmVXDBJomltkva492eoyhuYDOQ0n43UASEtOiVY1Wa5pAaip7mHhS6yWwIh266WlG5tKRtYVIgjUNUiMkVZZWWg11SCXQNIFuoqymDZbdvWybl3A+7pEsInMm4DGetTN5jadbJdm5H61a9lhOtXu0qqKVMghOZNb3OVB1UHPKq5xnW9/kT7YukBVbC7tpNa2Uu7L7DTZMk/PpMdqNE/TsKbzeC8RQUiBeJk6ndCsx+6IMgZHbiDnQs05yi2iHm2NMk0ayylrSqgbyKrE7yN4sxKUUzjQ6DIqPMEgvkpUQFjSN7q1Ik2V6zTMcbxpWG7GpF4CvNrgTgZVUub/TZDj6IQaNC0yKo2mLxPZoN06w0jAspTrbB9Bq50CAng3c07E8sfsEG83nGeUB7/uHuu53y9wJOAScHzLILvJy59ssm3tY0eMMy24cXlr1+Z1O2vAa24yOXOx7XDCh8NM+yfPRx+G9EdBQkJc7pLZLsJ7CK5lXJl6q1qpaNXCK/T+KEPslp6ZNYlo4tll8O1yLNR1ajWM0TczxzYNUxSqBMuyNXUuMxIEFjTV8yMmPoGiltkvHnpj1gcmCwauQmkARzOwEnxroSZN2uobzkyNpi8S06djjs/7p2b5+4EElTubRqCFta2klqxRuUAv4nUMJYUwecjYGfcruBBIYlKewYrAiNG3Mwph6nRNP1akuUJLhtNjndS46Wm3ijfChIFCS+UDuDWMHwVI18Y9hk7jzHnAay4zTpEtWHC3xhkq0kmMk1W52MUsxNXFyT1VKnxroSoO2XaNIBuQHNzD8wEiTLj05PhszN3aWpg2VcAu5j611EyN0bUKLNbJk94jaokNH1oD4uRg7Q7DQrHLj8xArmK1AJzG0wxpZskO2oPu6cyjge+i0sxgl/dIBCdZYnMa4b9Dp0FWlXwZ2XYkRi5Gkahm5ACXCt4lIijUNElpkwTMqDZ3E4jWmibkQs1DGXmM/Zy1JTJPEtBGSIxRCifNAG+kSbVmr2ovC2AtOJY4bq5pZFWEv5bQpPRrJ0tGzhaQ0TJeWdsl8HGAxcKYSqaqxV4iVqa4yxpG1BZbUtGlqm8wrg1IYJPEKJKtCKUUQxj6OmfMBmvRo2O7C2mKOBbxq0kMg4KuCDUvhY+VZGqum942ZChTqQMDcdrvIxWJanVkGTx4sLhjcdT2wXYgFLI4Qx3EdwZMZYqpRYymqcUZk9MAVra7M3L59AdgrjDpytqgGHxtkZodpHnhaVenunyXqVeOukqRHzxaSIGIxmKgJIZ68ANTEnkFlmA69iBW8BoNZYFvMsYADBEIoCeoYkiFeaLkY35yOmHJByb0yDo4RGZ6ScEOhw9Eg9jkuGZOTesvYRQ88MdNxYbHlahECJa46+R09O0xRAnkITJwwKIVGNSptmpE0dnFqeaEelSNsBw3kkjEqY8fGmI0is72R2BcGRi6QBY+XcuFqAW4XVU8uE0ZOSQulbeMgh4bRWa+goatCrt5Tklfl+IvJHAu4ourxoaDQMbnkpJowdgneGBo2ehM+xJhvoY5SSkod40NevUGPyod2zxaljikoydUx8Wn0tFRm48IyH8jVUUpxxG2RUzAm857MG0YOSsOegPu4H5CHQI6j0DFlmBw5W6COIozI7IiR86QGBtZcs080zdDKnJJpSc7wSNpCtaTwI7JkyMg5GiZlWEYbNIzspdc6mHglC45SJrjZbMvFs8UcCziAx/kJRRgyScZYTRj5FK8Wr7EXdqmBiXdMoqxRhgnOZyzapdDbE21RhgmFzci1ycg5VC1lNdvRB2XiPTmOkuKI2yIjC7tMKEmdoV9EzzMxEqvwAozK6XsjowjDI2kLxVG6EeNkm6ErMUg1di3mhPsqxDQolbGPV6lH1RbgKf2ISdhmREFaGppWyAM0q83uUmM188R5xuRkYZfSLcaEqpsx5wIOLoyZFJsM7AZqAmlIKX1CFixeA57AhIK+2WHENuNik9ztHvay7wpTW+x2NggEkmAoQkIipkqpDAzJ6MsuA9k60rYo3C6jIuFKa53cr0Heo2kMDSv4KsbZdwVbssu2rDPOjq4tcrdLvzjP5cZZJm4VP+nQspUtqk3eflmyobvsmI0jbYvSDRkV61xuXSb3xwlZj5aN1boQQ64D59gMQ7aPgC3mXsBDKCj9mCzsItbQki6lpqSaogQcnkzGTOgzCduUfkQIi1sa+1ZMbTEJ2xhjaWkHR5OGWgJKiWMkY0ayQxb6R9sWWlC4AQPdAANNn9IKKamPPUFKDfSZ0JdtJrpN4QZH1xYhoygH7DSvEAik3pKHhKaP74siBHZ1zMBsMz7qtqjeF0PdQoyh6VOKkNIy1WckBAaa0zc7DHWL0i32Z0RUDy5tRkSUWw5ZfcvfJE2O00iW6DZO0TIrpNKOG5bqyMIu43KTvNyldNss6uXQ7bFni17zDC2zQoMOAY+rLgnH5SZZsf0DYYtmeppmusJy437askJLu3hxlOSMdJNhcYVJsfUDZ4uOrNHWLm5miy0GxcUfKFs00iVWGmdntggESskZ6haDMtqiKLdYDFv4b6vqJ68/+p48cBH5GeBvEVX5f1TV33gvj3drFB/G5KWLO+5JH2uas4bupRtT+KlXsXgbEe+Ma22RJTukph3zTdRRlIMfKFuUPg5oCOoYJ0s0TI8QSpzmFG5AXu5Wzcx+0GyxNbNF0JLM7f5A2kI1zGwRm+XlZG6Xwg1wfsii2+Jde+AiYoEXgZ8EzgPfBP6Kqn7/LX7nXXrg1z9OA5EGsT+yOwLtMd89tS32EGlgTKuyQ22L2hYRkRbGNKrCQLeg3QfvvAf+aeBlVX0VQET+IfAF4JYCfqeIE7Snlz2LfQZ9r9S22EO1xPvaFlDbYj+qOd5Pc72Pli3ealT123Ef8Oa+f5+vjl2DiHxRRL4lIt96D3/rOqYVU4tX+nrnqW2xR22LPWpb7HF0bXHXs1BU9cvAlwFEZAP8CNi823/3DnCC+V/nIqwR6nXeaep13lkWYZ0P3uzgexHwC8AD+/59f3XslqjqSRH51s1iOfPGIqxzEdYI9TrvNPU67yyLss6b8V5CKN8EHhORhyXuov0C8JU7s6yampqamrfjXXvgqupE5K8Bf0hMLfltVf3eHVtZTU1NTc1b8p5i4Kr6T4F/+g5/7cvv5W8eIIuwzkVYI9TrvNPU67yzLMo6b+BAKzFrampqau4c7yUGXlNTU1NziNQCXlNTU7OgHJiAi8jPiMgLIvKyiHzpoP7u2yEiD4jIn4jI90XkeyLyH1TH/7qIXBCRp6qvn52DtZ4TkWeq9XyrOnZMRP5YRF6qvq8d8hqf2Gezp0SkLyK/Ng/2FJHfFpF1EXl237Gb2k8if7t6v35XRD5xiGv8r0Xk+Wod/1hEVqvjD4nIZJ9N/+5BrPEt1nnL11hE/rPKli+IyE8f8jp/b98az4nIU9XxQ7Pnu0ZV7/oXMUvlFeARoAE8DXzgIP72baztHuAT1e0lYn+XDwB/HfiPD3t91631HHDiumP/FfCl6vaXgL9x2Ou87nW/TCxCOHR7Aj8GfAJ49u3sB/ws8H8AAnwG+PohrvGngKS6/Tf2rfGh/febA1ve9DWuPk9PA03g4UoL7GGt87qf/7fAf3HY9ny3Xwflgc/6pmjsqjPtm3LoqOolVf1OdXsAPMdNWgLMMV8Afqe6/TvAXzy8pdzAnwFeUdXXD3shAKr6z4Gr1x2+lf2+APx9jXwNWBWRew5jjar6R6rqqn9+jVg0d6jcwpa34gvAP1TVXFVfA14masJd563WKSIC/DzwPx/EWu4GByXgt9U35bARkYeAjwNfrw79teqy9bcPOzRRocAfici3ReSL1bHTqnqpun0ZOH04S7spv8C1H455syfc2n7z+p79q8QrgykPi8i/FJE/FZHPHdai9nGz13hebfk54IqqvrTv2LzZ8y2pNzErRKQH/CPg11S1D/wm8CjwMeAS8VLrsPmsqn4C+Dzw74vIj+3/ocbrwLnIC62qc/8C8L9Wh+bRntcwT/a7GSLy64ADfrc6dAk4q6ofB/5D4H8SkeXDWh8L8Bpfx1/hWgdj3uz5thyUgL/jvikHiYikRPH+XVX9fQBVvaKqXlUD8D9wQJd8b4WqXqi+rwP/mLimK9NL++r7+uGt8Bo+D3xHVa/AfNqz4lb2m6v3rIj8CvDngV+sTjRUIYmt6va3ibHlxw9rjW/xGs+VLQFEJAH+NeD3psfmzZ63w0EJ+Nz2TaniYL8FPKeqf3Pf8f3xzr8EPHv97x4kItIVkaXpbeLG1rNEO/5ydbdfBv7gcFZ4A9d4N/Nmz33cyn5fAf6tKhvlM8DuvlDLgSJx8tV/AvwFVR3vO35S4mAVROQR4DHg1cNYY7WGW73GXwF+QUSaIvIwcZ3fOOj1XcefBZ5X1fPTA/Nmz9vioHZLibv6LxLPar9+2Lu3+9b1WeJl83eBp6qvnwX+AfBMdfwrwD2HvM5HiDv5TwPfm9oQOA78M+Al4P8Ejs2BTbvAFrCy79ih25N4QrkElMQ47K/eyn7E7JP/vnq/PgN88hDX+DIxhjx9f/7d6r5/uXovPAV8B/i5Q7blLV9j4NcrW74AfP4w11kd/3vAv3vdfQ/Nnu/2qy6lr6mpqVlQ6k3MmpqamgWlFvCampqaBaUW8JqampoFpRbwmpqamgWlFvCampqaBaUW8JqampoFpRbwmpqamgXl/wdDFg1OSqcuSQAAAABJRU5ErkJggg==\n", - "text/plain": [ - "
" + "
" ] }, "metadata": { @@ -341,14 +338,17 @@ } ], "source": [ - "frames = [smoke.values.numpy('y,x')]\n", - "for frame in range(20):\n", - " print('Computing frame %d' % frame)\n", + "steps = [smoke.values.numpy('y,x')]\n", + "for time_step in range(20):\n", + " if time_step<3 or time_step%10==0: \n", + " print('Computing time step %d' % time_step)\n", " velocity, smoke, pressure = step(velocity, smoke, pressure)\n", - " if frame%5==0:\n", - " frames.append(smoke.values.numpy('y,x'))\n", + " if time_step%5==0:\n", + " steps.append(smoke.values.numpy('y,x'))\n", "\n", - "pylab.imshow(np.concatenate(frames,axis=1), origin='lower', cmap='magma')" + "fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))\n", + "for i in range(len(steps)):\n", + " axes[i].imshow(steps[i], origin='lower', cmap='magma')" ] }, { diff --git a/supervised-airfoils.ipynb b/supervised-airfoils.ipynb index 3ee9edc..4fc9a73 100644 --- a/supervised-airfoils.ipynb +++ b/supervised-airfoils.ipynb @@ -135,9 +135,9 @@ " c = np.concatenate(c,axis=1)\n", " display(Image.fromarray( cm.magma(c, bytes=True) ))\n", "\n", - "num=72\n", + "NUM=72\n", "print(\"\\nHere are all 3 inputs are shown at the top (mask,in x, in y) \\nSide by side with the 3 output channels (p,vx,vy) at the bottom:\")\n", - "showSbs(npfile[\"inputs\"][num],npfile[\"targets\"][num])\n" + "showSbs(npfile[\"inputs\"][NUM],npfile[\"targets\"][NUM])\n" ] }, { @@ -173,14 +173,14 @@ } ], "source": [ - "# some global training parameters\n", + "# some global training constants\n", "\n", "# number of training epochs\n", - "epochs = 100\n", + "EPOCHS = 100\n", "# batch size\n", - "batch_size = 10\n", + "BATCH_SIZE = 10\n", "# learning rate\n", - "lr = 0.00002\n", + "LR = 0.00002\n", "\n", "class DfpDataset():\n", " def __init__(self, inputs,targets): \n", @@ -196,8 +196,8 @@ "tdata = DfpDataset(npfile[\"inputs\"],npfile[\"targets\"])\n", "vdata = DfpDataset(npfile[\"vinputs\"],npfile[\"vtargets\"])\n", "\n", - "trainLoader = torch.utils.data.DataLoader(tdata, batch_size=batch_size, shuffle=True , drop_last=True) \n", - "valiLoader = torch.utils.data.DataLoader(vdata, batch_size=batch_size, shuffle=False, drop_last=True) \n", + "trainLoader = torch.utils.data.DataLoader(tdata, batch_size=BATCH_SIZE, shuffle=True , drop_last=True) \n", + "valiLoader = torch.utils.data.DataLoader(vdata, batch_size=BATCH_SIZE, shuffle=False, drop_last=True) \n", "\n", "print(\"Training & validation batches: {} , {}\".format(len(trainLoader),len(valiLoader) ))" ] @@ -326,9 +326,9 @@ "id": "QAl3VgKVQSI3" }, "source": [ - "Initialize net...\n", + "Next, we can initialize an instance of the `DfpNet`.\n", "\n", - "The `expo` parameter here controls the exponent for the feature maps of our Unet: this directly scales the network size (3 gives a model with ca. 150k parameters). This is relatively small for a generative model for $3 \\times 128^2 = \\text{ca. }49k$ outputs, but yields fast training times and prevents overfitting given the relatively small data set we're using here. Hence it's a good starting point." + "Here, the `EXPO` parameter here controls the exponent for the feature maps of our Unet: this directly scales the network size (3 gives a model with ca. 150k parameters). This is relatively small for a generative model for $3 \\times 128^2 = \\text{ca. }49k$ outputs, but yields fast training times and prevents overfitting given the relatively small data set we're using here. Hence it's a good starting point." ] }, { @@ -426,10 +426,10 @@ ], "source": [ "# channel exponent to control network size\n", - "expo = 3\n", + "EXPO = 3\n", "\n", "# setup network\n", - "net = DfpNet(channelExponent=expo)\n", + "net = DfpNet(channelExponent=EXPO)\n", "#print(net) # to double check the details...\n", "\n", "model_parameters = filter(lambda p: p.requires_grad, net.parameters())\n", @@ -441,10 +441,10 @@ "net.apply(weights_init)\n", "\n", "criterionL1 = nn.L1Loss()\n", - "optimizerG = optim.Adam(net.parameters(), lr=lr, betas=(0.5, 0.999), weight_decay=0.0)\n", + "optimizerG = optim.Adam(net.parameters(), lr=LR, betas=(0.5, 0.999), weight_decay=0.0)\n", "\n", - "targets = torch.autograd.Variable(torch.FloatTensor(batch_size, 3, 128, 128))\n", - "inputs = torch.autograd.Variable(torch.FloatTensor(batch_size, 3, 128, 128))\n" + "targets = torch.autograd.Variable(torch.FloatTensor(BATCH_SIZE, 3, 128, 128))\n", + "inputs = torch.autograd.Variable(torch.FloatTensor(BATCH_SIZE, 3, 128, 128))\n" ] }, { @@ -455,7 +455,7 @@ "source": [ "## Training\n", "\n", - "Finally, we can train" + "Finally, we can train the model. This step can take a while, as we'll go over all 320 samples 100 times, and continually evaluate the validation samples to keep track of how well we're doing." ] }, { @@ -584,11 +584,11 @@ "\n", "if os.path.isfile(\"model\"):\n", " print(\"Found existing model, loading & skipping training\")\n", - " net.load_state_dict(torch.load(doLoad)) # optionally, load existing model\n", + " net.load_state_dict(torch.load(\"model\")) # optionally, load existing model\n", "\n", "else:\n", " print(\"Training from scratch\")\n", - " for epoch in range(epochs):\n", + " for epoch in range(EPOCHS):\n", " net.train()\n", " L1_accum = 0.0\n", " for i, traindata in enumerate(trainLoader, 0):\n", @@ -622,7 +622,7 @@ " history_L1.append( L1_accum / len(trainLoader) )\n", " history_L1val.append( L1val_accum / len(valiLoader) )\n", "\n", - " if i<3 or i%20==0:\n", + " if epoch<3 or epoch%20==0:\n", " print( \"Epoch: {}, L1 train: {:7.5f}, L1 vali: {:7.5f}\".format(epoch, history_L1[-1], history_L1val[-1]) )\n", "\n", " torch.save(net.state_dict(), \"model\" )\n", @@ -635,11 +635,9 @@ "id": "4KuUpJsSL3Jv" }, "source": [ - "Yay, the model is trained...! \n", + "Yay, the model is trained...! The losses nicely went down in terms of absolute values: With the standard settings from an initial value of around 0.2 for the validation loss, to ca. 0.02 after 100 epochs. \n", "\n", - "The losses nicely went down in terms of absolute values, let's look at the graphs.\n", - "\n", - "This is typically important to identify longer-term trends in the data. In practice it's tricky to spot whether the overall trend of 100 or so noisy numbers is going slightly up or down." + "Let's look at the graphs to get some intuition for how the trained progressed over time. This is typically important to identify longer-term trends in the training. In practice it's tricky to spot whether the overall trend of 100 or so noisy numbers in a command line log is going slightly up or down - this is much easier to spot in a visualization." ] }, { @@ -912,11 +910,11 @@ "\n", "## Nex steps\n", "\n", - "* Experiment with learning rate, dropout, and model size to reduce the error on the test set. How low can you get it with the given training data?\n", + "* Experiment with learning rate, dropout, and model size to reduce the error on the test set. How small can you make it with the given training data?\n", "\n", "* As you'll see, it's a bit limited here what you can get out of this dataset, head over to [the main github repo of this project](https://github.com/thunil/Deep-Flow-Prediction) to download larger data sets, or generate own data\n", "\n", - "**TODO us: provide data with \"errors\" (nan & huge neg number in 1 cell), filter out to make model train...**\n" + "**(TODO, us: for exercise, provide data with \"errors\" (nan & huge neg number in 1 cell), filter out to make model train...)**\n" ] }, { @@ -947,7 +945,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.6" + "version": "3.8.5" } }, "nbformat": 4,