{ "cells": [ { "cell_type": "markdown", "id": "f683dd8a", "metadata": {}, "source": [ "# Mean Curvature Flow (revisited)\n", "\n", "[As discussed before](mcf_nb.ipynb)\n", "we can simulate the shrinking of a sphere under mean curvature flow\n", "using a finite element approach based on\n", "the following time discrete approximation:\n", "\\begin{align}\n", "\\int_{\\Gamma^n} \\big( U^{n+1} - {\\rm id}\\big) \\cdot \\varphi +\n", "\\tau \\int_{\\Gamma^n} \\big(\n", "\\theta\\nabla_{\\Gamma^n} U^{n+1} + (1-\\theta) I \\big)\n", "\\colon\\nabla_{\\Gamma^n}\\varphi\n", "=0~.\n", "\\end{align}\n", "Here $U^n$ parametrizes $\\Gamma(t^{n+1})$ over\n", "$\\Gamma^n:=\\Gamma(t^{n})$,\n", "$I$ is the identity matrix, $\\tau$ is the time step and\n", "$\\theta\\in[0,1]$ is a discretization parameter.\n", "\n", "If the initial surface $\\Gamma^0$ is a sphere of radius $R_0$,\n", "the surface remains sphere and we have an exact formula for the evolution\n", "of the radius of the surface\n", "\n", "$$R(t) = \\sqrt{R_0^2 - 4t}.$$\n", "\n", "To compare the accuracy of the surface approximation we compute an\n", "average radius of the discrete surface in each time step $t^n$ using\n", "\n", "$$R_h^n = \\frac{ \\int_{\\Gamma^n} |x| }{ |\\Gamma^n| }.$$\n", "\n", "Computing $R_h^n$ requires a grid traversal and a number of calls to\n", "interface methods on each element. Doing this on the Python side has a\n", "potential performance impact which we investigate here by comparing a\n", "pure python implementation with a hybrid approach where computing $R_h^n$\n", "is implemented in C++ using the `dune.generator.algorithm` functionality." ] }, { "cell_type": "code", "execution_count": 1, "id": "74eee58d", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:12:10.339887Z", "iopub.status.busy": "2024-02-29T12:12:10.339594Z", "iopub.status.idle": "2024-02-29T12:12:11.379512Z", "shell.execute_reply": "2024-02-29T12:12:11.378069Z" } }, "outputs": [], "source": [ "import time, io\n", "import pickle\n", "import numpy\n", "from matplotlib import pyplot as plt\n", "from matplotlib.ticker import ScalarFormatter\n", "\n", "from ufl import *\n", "import dune.ufl\n", "from dune.generator import algorithm\n", "import dune.geometry as geometry\n", "import dune.fem as fem\n", "\n", "def plot(ct, ct2):\n", " fig, ax = plt.subplots()\n", " plt.loglog(ct[0][:], ct[1][:],'*-', markersize=15, label='hybrid')\n", " plt.loglog(ct2[0][:], ct2[1][:],'*-', markersize=15, label='python')\n", " plt.grid(True)\n", " for axis in [ax.xaxis, ax.yaxis]:\n", " axis.set_major_formatter(ScalarFormatter())\n", " axis.set_minor_formatter(ScalarFormatter())\n", " yticks = ax.yaxis.get_minor_ticks()\n", " for t in yticks:\n", " t.label1.set_visible(False)\n", " xticks = ax.xaxis.get_minor_ticks()\n", " for t in xticks:\n", " t.label1.set_visible(False)\n", " plt.yticks(numpy.append(ct[1], ct2[1]))\n", " plt.xticks(ct[0])\n", " plt.legend(loc=\"upper left\")\n", " plt.xlabel('Number of Grid Elements (log)',fontsize=20)\n", " plt.gcf().subplots_adjust(bottom=0.17, left=0.16)\n", " plt.ylabel('Runtime in s (log)',fontsize=20)\n", "\n", "# polynomial order of surface approximation\n", "order = 2\n", "\n", "# initial radius\n", "R0 = 2.\n", "\n", "# end time\n", "endTime = 0.1" ] }, { "cell_type": "markdown", "id": "345ed235", "metadata": { "lines_to_next_cell": 0 }, "source": [ "Main function for calculating the mean curvature flow of a given surface.\n", "If first argument is `True` the radius of the computed surface is\n", "computed using an algorithm implemented in C++ otherwise the computation\n", "is done in Python." ] }, { "cell_type": "code", "execution_count": 2, "id": "0ca5e669", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:12:11.384844Z", "iopub.status.busy": "2024-02-29T12:12:11.384557Z", "iopub.status.idle": "2024-02-29T12:12:11.391137Z", "shell.execute_reply": "2024-02-29T12:12:11.390152Z" } }, "outputs": [], "source": [ "def calcRadius(surface):\n", " R,vol = 0, 0\n", " for e in surface.elements:\n", " rule = geometry.quadratureRule(e.type, 4)\n", " for p in rule:\n", " geo = e.geometry\n", " weight = geo.volume * p.weight\n", " R += geo.toGlobal(p.position).two_norm * weight\n", " vol += weight\n", " return R/vol\n", "\n", "code = \"\"\"\n", "#include \n", "template< class Surface >\n", "double calcRadius( const Surface &surface ) {\n", " double R = 0, vol = 0.;\n", " for( const auto &entity : elements( surface ) ) {\n", " const auto& rule = Dune::QuadratureRules::rule(entity.type(), 4);\n", " for ( const auto &p : rule ) {\n", " const auto geo = entity.geometry();\n", " const double weight = geo.volume() * p.weight();\n", " R += geo.global(p.position()).two_norm() * weight;\n", " vol += weight;\n", " }\n", " }\n", " return R/vol;\n", "}\n", "\"\"\"\n", "switchCalcRadius = lambda use_cpp,surface: \\\n", " algorithm.load('calcRadius', io.StringIO(code), surface) \\\n", " if use_cpp else calcRadius" ] }, { "cell_type": "markdown", "id": "8f2cc4cb", "metadata": { "lines_to_next_cell": 0 }, "source": [ "Timings for a number of different grid refinements is dumped to disk" ] }, { "cell_type": "code", "execution_count": 3, "id": "c23d371b", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:12:11.395700Z", "iopub.status.busy": "2024-02-29T12:12:11.395480Z", "iopub.status.idle": "2024-02-29T12:12:11.407070Z", "shell.execute_reply": "2024-02-29T12:12:11.406106Z" } }, "outputs": [], "source": [ "from dune.fem.view import geometryGridView as geoGridView\n", "from dune.fem.space import lagrange as solutionSpace\n", "from dune.fem.scheme import galerkin as solutionScheme\n", "def calculate(use_cpp, gridView):\n", " # space on Gamma_0 to describe position of Gamma(t)\n", " space = solutionSpace(gridView, dimRange=gridView.dimWorld, order=order)\n", " u = TrialFunction(space)\n", " v = TestFunction(space)\n", " x = SpatialCoordinate(space.cell())\n", " positions = space.interpolate(x, name=\"position\")\n", "\n", " # space for discrete solution on Gamma(t)\n", " surface = geoGridView(positions)\n", " space = solutionSpace(surface, dimRange=surface.dimWorld, order=order)\n", " solution = space.interpolate(x, name=\"solution\")\n", "\n", " # set up model using theta scheme\n", " theta = 0.5 # Crank-Nicolson\n", "\n", " I = Identity(3)\n", " dt = dune.ufl.Constant(0,\"dt\")\n", "\n", " a = (inner(u - x, v) + dt * inner(theta*grad(u)\n", " + (1 - theta)*I, grad(v))) * dx\n", "\n", " scheme = solutionScheme(a == 0, space, solver=\"cg\")\n", "\n", " Rexact = lambda t: numpy.sqrt(R0*R0 - 4.*t)\n", " radius = switchCalcRadius(use_cpp,surface)\n", "\n", " dt.value = 0.02\n", "\n", " numberOfLoops = 3\n", " times = numpy.zeros(numberOfLoops)\n", " errors = numpy.zeros(numberOfLoops)\n", " totalIterations = numpy.zeros(numberOfLoops, numpy.dtype(numpy.uint32))\n", " gridSizes = numpy.zeros(numberOfLoops, numpy.dtype(numpy.uint32))\n", " for i in range(numberOfLoops):\n", " positions.interpolate(x * (R0/sqrt(dot(x,x))))\n", " solution.interpolate(x)\n", " t = 0.\n", " error = abs(radius(surface)-Rexact(t))\n", " iterations = 0\n", " start = time.time()\n", " while t < endTime:\n", " info = scheme.solve(target=solution)\n", " # move the surface\n", " positions.assign(solution)\n", " # store some information about the solution process\n", " iterations += int( info[\"linear_iterations\"] )\n", " t += scheme.model.dt\n", " error = max(error, abs(radius(surface)-Rexact(t)))\n", " print(\"time used:\", time.time() - start)\n", " times[i] = time.time() - start\n", " errors[i] = error\n", " totalIterations[i] = iterations\n", " gridSizes[i] = gridView.size(2)\n", " if i < numberOfLoops - 1:\n", " gridView.hierarchicalGrid.globalRefine(1)\n", " dt.value /= 2\n", " eocs = numpy.log(errors[0:][:numberOfLoops-1] / errors[1:]) / numpy.log(numpy.sqrt(2))\n", " try:\n", " import pandas as pd\n", " keys = {'size': gridSizes, 'error': errors, \"eoc\": numpy.insert(eocs, 0, None), 'iterations': totalIterations}\n", " table = pd.DataFrame(keys, index=range(numberOfLoops),columns=['size', 'error', 'eoc', 'iterations'])\n", " print(table)\n", " except ImportError:\n", " print(\"pandas module not found so not showing table - ignored\")\n", " pass\n", " return gridSizes, times" ] }, { "cell_type": "markdown", "id": "407e8aa0", "metadata": {}, "source": [ "Compute the mean curvature flow evolution of a spherical surface. Compare\n", "computational time of a pure Python implementation and using a C++\n", "algorithm to compute the radius of the surface for verifying the\n", "algorithm." ] }, { "cell_type": "code", "execution_count": 4, "id": "234ef6f4", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:12:11.411462Z", "iopub.status.busy": "2024-02-29T12:12:11.411226Z", "iopub.status.idle": "2024-02-29T12:26:32.496558Z", "shell.execute_reply": "2024-02-29T12:26:32.495156Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time used: 0.48349428176879883\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "time used: 2.5196492671966553\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "time used: 11.43457579612732\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " size error eoc iterations\n", "0 318 0.001060 NaN 76\n", "1 766 0.000605 1.619153 208\n", "2 1745 0.000275 2.275142 420\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "time used: 0.6419370174407959\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "time used: 3.4394569396972656\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "time used: 14.355547189712524\n", " size error eoc iterations\n", "0 318 0.001060 NaN 76\n", "1 766 0.000605 1.619153 208\n", "2 1745 0.000275 2.275142 420\n" ] } ], "source": [ "# set up reference domain Gamma_0\n", "results = []\n", "from dune.alugrid import aluConformGrid as leafGridView\n", "gridView = leafGridView(\"sphere.dgf\", dimgrid=2, dimworld=3)\n", "results += [calculate(True, gridView)]\n", "\n", "gridView = leafGridView(\"sphere.dgf\", dimgrid=2, dimworld=3)\n", "results += [calculate(False, gridView)]" ] }, { "cell_type": "markdown", "id": "3e87f8ef", "metadata": {}, "source": [ "Compare the hybrid and pure Python versions" ] }, { "cell_type": "code", "execution_count": 5, "id": "93b8dc81", "metadata": { "execution": { "iopub.execute_input": "2024-02-29T12:26:32.502881Z", "iopub.status.busy": "2024-02-29T12:26:32.502414Z", "iopub.status.idle": "2024-02-29T12:26:33.058363Z", "shell.execute_reply": "2024-02-29T12:26:33.057545Z" } }, "outputs": [ { "data": { "image/jpeg": "/9j/4AAQSkZJRgABAQEAMgAyAAD/2wBDAAgGBgcGBQgHBwcJCQgKDBQNDAsLDBkSEw8UHRofHh0aHBwgJC4nICIsIxwcKDcpLDAxNDQ0Hyc5PTgyPC4zNDL/2wBDAQkJCQwLDBgNDRgyIRwhMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjL/wAARCADQASADASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD3+iiigAopG+6eSOOorjrHVtekvNOD22ovZLdXEdxPJaqjSIPO2FlIVlACxHKr8xbHXggHZVBe3SWNhcXkisyQRNKypjJCgkgZ78VzWq3+trqNyumtdGAwKcvaFhCTJGGZR5YLEIZGC7nyRjAxg2L+5vD4CvJLmCeW4NnKr/KqMRtYbyp27cjDbeozjGaANL+07v8A6Aeof99wf/HaRdVumLAaHqGVOD88Hpn/AJ6+9Ou7yO2hlnmtZpG3bURYWlOduQCEDED3xXns/wAQbuPxBCh0KSIKhjeyb77u2Np+7nPAwMdz61tSoTqfCjGrXhStzM9B/tO7/wCgHqH/AH3B/wDHaP7Tu/8AoB6h/wB9wf8Ax2maTfR6nah306W2mCgyRSwOm0nPALKu7pzisO58Uy2N5JFNo8U0RmmjieF/mOwoigrjqzyRqCD/ABZOAKyaadmapqSujf8A7Tu/+gHqH/fcH/x2j+07v/oB6h/33B/8drn7fxTLLdmz/syyluDdywDy58BAszIPMG0lSQMjruwx46U/T/Ez319plqbG0j+0OFly3zkeQ0mUT+7kAbiTyGGO9IZu/wBp3f8A0A9Q/wC+4P8A47R/ad3/ANAPUP8AvuD/AOO1aP7qZfLhPKnITA9Pen+c/wDz7S/mv+NAFL+07v8A6Aeof99wf/HaP7Tu/wDoB6h/33B/8dq75z/8+0v5r/jR5z/8+0v5r/jQBS/tO7/6Aeof99wf/HaP7Tu/+gHqH/fcH/x2rvnP/wA+0v5r/jR5z/8APtL+a/40AUv7Tu/+gHqH/fcH/wAdo/tO7/6Aeof99wf/AB2rvnP/AM+0v5r/AI0ec/8Az7S/mv8AjQBS/tO7/wCgHqH/AH3B/wDHaRNVunGV0PUMZI+/B2OP+etXvOf/AJ9pfzX/ABqG3lfyz/o8p+d+6/3j70AQf2nd/wDQD1D/AL7g/wDjtH9p3f8A0A9Q/wC+4P8A47V3zn/59pfzX/Gjzn/59pfzX/GgCl/ad3/0A9Q/77g/+O0f2nd/9APUP++4P/jtXfOf/n2l/Nf8aPOf/n2l/Nf8aAKX9p3f/QD1D/vuD/47R/ad3/0A9Q/77g/+O1d85/8An2l/Nf8AGjzn/wCfaX81/wAaAKX9p3f/AEA9Q/77g/8AjtH9p3f/AEA9Q/77g/8AjtXfOf8A59pfzX/Gjzn/AOfaX81/xoApf2nd/wDQD1D/AL7g/wDjtH9p3f8A0A9Q/wC+4P8A47V3zn/59pfzX/Gjzn/59pfzX/GgCl/ad3/0A9Q/77g/+O0f2nd/9APUP++4P/jtXfOf/n2l/Nf8aydb1i50+FWtoUMvLNHJydvr8rcY9+Kzq1I0oOctkaUqcqk1CO7LCardOMroeoYyR9+Dscf89afb6o01+lnNp91bO8Tyq0pjIIUqCPkdufnFc74Y1fVLyR4RFHJAjmSVz94biTgcjvmt53Z/EtlujZMWdx94jn54fQms8NiI4imqkVoXicPLD1PZyepqUUjMqKWYhVAySTgAVVtdV06+kaO0v7W4kX7yxTK5H1ANdBgWicAnGfYVlWWvRXtzbwCzuommE2DJswrRSbHU4Y856YyMd61jyMVk2/hvTbWaKWFboPEZGQtezNgyHL9X5yeee9ACz69bQ39xZrFNLJBAZnZSiqANuV3MwAIDKTnAwetVdW1W0uvA95qBkEMNxZSFPOZRyyHAznBJ9ic9s1cPh/S/tJultBHctGYmuInZJWUkE5dSGJ+Uc5z+ZqLVrWGy8I6lb26bIks5sAkk/dYkknkknnJoAuPfWtvHLNLPEqb8AlwAx2jgEnGeK4O6uvD1zqr3L2lw3muZXuCf3yyDG3ac8AcjH88V6H8yO/7tmDHIII9AO59qxJdAt7q+mvG8/wC1iYOsg24TABAxnkYxXJifras8NK3f+r/gdWHWEd1iI37f1+pq2OpWuowCWCVDwCybwWXPTIBOKz7rxDoIuPJu508yB3dfNgYgNGu5ipK4JUenrWv5jf8APF/zX/GsOXwtp090bmWC5eQ3P2g7pARncrbQCcAbkQnGCdoBJGQemKklaTuznk03eKsh8Xi7TJLiOFhdRmSWWJWe3YLuSURHJxwCzDB6epB4p58WaKqbzdvtKlwfIk5QDdv+79zGTu6HB5qT+xNP3s/2KUlnaTmYkBmkWQ4G7jLqGwOOvqajXw5pSxPGNPkKvE0JBmY4jK7Sgy3C4JAA4GeMVRJqSSxpOheRFG1uScdxS/arf/nvF/32KUbmlDFCoCkckd8en0qSgCL7Vb/894v++xR9qt/+e8X/AH2KlooAi+1W/wDz3i/77FH2q3/57xf99ipaKAIvtVv/AM94v++xR9qt/wDnvF/32KlooAi+1W//AD3i/wC+xUNvc24jOZ4vvv8Axj+8at1Fbf6pv+uj/wDoRoAPtVv/AM94v++xR9qt/wDnvF/32KlooAi+1W//AD3i/wC+xR9qt/8AnvF/32KlooAi+1W//PeL/vsUfarf/nvF/wB9ipaKAIvtVv8A894v++xR9qt/+e8X/fYqWigCL7Vb/wDPeL/vsU1721jQs9zCqjkkuOKinvwspt7aM3FwOqqcKn+83b+ftTY7AySLPfSCeVTlVxiOM+w7n3PP0rJ1Lu0Nfy/r+tDRU7K89PzIDqDXx220yW8HeeQjc3+6p/mfyNUdY1LR9B00vcQS3UErATPFiRiexckg4/TtWrfapDZyLAqvcXbjMdtFy59z2Ue5wKzb3QLnXbORNVuthkGEgg5ji9+fvtjueB1ArWnQjdSq6/10X9fMzqVpWcaWj/rd/wBfIwvCnjLT9QvPINjIl67uPMjRQixbiV3HI6DA6dveupeWOXxLZeXIr4s7jO05x88NUvDvhiz0G4nmsZJVjlHlvE53AlWIDZ69M/nV+f8A5GWx/wCvO4/9DhrSr7Lm/dKyMqXteX967sreJtN0u702S81ay+3QWMUk/wBmeTEcmBn5lYhGPy8F+Bk9MmuGutQ0SW1sprj4X6jbx3DqtvPFBbQyIzfd2ssisjHoOQSSAOTXpGrxWVxot/DqRUWElvIlyWOAIypDZPbjNec6Zcx3uoaVb6rq/iSXShcRPp5vtOSCK4kU5i3yKoY8gEbgu4gdTxWZqeh6LqNrq2jWt7ZtKYJFwPOBEilTtZWzzuBBB9wav1T0vTYNJsRaWxcx+ZJKS5yS0js7E/8AAmNXKACsrxMVHhfVNzFR9lk5HrtOB+PStWszxG2zwxqh2k/6JKMAeqkUAadRQ/624/66D/0FalqKH/W3H/XQf+grQBLRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABUVt/qm/66P/6EalqK2/1Tf9dH/wDQjQBLRRRQAUUUUAFFFUGvnuWMenoshBw07f6tf/ij7D8xUSmo7lRg5bFq4uYbWPzJnCjOAOpY+gHUn2FVNt3f/f32lsf4Qf3rj3P8I+nPuKlt7BIZPPkdp7kjBlfqPZR0UfSo77VYbORbdEe4vHGUt4uWI9T2Vfc4FJQnUfvbdv8AN/18ynOFNab9/wDJf18iyiW9jbEKI4YUBJPQD1JP9azftt5q/wAum5t7Q9b2ReWH/TNT1/3jx6A06PS5r2RbjWHSUqdyWif6mM9ic/fb3PHoBWtW6UYKy/4Bg3Kbu/8AglWx0620+NlgQ73O6SVzueQ+rMeTVqiipbb1ZSSSsiK2/wBU3/XR/wD0I1Sn/wCRlsf+vO4/9Dhq7bf6pv8Aro//AKEapT/8jLY/9edx/wChw0hlPxXqmiWmlyadrc00dvqMMsBEUEkhZSu1uUU44bvXF2GqW+p6hpmn6h4vlvrSG6heGBNEmgkmkRgY/MkIK4DBScBc45wK7/V7fWJRFLo99bwSx53Q3MHmRTA46kEMpGOCDjk5B4xlKvjm6PlTPoOnxngzwGW5fHsrBAD9SR7GgDqKKbGrJGqs5dgAC7AZY+pxxTqACszxGWHhjVCq7j9klGM442nJrTrM8Rlh4Y1Qqu4/ZJRjOONpyaANOoof9bcf9dB/6CtS1FD/AK24/wCug/8AQVoAlooooAKKKKACiiigAooooAKKKKACiiigAqK2/wBU3/XR/wD0I1LUVt/qm/66P/6EaAJaKKbJIkUbSSOqIoyWY4AFGwbjqrXN7FalUO55m+5FGMs34enueKg+0XN9xZgwwHrcSLy3+4p/mePY1ZtrOG1DeWCXbl5HOWc+5rLnlP4Nu/8Al3/L1NeRQ+Pft/n/AFf0K32Se9+a+YLF2tozwf8AfP8AF9On1q4zRW0BZ2SKKNckkhVUD+QqpfarFaSrbxo9zeOMpbxctj1Y9FX3P61BHpUt5KtxrDpMyndHap/qYz2PP329z+AFawpKPvP/AIL/AK+4ynVcvdX/AAF/X3jftl5q/wAunZtrM9byRfmcf9M1P/oTcegNXrHT7bT42WBDuc7pJGO55D6sx5Jq1RVOXRbEqOt3qwoooqSgooooAitv9U3/AF0f/wBCNUp/+Rlsf+vO4/8AQ4au23+qb/ro/wD6EapT/wDIy2P/AF53H/ocNAEmt6gdI0HUdSWIzNaWstwIwcbyilsfjjFcxL40vG1eMWcFlPpUU9paXUyykuZrggAR4GCF3xk56hu2K7QgMpDAEEYINc/a3Pg6Aw6baT6FGVuA8VrC8K4mzwVQfx59BmgDoaKKzdU1OXTZrMC2WSG4nSF33sChZgoPCle/8TLnoMkgUAaVZniN1TwxqhY4H2SUfiVIFZN144sba71aAeSTYogTdcKplkLOpUjHyBSnLH3OMYzb1fUYrjwNdXwVilzYM6iJTJjdGSPujpz16UAb1RQ/624/66D/ANBWj7Sn92X/AL9N/hUMVwnmT/LL98f8sm/ur7UAW6Ki+0p/dl/79N/hR9pT+7L/AN+m/wAKAJaKi+0p/dl/79N/hR9pT+7L/wB+m/woAloqL7Sn92X/AL9N/hR9pT+7L/36b/CgCWiovtKf3Zf+/Tf4UfaU/uy/9+m/woAloqL7Sn92X/v03+FH2lP7sv8A36b/AAoAloqL7Sn92X/v03+FH2lP7sv/AH6b/CgCWorb/VN/10f/ANCNZ+r61/ZlqLhbZ5l3YYHKY9Oo5rA0fW7/AFWT7L5ciwb2eaWLO4BiTjIBx6cc+4rkq42nTqqi7uT8jqpYOpUpuqvhR1E9+qSmC3Q3Fx3RTwvux6D+ftTY7BpZFmvnE0gOVjAxGh9h3Puf0qSA21tEIoYnRB2ELdfXpyarXmuQW0n2eGKa5vCMrbxxndj1ORhR7n9a3jSc37+vl0/r1/AxlVUF7unn1/r0NCWWOCJpZZFjjQZZ2OAB6k1lfa73V/l0/da2Z63ci/O4/wCman/0JvwB60yKzN3Ktzq7NO6ndHbJE/kxHscEfO3+0fwArV+0p/dl/wC/Tf4VvpHbVmGst9ER2On22nxMlumC53SOxLPI3qzHkn61aqL7Sn92X/v03+FH2lP7sv8A36b/AAqW29WUkkrIloqL7Sn92X/v03+FH2lP7sv/AH6b/CkMloqL7Sn92X/v03+FH2lP7sv/AH6b/CgCWiovtKf3Zf8Av03+FH2lP7sv/fpv8KAC2/1Tf9dH/wDQjVKf/kZbH/rzuP8A0OGp7e4QRn5Zfvv/AMsm/vH2qs8iyeJbLaGGLO4+8hX+OH1oAfr9pLf+HNUs4A5luLSWJAjhGLMhAwx4ByetcpodrqdsdPhn+HmmWnlmNXuILmEiLBGXUbc8dcZz712N5pllqGftVusu6CS3O7P+rkxvX6Hav5Vi2vgHwpaXcV1baNbpPBIsiOpbKMCCD1+hoA6Wqt1pthfMrXllbXDKCqmaJXIB6gZFWqw9Zvri11OxSB78KSXm8q0MsXljqGKoW3HoACPXtggGwYITI8hiQu6hHYqMsozgH1HJ49zWd4gTy/CupRwqiqtnIAuMAKEPAx7dKxbq+1a1m1UG51S5iDKbc2+nbTExdgIxmNt64wS+DgD3GdPXGuh4LvTiOac2L+YXzED+7O4gYJB9j9OKANyoof8AW3H/AF0H/oK0brj/AJ5Rf9/D/wDE1DE1x5k/7qL74/5aH+6v+zQBboqLdcf88ov+/h/+Jo3XH/PKL/v4f/iaAJaKi3XH/PKL/v4f/iaN1x/zyi/7+H/4mgCWiot1x/zyi/7+H/4mjdcf88ov+/h/+JoAloqLdcf88ov+/h/+Jo3XH/PKL/v4f/iaAJaKi3XH/PKL/v4f/iaZLPJBE0sqwIijJZpSAP8Ax2k3bVglfRFiqc99iU29qnn3I6qDhU92Pb6dfaq3mahqKjy4xbW56sXIkkHtx8o9yM/SrcEL20Qiht4UQdhIevqfl5NZc0p/Dou/+X+Zryxh8Wr7f5/5EDaQlzG/26Rp5XGCw+UIO4Udh2z1PrUVjZ6dpqTXkRSCMb0kbfhAFYjJz3GMVkeMtZv9PsCdPu44LyIeY0SgSZT1bK/KPQ5GenNcn4Ih1bVdV8u8uf8AR7RvtJtZx/rGck7gCD35z2zxjNdNPAU3D2ztp95zVMfUU/Yq+v3HoP2q91f5bHdaWR63br+8kH/TNT0H+034A9av2Wn22nxGO3jxuO53YlmdvVmPJP1qTdcf88ov+/h/+Jo3XH/PKL/v4f8A4mk5dFsUo2d3uS0VFuuP+eUX/fw//E0brj/nlF/38P8A8TUlEtFRbrj/AJ5Rf9/D/wDE0brj/nlF/wB/D/8AE0AS0VFuuP8AnlF/38P/AMTRuuP+eUX/AH8P/wATQBLRUW64/wCeUX/fw/8AxNG64/55Rf8Afw//ABNAEtFRbrj/AJ5Rf9/D/wDE0brj/nlF/wB/D/8AE0AFt/qm/wCuj/8AoRqlP/yMtj/153H/AKHDU9u1x5ZxFF99/wDlof7x/wBmqzmQ+JbLzFVf9DuMbWz/ABw+woAm1mC7udC1C30+Xyb2W2kS3kzjZIVIU59jg15loOnXa3kuk2PhnUdNRtZtL7zJ4wscKRxQiU78kOzMki8Zzvya9L12SaLw9qUlvdR2kyWsrR3Ev3YWCHDt7A8n6V5/4ek0G6n0y4Fr44e6d4mEt19vaIvkYZznyymeT/Dj2oA9QooooAKzPEbqnhjVCxwPsko/EqQK06zfEP8AyLWq/wDXnN/6AaANKoof9bcf9dB/6CtS1FD/AK24/wCug/8AQVoAlooooAKKKKACiiigAopGZUQu7BVAySTgAVQ+0XGocWhMNv3uGHLf7gP8z+ANRKajp17Fxg5a9Ca5vlhkEMSGa5IyIk7D1Y/wj3P61HFYtJKtxfOJZVOUQf6uP6Dufc/pVi2tYbSMpEuMnLMTlmPqT1JqvfapDZOsCq893IMx28XLt7nso9zgUo05Td5/d0/4I5VIwXu/f1LkkiQxtJK6oijLMxwAPUmsj7Zeax8unbrazPW8dfmcf9M1P/oR49AetPj0ua9kW41hkkKndHaJzDGexP8Afb3PHoBWtW+kfNmGsvJFO10yztLdoY4VZXbdI0nztI395ieSfc1JbxpgybF3hnXdjnG88ZqxUVt/qm/66P8A+hGpbb3KSS0RLRRRSGFFFFABRRRQAUUUUAFFFFABRRRQBFbf6pv+uj/+hGqU/wDyMtj/ANedx/6HDV22/wBU3/XR/wD0I1Sn/wCRlsf+vO4/9DhoAi8RSWUmj3+n3dzHCLmxuCxkjLqIwoDsVGMgbxkZGc1xGh+I3+26dZf8LEs74ebHF5R0co8wyBt3bsAnpnHfNdj4wbXf+EYvk8O2yT6jLDJHGWuPJaPKNh0OCCwbbgHA9xWN/anijUbm0sp4dG0xRcwvNLFqxmlKq6s0YTylyWA29ejGgDtqKKKACs3xD/yLWq/9ec3/AKAa0qzfEP8AyLWq/wDXnN/6AaANKoof9bcf9dB/6CtS1FD/AK24/wCug/8AQVoAlooooAKKKCcDJ6UAFVrq9jtiqYaSZ/uRIMs3+A9zxUBu5r0lLDAj6NcsMqP90fxH36fXpVi1s4rUMU3NI/Lyucs59z/TpWXO56Q27/5d/wAjXkUNZ79v8+35kC2Ul04l1AqwBytuvKL9f7x+vHtV5mVELMwVVGSScACqd9qcFgUjYPLcSf6q3iG6R/oOw9zgD1qqum3GpOJtYKmMHKWUZzGv++f4z+g9O9awpKKu/wDgv+vuMp1XJ2X/AAF/X3iG/utWOzSv3Vt/FfOuQf8Armp+9/vH5f8Aeq7Y6db6ejCFWMkhzJK53PIfVmPX+narYAAwOBRVOWllsQo63erCiiipLCorb/VN/wBdH/8AQjUtRW3+qb/ro/8A6EaAJaKKKACiiigAooooAKKKKACiiigAooooAitv9U3/AF0f/wBCNUp/+Rlsf+vO4/8AQ4au23+qb/ro/wD6EapT/wDIy2P/AF53H/ocNAC6/bz3fhzVLa1837RNaSxxeSwV95QgbSSADnGMkD3FcVoOlx20mmrL8KLexnjaMNdo1k/ksCP3gYPvOOuQN3HrXeamIW0m8W4immgMDiSKAMZHXachQvO4jgY5z0rzNNK0m41fSX0Hw74ktruG+hkeW8FzHCsQYF9/mPg/LnAHJOO2aAPVqKKKACs3xD/yLWq/9ec3/oBrSrM8Rhj4Y1QK20/ZJTnGeNpyKANOoof9bcf9dB/6CtS1FD/rbj/roP8A0FaAJaKKoyXrzyNDYKsjA4eZv9Wh/wDZj7D8SKmU1HcqMHLYnubuG0UGQksxwiKMs59AO9Vhaz3x3X3yQ9rZTkH/AHz3+nT61PbWSW7GVmaW4YYaZ/vH2HoPYU2+1K309U80s0shxFDGN0kh9FH9eg71CpyqP3vu/wA/6t6lOpGmvd+//L+r+ha+VE7Kqj6ACshtSuNTYxaQFEOcPfSLlB/1zH8Z9/u+56UDT7rVSJNWwlv1WxjbKn/ro38R9vu/XrWuqhVCqAFAwAOgro0j5sw96XkvxKljplvYB2TdJPJzLPKd0kh9z6ew4HYVcooqW23dlJJKyCiiikMKKKKACorb/VN/10f/ANCNS1Fbf6pv+uj/APoRoAlooooAKKKKACiiigAooooAKKKKACiiigCK2/1Tf9dH/wDQjVKf/kZbH/rzuP8A0OGrtt/qm/66P/6EapT/APIy2P8A153H/ocNAGL4+1A2OnabG+oXWn2l1frDd3VoD5kcflSP8pAJGWRVyB3rkvCep6JqM1rFrviHUdTu7e/ePToblJVXCzEQyOFUB3I2nc+cZHA5Nes0UAFFFYGuaNd32o213arbSNDsKee5QxlZFc7SFP3wNhPGB69KAN+srxN5f/CL6p5n3fssmPrtOP1xWM/hW/a71J0u4I1vC/mSvulaZS4ZUZfl2BVynDHIPbAq1qVreWHw9ltDLEbiDTvKkZYiynEe04GQR9e3pQB0lUvttvA100koG2UKQOSW2LwB1J9qZqNvcTWEo+2GDaN3mRKQwx/wKuO0TTdQl1lJLqS5gL5cSEHLng4JzkZHP0rhxOKqUqkacI3v17HbhsNTqU5TnK1unc7DyrnUebgNb2x6Qg4dx/tEdB7D8T2q+kaRRqkaqiKMBVGABVG+vv7PjVp7iPc52xxpCWeRvRVDZJql9g1LV1Dak6W9selmgJ3j/powbn/dBx6k1206SXvP7/6/4Y46lVv3V939f8OTPqc9+7QaOqOAdr3jjMSHuF/vt9OB3ParVjpcNizy7nmupB+8uJTl39vYewwKmSKaKNY43hRFGFVYiAB6D5qdtuP+esX/AH7P/wAVVuWlkZqOt3uS0VFtuP8AnrF/37P/AMVRtuP+esX/AH7P/wAVUlktFRbbj/nrF/37P/xVG24/56xf9+z/APFUAS0VFtuP+esX/fs//FUbbj/nrF/37P8A8VQBLRUW24/56xf9+z/8VRtuP+esX/fs/wDxVAEtRW3+qb/ro/8A6EaNtx/z1i/79n/4qobdbjyziWL77/8ALM/3j/tUAW6Ki23H/PWL/v2f/iqNtx/z1i/79n/4qgCWiottx/z1i/79n/4qjbcf89Yv+/Z/+KoAloqLbcf89Yv+/Z/+Ko23H/PWL/v2f/iqAJaKi23H/PWL/v2f/iqNtx/z1i/79n/4qgCWiottx/z1i/79n/4qjbcf89Yv+/Z/+KoAloqLbcf89Yv+/Z/+Ko23H/PWL/v2f/iqAC2/1Tf9dH/9CNUp/wDkZbH/AK87j/0OGp7dbjyziWL77/8ALM/3j/tVWcSDxLZeYyt/odxjauP44fc0AalFFFABWde6zBY3sNq8UztIUDOgG2Pe4RC2SDyxxwD0OcVo1RvtHstSkSS5jcugwGSZ4z1BGdpGcEAjPQ8jBoAo3viuwspL9XSZhZeV5kg2KhMjlAAzMBwVOScAdM5BAdq94l54KvLyJJPLuLB5FUr8wDJkZHbGeall8N6PLczXP2COO4m2mSaEmJyVYsDuUgg5YkkcnvmodesLOLwZfWn2eM21vZOI0kG4IFQ4PPcY69aANDUdStNJs2u72UxQKQGfYzYz64Brj4/iFaXOoXFlZFRJLcBYJ5EYoV2qM7QNxOc4HHuRW9c6p4dntJ4o9Q0ZmZGUB5oyucd+elcn4d0PQNE1dLv+39LuQimOVZZI+pAO5OePT6Zropex5W579Dmre25kobdTurHSorSRriR3ubxxh7iXliPQdlX2FX6yP7X8Of8AQQ0r/v8AR/40f2v4c/6CGlf9/o/8awbbd2dCSSsjXorHOseG1GTqOlAepnj/AMaX+1/Dn/QQ0r/v9H/jSGa9FZH9r+HP+ghpX/f6P/Gj+1/Dn/QQ0r/v9H/jQBr0VjjWPDZJA1HSsjqPPj4/Wl/tfw5/0ENK/wC/0f8AjQBr0Vkf2v4c/wCghpX/AH+j/wAaQax4bOcajpRwcHE8f+NAGxRWR/a/hz/oIaV/3+j/AMaQ6x4bBAOo6Vk9B58fP60AbFRW3+qb/ro//oRrN/tfw5/0ENK/7/R/41DBrHhwRndqWlDMjgZnj/vH3oA3aKyP7X8Of9BDSv8Av9H/AI0h1jw2CAdR0rJ6Dz4+f1oA2KKyP7X8Of8AQQ0r/v8AR/40f2v4c/6CGlf9/o/8aANeisj+1/Dn/QQ0r/v9H/jR/a/hz/oIaV/3+j/xoA16KyP7X8Of9BDSv+/0f+NH9r+HP+ghpX/f6P8AxoA16KyP7X8Of9BDSv8Av9H/AI0f2v4c/wCghpX/AH+j/wAaANeisj+1/Dn/AEENK/7/AEf+NH9r+HP+ghpX/f6P/GgDStv9U3/XR/8A0I1Sn/5GWx/687j/ANDhqrb6v4d8s51HS/vv/wAt4/7x96bBeabdeJrQafc2kxWzuN4t3VsfPDjOPxoA3qKKKACiiigAooooAy9X1S40x7XyrVJknlWEs0jJtZmCrkhCAMnuR6DJIByo/EwifTGXTVD6oI5ZysrERszJGMkJj0wW2Z24GTxXQXGn2V5IklzZ288icI0sSsV+hI4rHv5vCugXej2d5b2FpLcTGLT1+yjAkyDhSFwhzj0ycUAdDRVWHUrSfU7rTo5d13axxyTR7SNqybthzjBzsboe3NWqACiis/UdXt9MvNNtpklZ9QuDbRFACFYRvJlsnphD0zzigDQooooAKKKKACisnTNfi1bUr61trO78mzkaF7xwgiaRSAyL824kZ67ccHmm6p4o0jR7yOzu7lzdunmLb28Ek8m3ONxSNWIGe5GKANiiqun6jaarZrdWU3mwsSM7SpBBwQQQCCPQipEmdruWE28qoiqwmJXY5OcgYOcjAzkAcjBPOADK8Sa9N4ftkuVs0uYTuLjzirgKjOxVdpBwqMeSOQB3rPtvE7Q3FvZHSfIdpiLsCbK27PKFHOPnLM4PHTPX0159QsX8RwaPNbGS7ezluUkZFKiMMiMuSc5JZeMYIFUdVv8Aw14euNOtruyVZm8ySzitdNknZdpUuyiJGK8suTx1oA6KisVPFFhLcaVFHFeZ1OaSGEy2zwlWRGc7lkCsBhTjjnjtzW1QAUUUUAFFFNkdYo2kc4VQWJ9AKAHUVX0++ttU0621Czk821uYlmhk2ldyMMg4OCOD3qxQAUVg3viywsrXVLporiSHTZ47eZ41UhpH28LkjO3euenpyQa3qACsHVNevLB9SSOwgma0himjDXJUzCQuoUfIcNuTAHIO4cjpW9WNOdHg8RwW82nQi+vozIlyYE/eGIqdpbruGQwB/ukjpQBQHiQrf3+l6RpsN3dWyiUwJdqhd2kxL1HAUsTk9TkYHBO/p14uoaba3qDCzxLIBnpkZx0H8qW5sLO8DC6tIJwyhWEsYbIByAc9s81NHGkUaRxoqRoAqqowFA6AD0oAdRRRQAUUUUAFFFFABXDeONJt9d8SeG9MutwiuFvFLLwyHygVYehBAIPqBXc0UAeKtc6lrNj8Qbe7t3bV7PTLK3ukjQ/vJI2nJdB3DLhwP9rFbeqa7puu+KdQm0u6S6gTwteAzR8pkunAPqO47ZGa9PooA8n0/wAMaSl/4HiFqPL1DTJPt6liRebYo2Xzf7+GOQD9OnFULwWEOlaLZ6gXGk2fiy8twgLERwKtxhTjnYBwe23OeK9nooA8dutx0DXB4bZV8LnVrbayRvJAsO1ftBRUIYxb8bgpA+/jjNQXVnaDwX4m/s3W9LubOX7Gpt9GgeGGB/OA3rmRwGYYyFI+6CRnr7TRQBS03SbDRrRrbTbSK2iLFysa43MQMsfUnAyTya8q0D7B5nhr7Ju/4TT7cP7Y6+ftw3n+d/sf3c8fd217FRQBw/w50bS9OTXZrLTrW2k/te6g3wwqp8tZPlTIH3R2HaqqahaeG/GHisa3ff2W2qGGSy1GUAIyCEJtVmBTcjBjtb+8Dg5NehUUAeTTPceI9N8PRardz39s/il4obiSMRG4t1gn2n5FUFSARkDBBNTa5Ztp954ztNFga3SLSNPVYrNdrLF5s/mBAvQ7N+Mc+lep0UAeY+GP+Ea/4WlD/wAIt5P2H+w5d32bPk7/ADos47b8Y3d+meau+Orq3s/G/hiW58QDQ4/sl+PtZaIc5t/l/eqy8/TPFeg0UAcC93b3mt+CJLXW11mMX12PtoaI7z9ml4/dgLx04Fcn4atmm1PTJNQ1vSLLxIl8Gu4vscov5DvO+NmMvzRsMgHbsAIIAxXtVFAHjU+iWR8Pajqyo6al/wAJZJEl2khWSNHv/LZVYHKqVZsgdznrWhrFlPo7eNtM8NwtbRiysLgQW6thd8kizMqqQcmNOdpBJHBzzXqtFAHjKW1hH4d8WTaXreiTwN4eu1lstItHhTdsOJHzI4Dj5h2J3c9K2/7A0/T9e0SzsrRAmq6LeR3yN832sqsJUyZ+83zNyefmIr0uigDxnT7Xwo/w20GOG/0PT54Ps7ajDcxjyZ5xCymO6ClSDklhu7r0Neh+BbqG78LQvb6dDYQJLLGkduSYXAcjfGSAdjdRwOtdHRQB5begx/DLxFBL/wAfMWtyif1Ja9V1P4o6EexFVfE/9k/b/Fv9t7/+EgyP7C+95u3yV8r7Nj+Lzd+7b368V6rFZ28FzcXMUSpNcFTM46uVGAT744/AelT0AeTa7/Zv9ra3/wAJfj+0f7Nh/svOd2fKPmeRj/lp5uc7ecbe1adiDLoHwugh/wCPkCCbI6iJbJw5+nzKPqwrrNb8O/24WWTV9TtbaSLypra1kRUlXnOSULAkHGVI4q/BpdjbSW8kNsiPbQfZoSB/q4uPlHoPlX8h6UAW6KKKAP/Z", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot(results[0],results[1])" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all" }, "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.12" } }, "nbformat": 4, "nbformat_minor": 5 }