diff options
Count operations
| -rw-r--r-- | lbm_codegen.ipynb | 559 |
1 files changed, 364 insertions, 195 deletions
diff --git a/lbm_codegen.ipynb b/lbm_codegen.ipynb index 49b1740..e1593c0 100644 --- a/lbm_codegen.ipynb +++ b/lbm_codegen.ipynb @@ -23,13 +23,23 @@ "metadata": {}, "outputs": [], "source": [ - "c = [Matrix(x) for x in [(-1, 1), ( 0, 1), ( 1, 1), (-1, 0), ( 0, 0), ( 1, 0), (-1,-1), ( 0, -1), ( 1, -1)]]" + "q = 9\n", + "d = 2" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, + "outputs": [], + "source": [ + "c = [Matrix(x) for x in [(-1, 1), ( 0, 1), ( 1, 1), (-1, 0), ( 0, 0), ( 1, 0), (-1,-1), ( 0, -1), ( 1, -1)]]" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, "outputs": [ { "data": { @@ -42,7 +52,7 @@ "⎣⎣1 ⎦ ⎣1⎦ ⎣1⎦ ⎣0 ⎦ ⎣0⎦ ⎣0⎦ ⎣-1⎦ ⎣-1⎦ ⎣-1⎦⎦" ] }, - "execution_count": 3, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -53,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -62,7 +72,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -75,7 +85,7 @@ "[1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36]" ] }, - "execution_count": 5, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -86,7 +96,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -99,7 +109,7 @@ "1" ] }, - "execution_count": 6, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -110,7 +120,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": {}, "outputs": [ { @@ -125,7 +135,7 @@ "3 " ] }, - "execution_count": 7, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } @@ -144,16 +154,16 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ - "u_x, u_y, rho, tau = symbols('u_x u_y rho tau')" + "rho, tau = symbols('rho tau')" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 10, "metadata": {}, "outputs": [ { @@ -163,19 +173,19 @@ " f_next_6, f_next_7, f_next_8], dtype=object)" ] }, - "execution_count": 9, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "f_next = symarray('f_next', 9)\n", + "f_next = symarray('f_next', q)\n", "f_next" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 11, "metadata": {}, "outputs": [ { @@ -185,45 +195,45 @@ " f_curr_6, f_curr_7, f_curr_8], dtype=object)" ] }, - "execution_count": 10, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "f_curr = symarray('f_curr', 9)\n", + "f_curr = symarray('f_curr', q)\n", "f_curr" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/latex": [ - "$$\\left[\\begin{matrix}u_{x}\\\\u_{y}\\end{matrix}\\right]$$" + "$$\\left[\\begin{matrix}u_{0}\\\\u_{1}\\end{matrix}\\right]$$" ], "text/plain": [ - "⎡uₓ ⎤\n", - "⎢ ⎥\n", - "⎣u_y⎦" + "⎡u₀⎤\n", + "⎢ ⎥\n", + "⎣u₁⎦" ] }, - "execution_count": 11, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "u = Matrix([u_x, u_y])\n", + "u = Matrix(symarray('u', d))\n", "u" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ @@ -232,14 +242,105 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 14, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAACmUAAAAaCAYAAABMp8pfAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXu4blO5wH9sti0bWxQ5nI2EitouJRV9FCpU0v1UZhedDqXL6SLdVhe3qFzilFPZOaUbXXmE0iKVREjZLmHtsN0VKmzsdf54x3y++c1vjjnHvI4x13p/z7OetdYYY8455jvHO9533EFRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVRFEVpnXnm9wCYTvxck3PNqcCdwBqt5ix81gKOB6aA5YjcPuwzQ4HQB7lsj+Trbb4zkoPqmdCH8tQ1fZFJ6HqmOjakL2WqS1Qm46hMxlGZjKMyUZpG/QmlTUIvX6BlLEbtyzh9kUmRnq3HaH/kdEf5cuEi4ATgSGC+57zYeD8is9f7zsgMok69G3+PNzSao37RB9saAqq7zaL+0sygi+/YF/8pNFRuSmgU+Rvqj7ih/kizqD/SPuorhIvKTQmNqr7C65F+wLOAV2dct3rZjLwTWGD+HpiHTgITwLss1+wArEAM9Wznp4jMzgI+i8jtqT4zFAh9kcsPgdsIs2Nd9WxIX8pTl/RJJqHqmerYKH0qU12hMhlHZTKOymQclYnSBupPhM2rkElLvwLuR+qAbzpeuz7wGNJp54tQyxdoGUui9mWcPskkT88eh+R9Aum8D2lS5jSwyCHdfOAWk/5HreZonNPMc7fo+Ll9pWjSZN1695vm/ltVvD401gXejujwX4AHgfuQCctvA1a2XBeybU2iutsP2tZbZchRwC+AmxF9vxe4HPgkUh+UYQJYXCJ9V9+xT/5TSKjc/LMR8HVgGfAw4jcfC6zjMU9t8iaGC7bebklT5G/0wR/x6YuA+iNlUH9EUF9BsaFy80ud/vm+0ravMAm8NyP8NcAmrpl8BaOGY4BkeKLgunOBv1NhBugMYytEXj/znZHA6JNcno3k9VDfGclA9UzoU3nqir7JJFQ9Ux0b0rcy1QUqk3FUJuOoTMZRmShtof5E2FyBfJ8HgCWU6/R5h0m/aztZcyLU8gVaxmLUvozTN5m46tkk/ZyUeYxJuwL4a6s5GuffkfKwUsfP7SvxpMktLfF1692rgX9gn6zYN96JyGsZ8C3gCGRCyN9N+Olkl72QbWsS1d1+0LbeKkOWAxcjen4kMrD7e0T+twIbF1y/MPH3BKMTLTYmv27s4jv2zX8KBZWbf54M3MFw4t6RwPkMT70sOxEqdDZG6oMHyJ9oUeRv9MEf8emLgPojZVB/RFBfQclC5eafOv3zfaQLX2GS7EmZKwNfY3giuZW5wDmpsAHFkzK3QByDk4seMAs4EJHXAb4z0iAR8k6DGvfom1yWAEuBOb4zkkD1bEjfypMLEfX0rI8yCU3PVMdG6WOZyiNi9tmyIiJUJllEzL76uIgIlUmaiPr6ExoR/Xwn9SfCZVfgKUgn/oBynT5nA3fj/7uGVr5Ay1iSmWhf6tJHmbjo2ST9m5T5NGRg6myGneDrtZwvpTpLkE77rMG+uvXuGsjuzxdVvD5EdgP2YVxeGyATB6aB/SzXhmhbk6ju9oc29VYZxTaweBiiIyflXLsp8p2OBdZmONFiHvBxZMee51uu7eo79tF/CoHZKreIcPpOzkHy8u5U+BdM+Jcr3jcinHeMWQn4OXADcDT5Ey2g2N8I2R9RX6RfqD8iqK+gZDFb5RYRjh2t0z+fR0Q47xjTla8wSfakTICXAh8qyujBwCGpsAHFkzKPNGleaInf1cQfA2yHrNi5FznW5IdIhw2Io3EacKeJOxNZjdEH9mO4DWr6p+9H00RUV6q+yuWTSB739J2RBHl6Nht0DPpbnlyIqKZnfZZJaHqmtkzoc5nKI2L22bIiIlQmWUTMvvq4iAiVSZqI6voTqj2JqNeQ9vVe6k/0gwHunT5rI8eunZII0/I1RMvYzLYvVemzTFz0bJL+Tco8H3gEOZbrVHPN7jnpdwZ+gHQeP4To5iXA4RXSxfXA5zKeMw/p+73aXL8UWfk/BxkI+2Mq/e7mXkcAWyP1xh3Iro+/AXZMpW+rDtqP4YT95cD1iXzXye9R2HXnjSaNS3/c0chOCj827zsNPN2keZ75/3jz3t8Ebs/JE8hgQgRcANyDHAF4ObB/RtoydGEjDjXPOMESH6JtTVJGd1311jVtU7pbVg+Sz26ybLjqbdk819VbGJbTfTPiFpq4H1iu7QNd+YPPNM85ryDdRsiuMbciRyf+AfgLIuM8P6XoO9alrv9UpRwV2Q0Xu+KbJvzOsrJzkUtXsosIo+9kM3OfmxifDLYmUn/+E1kcUpaIMN4xyXuQiVe7IPMVpsmfaFHkb4Tsj/hsR4DdH2m7HZF8tg9/pOt2BPj1R9RXcEN9hWqorxCeHYWwJmX22VeYxD4pcw7wJ2BBHJA1Y/89yAuX5UXIit+LLfHbmd9bIKuCH0Mq1r8ix6V/HVldewlyLvs3gOuAvRCHow/cBnwKuAt41Pz9KeSDX+8vW97pq1x+bX7nObpdk6dns0HHoL/lqU36LJPQ9ExtmdDnMtUWKpNxVCbjqEzGUZlkM1Ptia/3Un9i5rEXcopHsmNNy9cQLWNqX7Los0xC1LO6vAHpYD4JWVl/lQnfzpL+UOBCYHvgF8juRj8FVgNeXCFd/JzLU89ZA/glMtD4L+A48/8nkJ1I5mdcs635vQVy/Fxcb1wA7AT8BBn4Tz+7qTpoDvBt5DjszYHvI3Jdgey68vWa+b3MxIMMtn4q8XOBCXfpj9saGTxcAXwFGTS4xsRtb34vBC5FFh+cmpOn1ZHj3U5BOvK/Yf7eAPlOdY7a7MJGPGJ+P2qJD1nny+iuqz6WSduU7pbVg+SzmygbZfW2bJ7r6m3yfS/LiNvB/E5/hz7RlT+4j/mdnoiT5hbgbchufjsh3/t44JUM68osir5jXer6T1XKUZHdcLErvmnC7ywrOxe59EF20Jx+7mZ+n4u8a5IHEHv7OOA5jb9BMU3XQU9FJl4dh9hTF4r8jVD9Ed/tiOSzkjrYRTsi+Wwf/kjX7Yjk+/rwR9RXcEN9hWqor1CP2dC3O1N9hcfMddYFtdsis/3TkzUHyAzQCct1ayDKdJUlHuBb5h63A89IhK8J/M1cfxtSycbMRYS+Aodz1wNhDuKMFBmWvhFRb6ZzH+WyNvLOl/jOiKFIz2aLjkE/y5MLEdX1rK8yCUnP1JaN0tcylUfE7LNlRUSoTLKImH31cRERKpM0EdVlEqo9iahXJ/h6L/Un+sEA95W430f6JpLvquVL0DI2ZKbalzr0VSYuejZp0oTCNPadMtdEdve4G1jHhO1prvluRvr1Eb38FaKLadYrmQ6G9cCWqTTfNOEfR3ZijHkhw10r3pO65jsm/E7GB4PPMHE7Zzy7qTroSwx3rFklEb4q0tE9jeycUDW/AO8w4e9gHNf+uAewT3pYbNLc4ZinH5qw9OTLJyC7Zj7EsGyVpW0bsQoiq7zdJkKzrTFldLeMPvrQ3Sp60GTZKKu3VfJcR29Bdpa7yxJ3hLn33jnXh05buv4BZJzwi0iZngauROqnPDYCvsro7lc3IN/WtluSy3dsgjr+U5VyVGQ3XOxKCNT1O8vKzkUuXckuIoy+k/hYzv+2PCeui/+rQh4jwnhHEDtyKXAtsnAF3Ha/KvI3QvRHQmhHQLY/0kU7IvlsH/5I1+0I8OuPqK/gjvoK1VBfIQw7mmRAODtl9tlXmMS+UybAfyL1WCYHIzP50wzIn5S5hYk/N+fBS0ya9GoLkBnM0wy3ck5ynolbN+feIbENkt/FnvPRNBH1lKqvcnkQqQhCoEjPZouOQX/LUxER1fWszzIJRc/Ulo3S5zJlI2J22rI8IlQmWUTMzvo4jwiVSZqI6jIJ1Z5E1KsTfL6X+hPhM8Ct02ce0jH2/VS4li9By9iQmWpf6tBnmRTp2STybqEwjX1S5hdM/EGJsA1N2HUZ6V9g4r5W8EzXdCD1wAOMLrx/jrnedkLSUrIHGq814XtlXPNZE5fcNaDJOmhHpLPdlud4EPQtNfIL8GUTvsPYFe79cW+1xIMMQk0DL3HI097m/+9Z7vUNE1/1eL62bcQxJt1ZBelCsq0xZXS3jD760N0qetBU2aiit1XyXEdvH2/if2aJj99pQ0t8H2hL129nOPlmGjkOdv2CvGyKlOvjkIHNCcRXmYfssHY/8PyM61z83iao6j9VLUdFdsPFroRAHb+ziuxc5NKV7CLC6Ds5mfzJBoeZ+I9UyGNEGO8I8Glk16nkpIwJiidaQLG/EZo/EkI7Asb9ka7aEfGzffkjXbcjfPsj6iu4o75CNdRXCMOOJhng1j/vQkQ479i1rzBJ/qTMHYGHkd2ox/gyshI3zQDJ8ITlpjuZ+KxVGiCz1B8DbrTE34us8l01I+46pKJNciAyM/ohZLvatJH3yf6ILA72nI+DkA6/+83Pb8l2IrKYYtRYFv0sdrhnKHIBeBLSgXkXokRXIQ5pFrdiP2qna/L0rGkd2wXZtn6ZeearKuS3TUIpT3XkNEWzehaKTD6CbOt/P6JjP0W2AM8jFD3r2pZVkVWXhFKmoJrdn2Jm27KYQ5E8fckh7RQzVyYTjOfftUNtKuPamVAfT5Gd/xNrXNt3mcwBPsOwPrkJ6VBbJe8iwxTNyaRpewJh1JNNvlcVH2u2+hN12n1dM8Ct0+dlJt0bEmG+/bBQyhd0X8ZiyvgcXRGKfYHq/VRTlKuLi/QnJJmUtbtFejaJvFsoTJM9KfPpyNHNf0JkkORuZGBwrVT4esDfzT1/AryW7F0QXdPF9UB6Vf7/mWuflfVCyE4kKxg9km++uddSRnfEiTnN3HOz1LObqoPiHXlOQ3z+9M/pJj4eUCmb35hLgOVkd5i79MfdwfjJUzHzkDJxgyU+naefmP+3t6SPB+v3sMTn0ZaNiDnY5G0JMoiWR0i2Fcrrrqs+lknblO5W0YMmy0ZZva2a56p6C3LM5TRil7K4B9HrLpmiOZ+gbV0HmVyxLzKBZRn2Y3VjFib+nmC0fbkx2XVo0XeEZuRW1X+qUo6K7EZRfJ3+5CnC8TvLys7F3halmWD8HV36Dqcyrgul76RoUubhJv6QgjxNEe47PhvxFz6XSjdB/rvHFPkbIfkjIbQjINsf6aIdkXy2D3+k63YE+PVH1FeYub4ClJsLk2SK2e0rTJH9nkVjarbrQrCjaQa4fbs0U4T7jj58hUnyJ2VuYp69FYx3hG6IvGBZHjS/bdugPhMp2OdZMrQO8APE2UgyH3gyw7PaQRyF45AO74uQbdfPRraX/mv5rDdObFAu7+h5qzIuN4BbEEf7esR52B9ZDbI9xVsEHwssSIUtAl6OVOBTqbgrHPLZpVxsMlkVUfBfI2VnL8QYbYZsRZ7F6gzLt2/y9KxJHQOR05XAKci24KERgp7F5amqnJrWsxBkAuIknYR0Dq2ErEz4OVJH2+xLKHrWpS0DcXzKyqpLQqm3X0k1u993Wwb2uicOew5wAO5b/89kmYA07geJ8Mcc7zlT6+NnMdqBtzVSf6V3nstipsrkY8gEtv2RjohnII3Dh5FJI3k0KZOm7UnV9lHT37nJ96riY81Wf6JOuy9U9kU6tJM7bPn2w0IpX9B9GYPyPkdXhGJf6vRT3YBMWHRlWUF8KDJZFfgg5exuSHpWhxOR/tb3Me6PXoXUP4uACxPhdyO7fnwSeCmwj7n2POCjyCBnmXRxPRD/H7MH0ql9qSXvT0LsyQOJsEXmXuciHcpptgPuQybdJp/dVB0UTzx8vSXPMbGulc0vyPfaBrgaKZ9pXPrjzkIGorN4pnmGbfeWdJ5egAwwXGZJ/yTzu0o/eBs2IuYgpC68GtnFs6hfIzSdL6u7rvpYJm1TultFD5osG2X1tkqe6+gtDCc9Z+nZZsikYtuuPG3RpE/Qpq7H3IFs8PIHZGD0VPInBy7NibvZEl70HaEZuVX1n6qUoyK7URQ/oHp/ckh+Z1nZudrbojRV+g5D7ju5z/xe2/KstVLpbIT6jqsgkwGvQ47MrkKRvxGSPxJCOwKy/ZEu2hHJZ/vwR7puR4Bff0R9BTt99xUWUG4uTJLZ7itUHVML1Y42SajvGKqvEPteGwLXpCN/AZyQcdEAMUATlpvGW2dfZIl/l4k/ICNuPxOXtX36zibu2ETY74D/TaW7HjjC8uyu+RXjq0FiNkI6oG9HKrQ/Mdzmendk4Ck523dj5P2fnLh+Gtkx5ELEwL85JzzNvcj59VWIzDMGFa+vIpc2ZHI47hXRyibPtpXsXZOnZ03qWJppwtsp04eeFelYE3KKqK5nodY98eqyfSz5DknPurRlWRTJqmvaKlNly1OTdj+ie1sG7dQ9ayN6sxuyIqfqrlURM0MmE+b+TREx8+rjY4G/kL3S14WI/stkBeMr9L4BnFnhnaC6TJq2J6HUk23ZSRcfS/2JUeq0+9pkQPFK3DnIYMHZqXCfflhI5Qu6L2NN+RxtEIp9uYlw+qnyZALd9r+cibvdddGzSbIHx3wxzfhOmW/EbceAvNX0c5HFjt8zae8me9eVvHRxPRAl0s8zYbbBkK1M/LdT4fHOh+/MuGZNpO78ZSKsyToozvMFljxnUTa/IBOGp4GvW+7p0h+XZ3cPxL47QzpPa5q0V1ruNQcZ+LqDan59W/b0vSb+KuCJDvkIzbbW1V1XvS1K25TuVtGDpspGFb2tkuc6egvwHRO/MCMulsXhqfBXIzYuec1xiM1/guU5vuiizZEkPl5wvZLXFVH0HZuiqk9ZpRwV2Q0Xu5LEZ39yVblBedm5yKUozQTN9R1GhNF38nYT9hXLs84x8S+skM8I/++4ADf7bKu3ivyNkPyRUNoRMO6PdNWOSD7bhz/SdTsCyteFtwDvT6VbhPgnT7M8w4b6CuXok69QZi5M2/TNV0hTZ0wtwr8dzWJg0vg+vrzvvsIk+TZ5FfPMveN/kiwnfya5jduQWdZbWuLjWdBZM5nzZjlva37HqzHmmvTHpNKdCzzXKaftshIyqze9GgSkM/u3yLu8CpmN/nxGV7Bezehs30Um/sbE/wAfQI4Oux74G8N3T4fHzAFegzTQflP15WpQVS4707xMDkEG+E5DnM5lwFeRFUjTqbxtafLuslNYF+TpWVM61gd86NnzMsKSOuabUOseEEdrZewrdUPSs65smY0iWXVJm2UqPlLVpTw9gDRoQ7D7odU9JyNHbJwPfKLWm1UnJJkcjKyOuxXxZy9GHPapWm9YnlDr47lI51581GGXhCST1yIda1shq9Oehkwy6nrSTJP2JKT2kU9/VP0JwXe7rwl2AdZFVtIn8emHhVS+oPsyFoLPkUUo9uUm5Jje9A54ofmr0H3/y0a4293Q9KwKawFHIzKMj/dLswkyKL5tRlzMcmTHq58jAxjPR459S++KmJcurgeS+vyo+bEdU/ihjGtI5DVrV5xtkbozWac0WQfFAx9lBu7K5heGZds20OzSH2fbNSiZxvbeyTzFC4ls3ylCdiI6imp+fRv29MPAkYj+7o4M7BcRks43obuueluUtindraIHTZWNKnqbvI9rnuvoLYhteoTxHZlWYzhImy7vpyPl/WPI4OEHgNchdvAuy3N80XXbbEPz2/XEEleKvmMT1PEpq5SjIrvhYleS+OpPriM3KC+7MvY2L00IfYdN6mc8QWwPhpMGYtZE+lEfRN61S5p6x4eBr+U8Y1tkIta1SJlLU+RvhOKPhNSOgHF/pKt2RPLZPvyRrtsRUL4u/C3jR8h/EVkkerXlGTbUV3Cnb77CK3CfC9MmffUVYnyOqc2GeT8z3VeI51wuz4r8FvDdjPABUtgmcm58ukmzeUbcFYhA5mbExSt1sgzjYhMXb2Ucz3bfJZXuE4ggbddHOfluki3N876TEXc2Mit9Zcu132J8B4FPMDqz/6PAv5CGAw7h2wD/QJylvyFbklclovpM56pyaUMmD5mfIxAlfAsio4Mynv8Wk+93ZcQtptuyFWPTs6Z0LItp7LsTxddHOdc3jQ89s5WnJHlyciWimp6FVvck+S7SIJljiQ9Nz7qwZTZssoqvjwqub5I2y1SZ8lTW7hcR0b0tg+brngMQZzQui5P42SkzJJm8BFk9tQ3S0PwF0iB/fME72IiYWfXxaxB/cENLvAsR/ZfJyogPuAJp0E8Dny16gRwiqsmkSXsSUj3Zlp108bFmuz/h0u6Lr4/yMt4yA4pX4p6AdNCunwr35YdBeOULuitjZXyO+PooJ99NEop9aboerkOeTKD7/pcydjdPz2Imye4IX4wfPZxmdKfMY03YUTnX7GDSXJUI25bhbqRJNkfq9KWILF3TgdQDDzK+GP/P5vnp3ZIOYrhzQDruSux1yvvMNW9IhDVdX19pwl+ZcQ3IwE2y7i6b32T4/pZnQLX+uJg/mDSrZsRl5eliE/aiVNoXIvb+JrKPKV1MsS40/X0+bsIvpVzbKyTbWkV3y+ijD92togdNlo2yelslz3X0FuASE7dFImwNxB7GMs36bnsgNu0Q4H6kbKRZjH+/u2ld3wrYICP9ysBh5pq2dmTK+45NUMenrFKOiuyGi11JUtT33hZ15AblZecil6I0TfYdRoTTdxJf8+5UeDyJ5MsV8ghhvWMWEyZt1k7kMUVtDFt8nI/IIR9NEFI7ArL9kS7aEfGzffkjXbcjoHxd+H5k176Y/ZBJ+eum7ruY7tsF6itk48NXKDMXpk366CskqTumFhGmHR1Q3D/vSkSY7xgzQXu+wiT5O2XG/cU7wHgD/wZgx5yL8zgDqfz3ZNQgrIashv8j2TNBt0MMftYq2u0Qw78kFZ7uhF0pIwyGivxoXsYbJGs1K8i2ui9GZLuCbLZl/LirRYweV7MIOebpxox0WeHXmrgFyLc5FVGKJo/3dKGqXNqQSbySJd7u9nKkMj8IWSGQZA9kQPDHGXnrumzFZOlZGzrmig85+NAzW3kKhdDqnpijkcHJnbGvfgpNz7q0ZUnyZDXT9KxKeXK1+20SSt2zJbIl/85YVth0SCgygdFjbq9CVjzdiDQKvmB7gRYItT5+GyKjZbaMt0hIMnkV8B9IB9qfTZrjkIFt22q6pmnLnviuJ336o6D+hEu7z1f75RXmB4adszsxnPx1N7LLUDL9b5BjWWN8+mEQXvmCbspYWZ9jtve/+K6HwS4T8NP/Usbu5ulZEb70MMk2SL/SzcCnc9L9CXnPpwKrIzp3MDKIdwkipzuBTYGXmWveinwz13RxPXAF4zI5Atl95yxkEsftyGDkU5DdTLdktPysZvJ6Fdn1QHr3gjbq6w+a/J6B7ObzR+Sb/5t5/qrAv1fML6n/D0M69v+JyPj7iTRV+uNABhi2RnTpkYz4rDx9FBmEOBM5VnIZssPInoj8dgfuy7hXkS40/X32R8r7Y8guSwdnXDfF+IRvCMe2VtVdV32kRNqmdLeKHjRdNsrobdU8V9XbmHOQnaYuQHZon49MJvkjMknrcWS3r88Ffo8sMtib7F12fNulNuriFyM+84XIGOI9yCKqFyCLMm4n++jBJsj7jk1Qx6csW46Kvo2LXUni0vfeFnXkBuVk5yIXlzQh9B22oZ8HIm3o4xEZLkHkvytwHeJXdInv/qEkRW0MW3wffJE22hFg90fabkckn+3DH+m6HRFT1o5cDHwemUz+T+BzSLm5J3XfrtsFoL5CSL5CmbkwbdJHXyGJrzG1NvSzbP9828wEX6GIBeZ35rHoewG3ZIQPkM7kiZwbz0Uq1N+lwrc3134l45qFJu6MjLjVkM6q5LbqcxED8upU2hMRhUxzObJi0baldtMcTfZK5pcj+bbN9l7dxO+aCr+ZUQN1HeJgp7GFp/k5sj1xFSKqz3SuIpe2ZLKUcRm8CXFekqyNKPSPMu4B3ZetmCw9a1LHspjGvjuRDzn40DMXHcuTkysR1fQsxLrn88hg+tOsuQ5Tz7qwZWmKZDXT9KxMeSpr94uI6NaWQfN1T2Ty8WjiZxppvDyKlLkyxPcblLwOwpGJjV8C/+OYNk3EzKmPFyKNhpdbc+1GRP9lcjPwnlTYx7A0jByIKC+Tpu1JKPVkm3ayyMdSf2KcrHafr/bLBMPVy1k/U4m0zzJh70/dw6cfFmL5gm7KWEQ5n2O29r80XQ/XwSYT8NP/4mp3i/QsZpLsia6+9HCa4U6ZF5r/93W4bolJGy+KfwUywHkt8h7Lkbrxq8ggJyXTxfWAbVekdyPfYDlyROGJSP1wD+ODJXl1SvwuDzAsV23V189Cdj+53aS7Gxkw/QqjO/KUzW+SdyGyfcjc47BUfNn+uLp5egFwPrKjyT+Q9/0UchSpjSJdaPr7TJBv46cRvU0Tkm2tqruu+kiJtE3pbpUy14buuupt1TxDNb2NmYcsFFiG7PZ8KXIM4gLE15m05GU3ZBxhBcOB5jS+7FJMG99za6TMXYF8y0eRyeG/R+qCqqeUuJD3HZugqk8J5ctRUVl3sSsxLn3vbVJHblBOdk3YWxtV+w4jwuo72Rg4BZmkshwZBz2OeroZEdY7ppkgf/erIn8jL74Pvgg0346AfH+kzXZE0TVt+yNdtyNiytqR1ZCd5fZAdu2+huzd+LtuF4D6CjZ8+Aquc2Haps++QhNjahHh2NEJ3PvnyxARzjtmMUF7vsIk+TtlvpqcE5TmIwr5hFT4gOJJmSAzrqcZnuXeBr8DTk6FXYes0kiyAFGWz7WYF1deilQOa1nit0bk9qRE2C4m7Nnm/zWQ99kpda0tPIvzaWYr2qbIk0tbMjkNWb2d5DPA1amwd5tn7ZxxD99lqws9S2IbCPcthzRt6ZmrjhVNGPCBr7rneNw6hULVsy51rEhWM03PqpQnV7vvi67rngXmnsmf3yP2bWtkVybf+K6PQRpstyHHboaAT19wApFFeod83/iQyT2MHzPwEapPygyF0OvJuhT5WOpPjJNu94XmT9g4HJHZph09z8VnDbV8QftlrIzP4VsWSXzYlz7Uwz76X1ztbp6eJZlkfFL0fHdFAAAFGklEQVSmz7I3zejx5X3mTcj7dLkzQh/xaduLCKkeLiJk29o3VHeLaVJvnwn8Hdlt7EeM7roXo2W4HXzUv0U+pU9c+959ELLc0oTWd6g0S1Ebwxav9Xg51Bdxo0k7djFwAjLpcu+M+NlahtVXGOI6F8YXocotyQRhjqkpzVLVV4DiSZmHA8fmPfwUxlf4DxidLXuN5dp5yOzrn+Y9oCavRVZhvB3ZSvpYZMXwwlS6fZBVCRvgn3WBe5FK8OnAVshuA/FZ9xsglU+8EmZHRMaPITsWgHRwP4Z0eCexhR+JFJBNkO3PjzDPeEkD79MUeXJpQyYgq2IeQbbs3xwp6/chWzbHrI7Muj/dkm/fZasLPZuPDC4sYuhQL2L0aBnfckjTlp7llScXOfnER91zEtIY2M3cP/6Zn0oXsp51oWPgJquZpmdV6m1Xu+8LH3VPmkngS1VfoAV8yOQYZDeZTc39zkT0a6aXk7xwkNWGSxG/MDR8yGQxcirAXoiPvC9wF7LDRJ8JvZ6sgquPpf6EW7svNH/CxhJkRX0XuPhhIZcv6K6MJZkk2+fwLYskPuxLH+phH/0viym2u0V6th7jK/iT+Cx7fZuUOQd4Ykb4i5Dy+lfG2+3KKD7qXVdCqofzCN22hojqbj2a0tuFwK0MjwF+BmI70wN0WobbwUf9W+RT+sK1790XocoNwu87VJqjyN/Ii9d6fBz1RerTpB37IuKDnGOJn61lWH2FIS5zYXwSqtxiQh5TU5qjjq8AxZMyL0LajFY2AX6QETaR+EmvdE+yC/BJ3CYRVOVAZBvVh4HLzDND53nArxEH5W/AuYw6MYcgFdCtwKmIDJck4t9J9mRYW/hipMJ4GLgTOcJuzzov0BJ5cmlaJjF7AVciTsl1wMGM7vLxVKScb1LmRTqmbT0bkL198eKWntcUbehZXnkaEL6cuq57bFtfT6TSha5nXdgyV1mFRp0yVbXeDt3ud133pJkkrEmZ0L1MvoM4yMvNPc8gvB0Duq6PQY4ymQa2qJPxFulaJmsiE2WWIscO3IisWJtX5yUCIfR6siwD3Hws9Sf60+4LDRc/LPTyBd2UsSSThOdzZOHD5vahHu66/8XF7hbp2eMY7Y+csKTzQd8mZW6D9IP9GPgCsrvKr5D3uIswd38Mka7r3ZlGH2xraKju1qeu3j4esYnp4/W+hwy2Kd3go/4t8il90If+5BDlBv3oO1SaocjfUH+kHOqLNENTdmx/5Ojnp9fO0cxDfYUhRXNhfBOq3CD8MTWlGer6CpPYJ2VuAXzbJRNHoQVNURRFURRFURRFURRFURQlTd8mZW6JLMK/FRmY+RfwZ6QPeH2P+VIUJR/VXUVRFEVRfKK+SFicB5zoOxOKoiiznEnskzKPRnapLWQu8BnCmjGtKIqiKIqiKIqiKIqiKIrim75NylQURVEURVEURVH6x8rIBNhDgduAdfxmR1EUZdYzSfakzB2A15W50XroCgdFURRFURRFURRFURRFUZQkFwMnI0e0r+k5L4qiKIqiKIqiKMrMZACsAJYAz/WbFUVRlFnNm5F+wHOB12bEb91tdhRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFUfrH/wPjymG5Nz1Z1QAAAABJRU5ErkJggg==\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2oAAAAXCAYAAAB6f7yJAAAABHNCSVQICAgIfAhkiAAAERBJREFUeJztnXmUHEUdxz8QSAJJICreRyJiEiKBcIkHgRUiyCWoqE/Ul0bFg2A8nlfwWhUJCPIARUVFF8SASFAEHhIRFkWUgEISBAIKG1QIIaACShJD1j++VW96e6qvmtntnp36vDdvd6p+1V31nZrp/tWvqhoCgUAgEAgEAoFAIFBbeoDB2OvuSmuTz8dRPd9RdUVGCRcAa4EJHmXtZ3FMW2vkz3bA2cAAsBHV7dNVVqiDCNr5EXTzJ2jnR9DNj6CbP0E7P4Ju/owW7fZEdX9vSv4ODPXBBl1GPSajH+gFTsg56UTg76bMz8vVty0sNueeVsG5O408R2ovYLOx8+FCc/wZnuXbzRWoPlcBJ6H+vHOVFeoggnZ+BN38Cdr5EXTzI+jmT9DOj6CbP6NJu58BDyH/Kcm2qG29yCnNdNR6C57wdGO/GXigeD3bxkuQY7BFBefuNKwjNT0lfynwL2Abz+PfCTwJbOlZvp3MQG39ZdUV6UCCdn4E3fwJ2vkRdPMj6OZP0M6PoJs/o027V6L2nJhj10/MUfO9sZ4JLEDirQBejMJ2I8kDaHqm0+sMDGFP5Ejd68ibBswFLgGe8jj2BOQA3o6c9qo5wPxdUmktRp4IfRd6WjhGN2oXEXTzJSJo50NE0M2XiNa0C7r5E7TzI+jmz2jTbhnyWz4AjPE5QA/FI2rXAf9D4ccLTLnXZ9jPAS4D/gqsR2uhlgEne9q9zpzza4n08cBnUIRnPbAaea5jgCeQUxnn9eY4i4Bd0HTKh5FTcxOwj6Mt9tynA3ugaZ+PAf9GYc3nGbuZ5nhrTd6VKAro4i3A1cA6NAf33li9k5Sp86k0z3m1r3cZm1PM+wNT6mbbexoaDbjctHcQeAXwWvP/2abNFwJrUuoTZwv0Rb4BeBQ5ibcB81Ls83gL6W0tOiXzRGP/JkfeFJN3WSwtT5u8/HYS4f+j2Kp2ZXWD+mgXEXTzJaJztCuiS+hzo7fP1fH6UNSmHUSEPudLRDV9rlt1g9Gt3RdNmYMzbPqNDeAXUTsGVfhbwF3ASpO+R4r9icBvUFTn18AZaM7pOOANHnbxc90WS5sAXI8cmP8CZ5n3XwD60JzQuD3A7ubvNOAWY3M+ch5eDfwCmJRy7mnAjcDTwHkowncU8APgCORg2uPdAxyGnNo4Y4CLgEuBnYCfIl03A181x0pSps5/NHkgp+lLsdcNJn2uacMfHOeKt3cX4LembuciB/Fu9HmBOv+twPamnVkaboOisT8EJps6/hA5uX3kh4VdPGTa9QiwKdbOXtyRRBe2rX905O1l/sb7UJ42efl1oVXtyuoWL9PJ2gXd/Blp7Yro0gnahT7nRx2vD0Vtqib0OT+Cbv6MZu1+Z/5mBbdS6UEeXG+GzSTgHyjy8wyTdrAp9xOH/XORyL8FxjrydyhpZ/mxOWd8zZVdh/V5hq5bO5CGJ/6RxHEuNulraXY0l5i8OSnnXgPsGkufBPzTtOMh5KRYxiJHbjOK+lm+aY61CNgqlr41+jAHUZSqlTq/36S9n2YmmPqudORZbHufAF7lyO8z+Q8XrA8o8jhIs0P2bBRdW0+jf5VhDHLSk5HTotyPfhhcLEJ1PjyWlqdNXn47ifAfvYLWtCurG9RHu4igmy8RnaNdEV1Cnxvdfa5u14eiNu0gIvQ5XyKq6XPdrBuMXu22N8dZlmHTTyyiFqeHfEftDGMzP5b2ApN2j8N+f5N3XsYxy9hZ7kJi2Yjgq0z5tN0nV+N2GFaZ9MMcZU4yeUmv9y6TnozygTz0QRrTCuP8yuQ9y7zfBzluaXW2DtaxLdb5OyZtL4f9NJO3NKUO0Gjve1LyV5j8QwrW53CTdknK8c43+WlTMbOYZcr2eZR9pimbtmjVfn4viKXlaZOX304iWvtR9NXORzeoj3YRQTdfIjpHuyK6hD43uvtc3a4PRW3aQUToc75EjHyf63bdYHRr9xQK9qTRj6ej9gq0Lu0OmtdOrUNOx3aJ9B3QboKDaArc23FHSoragaJAT6Pom+VHpuzeKWX+ZOoXn4I30RxnNe6dI+32/zs6zn1fynkeQxGhrR159wCPx97bCOBiGltyxl+X0twpfOq8DK17G+ewf7Wxd0VDodHeh3FPkx2P+sRfU8q76vMLk7ans0RjMOCglPws5pmyCzzKzjVlT0rJfxTpYMnTJi/fcjwa/VmPwvTJwQQXA6auRV99BY7pq11Z3aB17fZD06IfNOc+umBdB+hu3Rai6dKPo9HGK9A0jiIM0LnaFfku5tnMR4NSj5vX73EPliUZoHN1g/b9zkFj3cg3C9Z3gPZqV6frQxGbXprbmHVzZxlwlOu2Pvd8NOj7CLopXokCAnkMUI8+V5VuA7jbeU5OfdPKdUufGwN8hca93P3m/Fs5bEEzEzeltiDhqKUdxMU5xv5jpsJxViJHbzZaZ2ZZB+yLFs8ditZuPY082s8iB6qMHcBuSKh42kFI/FtT6v58NK/1iVjabHOcpcQEibEH2gTkfse5f+Wwn4qcy8uQ8xJnIvAyGnNTbZ0h/4Hd8UcflK3zVmhU4k5gg8Pe7vI43pEHjfZehXtHx93MOdIici4N90c3Oq65w6DPCvwe+eBau1gU6zi66rUjGqWJj9AU0SYrHzQgcRZy1m4EPoQ2lZlJdvvPRGv74swGjkQXp4FE3u0Zx7L4aldWN2hduwnAcrSuscxuUN2uWw9aA3sLGuj5MnAt6m+P5dS3k7Ur8l3Ms/k72qjqXqTdPDQbYk+yp+Z0sm7Qnt850KyX4yg3jand2tXp+lDUZhVDIxLJey8X3d7nJqN7rRvRYMoj5lxrC9S3Ln2uqu/q3gwNxOyC7nd/mlPfbu9zn0aDefOQP7QrckY3IAcuyTb47bKeGVF7F8W85I9mHH8sjW3gB5Fz5ory5NmdYNIj8368eZ/2YdrnMFyUSF9g0j/oKDMJ/SBen0i35z7OUcbuUrPQkTfH5J2ZqPMNDtssytZ5V2Pv2pQEGtNWb0zJt+39QEr+8Sb/fQXrM8nYL0853hgUJXkYv+fj2QWeyc1LAF6Evjhr0IjHHQydkmnX/k1xlLU6nOxIS9MmLx/gZuB7ibR70RzqskTmfD0eZcFfu7K6xdNb0c4ySPGImouI7tQNGhH6IwraJ4noDO2K6FJWO5BzW8beEtEZusXTW+lz26NZFwegkeKiETUXEf7a1en6UMSm19SjHUR0T587maGD4q0SMfJ9rg7XB9A961/wux+L6J4+dyXNUcLzTXqSLVG70maigUdEbTu0HeX/aEwxTDIVrSna3ZFn2YhGb69FH8C+aBORZOQgz8562Taitsm80qZKfiphb7F1dUXhdkdiJj3yrJ1ksrx4ey5bB9vpyz57rmydZ5u/aU7sQ2i0Ke1B2La9aZHKLD1c9dmAOmjaZxWhiJp9rEAZtkCjHsnIKeg5f79H+h+NRtb2TdjNQH18daLsOBpfzngfKqpNWv5Y1GdOT6QvBV6TUma4aEW7srpB69rVhU7XbRL6juZF04aDkdSuiC5ltBsDvA05ujcVsG8nndjnvoum8l+HdmGugrpdH4ra7IimSW1EOzMvpDk6Mdx0Wp87Cs1MWYwG/R8Evo9mhZW9r2iFTtMtyVgUpLHLUUaSTtPuJhRAmYF2gJyJBqZcg+7TUfuKRBWb6MEdUTvTpJ+aUXYvYxPfPXB3NN0vyU5od8TV6CahqJ3ldhQyjDuZfzbnT25AMZ9GtC+Ztxw5Dq5dJj9myhyTSL89o8w1pozL+eozefH1IMtN2psd9qCOl1wLWLbONm1eyjmgsRZuJ0deVntBnX0D7jV5aRr+waTPTaQfiJ69dj8agY3Tx9AoqovpxuZiR97VaIQtaw3FMlN+WixtAo2dfgYZ2k/ztMnLt9HM/RLpX0DTXcoS4T961Yp2ZXWD1rWLU2VErZN1A61NvQ3PB2/SOdoV0aWIzSz0G7UJXZsOzbDNIqIzdIPW+9xxaLDO5vdTTUStbteHIjaHoJk6s9D18tdocPWZGcdMI6J7+tx681qE7i+PRd/b+Sn2eUSMfJ+rw/Xhbei3LrnpRlEiuqfPbYn622bkJA6SvkbuWJN/Qko+lIyozUKd+29oPUMad6ApNDvTmHu5ADkIy5AjtRZ4KfBGU+Y9qFFF7UAe8UwkWnwh3iIU7bsK3XysQY7Oy5F3O52hHvQ4U9eVaKQqiSs6Zs+9IqXMHsipXJeS9xTaNcbySVPfJSh6uAJ92C8059+aoQ/I9qmz/f+ryEn8D9I4Pt94CboYHIxC3PHzZbV3rDnmcprX5KXVB7Tm8BoUEr4EjXbtZs6/GoWv/50oY7+QWYsvk5FWyxS0Q6fdZTONa9D87BvQ4wMmIudxBbo4bktjE5k8bfLy4yRHqrZwpA03rWhXRjdor3ZV08m6nYYGCeZQbN1Luxkp7YroUlS7VWiWwmT0m3kBuglp1/S0InRSn5uOphjNSckfSep0fShqc3Xs/5UoinAfugE+I6Ou7aaT+hw0ZvLYZSi3oZv2+eRvitFOOk23JO9FffDBgvbtpNO0Oxp4JwpM/BldJ85CgYfkbvYHoWvu5SnHyqSH5ojab0ya6+neSezWlfuY90ch52kV2jxiIwrZfx85UJS0A938D6It55N8GM353IimSZ6DPtRHaf6w7XHOzWhLfPv/vDJTTJ5rY4NxyJFxPVB6bxTRWmNs1qEf5HNpjgD61Bnkta9CI0yDyGmLM9ac/+aS5/OtD2hDkevQKNeTqM1fwj0XGfRD+zjZz1Y7DXek7kjk4OXtSDYefbEeRM/uuBWFySejH4T+mG2r2oB03wS8NZF+DuXXLkJro1etaFdGN2iPdnGqjKh1qm5fR+tAk89oLEtE/bUrokvZPme5Fl2nyhJRf92g9T4XmfxNsdegOdcm3GvU87DH7ClZrk7Xh6I2Lq4Hvl2yDHRPnwMN+Ca/l+9GA9U+RIx8n6v6ujoFORNHFrR3EdE9fe5vND+n+XM0r0PbHgVt0h7LZeknZcC+B/fUx07m3ahNn6i6IjVnIdIpa41hVUxGPxhf8yx/KO5HR9SBm9H6jTj34LeZyHBQZ+0srTpqw0GddTub9jhpw0WdtUtyHXrMSh2oo26T0ayL+OsWtHZoF/w2KGg3ddQtjfEoIlDVOr8kddVuMUMf3wTaee/OCurioq66xelFfa3MzvAjQV21e5TmqYwLaXbUPozuWfIew9RPjqNmX3eXq2dljAGe40ifi6I1D6DQZyCd8WgU6oqqK+LgCBQNfJ5n+WehzRIWo2cBzkDrJoo+P2o4eTuKAL8PTWs9E/XZKVVWKkZdtZuIphbMpjEQM5uhU4WrpK66fQtFpg9A3yf7qtPvY121OwVdXKeiJQF2PcIhFdYpTl11S9JPa2vU2k2ddTsdzT55KZqpdCX6/obrQzZ7oxlKn0Vr79+KllP4rlFrN3XVzbIluh88peqKOKirdn3oES6HoWvEm9BGfV+P2WyDonyXphxjB4b6YE5HbSpDH7ictdCtTsxCN/KXo3nb30CjKYNIqDpGierIfug5dhOqrsgw8Fq0+PRJtAnAUtzOfRUcj6b6bkDz6pObi1RNHbXrofkHbZBiD9EcKeqom0uzOs6iqKN2fejmZQNaR30tWldbJ+qoW5J+6uWoQX11uxjd2G1EOz8uoX6R8LpqdxhaP78ezVJZQD0iuJa66gZaQ5XcjKNO1FG7SWigfTWa2ngfWp8bf0bxzuhaOzXlGNsy1AfrbX81q2M6esj0P9CX8r9oMd+paFv/QCAQCAQCgUAgEAgEAoFAIBAIBAKBQCAQCAQCgUAgEAgEAoFAIBAIBAKBQCAQCAwP/we+rdDNnvBn4AAAAABJRU5ErkJggg==\n", "text/latex": [ - "$$\\left ( \\left [ \\left ( m_{0}, \\quad f_{curr 3} + f_{curr 6}\\right ), \\quad \\left ( m_{1}, \\quad f_{curr 1} + f_{curr 2}\\right ), \\quad \\left ( m_{2}, \\quad f_{curr 0} + f_{curr 4} + f_{curr 5} + f_{curr 7} + f_{curr 8} + m_{0} + m_{1}\\right ), \\quad \\left ( m_{3}, \\quad \\frac{1}{m_{2}}\\right ), \\quad \\left ( m_{4}, \\quad f_{curr 0} - f_{curr 8}\\right )\\right ], \\quad \\left [ Assignment(rho, m2), \\quad Assignment(u_x, -m3*(-f_curr_2 - f_curr_5 + m0 + m4)), \\quad Assignment(u_y, m3*(-f_curr_6 - f_curr_7 + m1 + m4))\\right ]\\right )$$" + "$$\\left [ Assignment(rho, f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8)\\right ]$$" + ], + "text/plain": [ + "[ρ := f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr\n", + "_6 + f_curr_7 + f_curr_8]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "moments = [ Assignment(rho, sum(f_curr)) ]\n", + "moments" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAADhwAAAAXCAYAAABKvARRAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXuUNldVp5/kg+QLSSBIZiLMSAJiEkKCSSCio4GWcBEiQgZHlwwuChVwCMbLmhkNzJBWkYAwWYDCiDNIgzEiEAy3FQ0IH1chgLkBuTBAJ2LIDZDLkIsf6fljn5qur7quu+qtOlX1e9Z6V3fXOafq1K/3qbP3fqtOgRBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCFEhjVgK/O5ZtTe1PNbWD9/YeyOzIQ3A7cABzvapv+LZ/TaIz/3Bl4DbAJ3YX377TE7NCGknQ/p5kfa+ZBuPqSbH2nnQ7r5mYt2j8D6/stjd2QCKL7rly7x3RCsun9zuYaMgbTzId38SDsf0s3PXLTry8+sm5PnotfQSDc/0s6HdPMh3fxIOx/SzY+08zEn3er8X+Vhm6EcbL9MPQc7p2vE0Eg7H9LNh3TzI+18zEm3IXKoc9JraKSdD+nmR9r5kG4+pJsfaedDuvmZk3ZV/m+db3w4+z5LuFVUaS0U7AHWgRfUdOgQ4CuhzUU1dVfBBeHYR49w7KlR90DgI4G7Qz0P54f9H+ts3zfvxvrzXuAlmD0/dMwOTQhp50O6+ZF2PqSbD+nmR9r5kG5+5qTdXwNfxeKnmFF8Nx1WHd+tmiH6N6dryNBIOx/SzY+08yHd/MxJu6Z+5hHA97AvSrI0mZPnpNeQSDc/0s6HdPMh3fxIOx/SzY+08zE33er83ynkYZWDnQ5LyMHO7RoxJNLOh3TzId38SDsfc9Nt1TnUuek1JNLOh3TzI+18SDcf0s2PtPMh3fzMTbsq/7eq7F7Yua9jD19WPnC43rAzrwz17wZuaNimTx6IPeC23wjHnhrpA4HHlJRfAvwzcJBz/58HvgPs72zfJ8di5/o3Y3dkgkg7H9LNj7TzId18SDc/0s6HdPMzN+1+BDufF47dkRoU302HVcd3q2bV/ZvbNWRIpJ0P6eZH2vmQbn7mpl1TP/O5od5P5rbXzclz02sopJsfaedDuvmQbn6knQ/p5kfa+ZijbnX+7xTysMrBToe552DneI0YCmnnQ7r5kG5+pJ2POeq2yhzqHPUaCmnnQ7r5kXY+pJsP6eZH2vmQbn7mqF2V/9vUN95D5oFD7wNixwFnYeJeCfwA9hrFIbkBuIaSpyfFPjwCeyDwCwVlRwOPA94K3O7Y98FYkvVyLDE+No8NPy8ctRfDk2BjYa3DPpaoXYJ085Ig7TwkSDcvCd20k25+pJ0P6eZnbtpdisUtzwN2jdyXMhTfTYtVxnerZoj+ze0a0pQEXX/HYqm6JcjmxmKpuiXI5vI09TPPAL4GfDizrcmcPDe9mpKgWNBDgsaohwTp5iVBY9VDgmzOS4JszkOCbM5DgnQros7/jT0PqxzstJh7DnaO14gmJOj66iFBuo2FdPMj7XzMUbdV5lDnqFcTEjQveElQbsFDgmzOQ4J085KgseohQTbnJUE25yFBNldElf/bOX+6RvM3HH4A+BfsdZFvDu0eX1H/VOAdwBeBO4BbQodf6qz3k+GYf5jbvhv4HeyNe3cA12NPYO4Cvo0lbrM8PuznXOB44ALgZixx+HHgUQXnkh77lcDJwEXA14FvYq+Z/P5Q77iwv1tC2XuwFeOKeDpwMXAbcBeWtEz7nadNn18e6hZ9nhnqvCz8fVpJ39LzfQX2VOs7w/luAQ8Dfjz8/ppwzucDN5X0J8t+2ED/EBas3Q5cBjyrpH4dT6f8XI9tuI8XhvpnFJQdGcrekdlWp01deZ8k+C+aXbVrqxvEo12CdPOSMB3tmugim5uvzcU4PzSt0wcJsjkvCePY3FJ1g3lrd05o88QWbYZkqfHdKmI7aB7fte1vH/GdZ5z0SV3/uqA5X9ffoZHNyeaGZum6JcjmiqjzM+8D3Am8Mbe9ak6es15NSPDZmnTTGPWQIN28JIwzVkE+nGyuHTHa3BR0g2XbXIJ0K+Mcqv3fuvIxWWoONnvspdxjE3MOdu7XiDoSpnN9lQ9pTN3mvMToQzatMzayOR9z163OR2ybQ136GE3QvOAlwaddjDY3Bd1g2TaXIN28JMxnrDat0wcJsjkvCfOxuSnoBvO3uXMo93+rylL2hDqA7w2Hz8BO6HXA1cBVYfvJJfVfiK388Qjg74DzgHcDBwI/5aiXPdZlmW0HAx/EkoTfBV4d/n4xsAEckqsPcFL4eTTwqVDnTdhDcD8GvAs4tOTYRwMfBb4HvAFbDe5pwJ8BT8GSuOn+rgNOxxLHWXYBfwm8HXgI8DZM17uBPwj7ytOmz58JZWCJ0t/NfD4Utj8unMMnCo6VPd/jgY+Evr0eS8Jeg/2/wAbHp7EA7M0l/Uk5CFu5743AYaGPb8QSyRvUv6aziK+G87oV2Js5z3WKV50rIj3XzxSUPTL8zNpQnTZ15bHQVbu2umXbTFk76eZnaO2a6DIF7WRzPmKcH5rWGRvZnA/p5mfO2n0s/Ky6gWQslhzf9RnbQfv4rm1/+4zv2oyTPqnrXxc05/uZ8/V3lcjm/MjmfEg3P3PWrs7PPB04gJ1faFTNyXPWa5VINz/Szod08yEfzo9szkeMNjcF3UA252XuutX5v7HmYZecg80ee2n32MSYg537NWKV6B4HH7I5HzH6kE3rjI1szsfcdes7h6ox6mfutrYqYrS5KegGsjkv0s1HjGO1aZ2xkc35iNHmpqAbzN/mqvzfTvnTNexJxPWKOocC/4StEnbfsO2Jod1fFdQ/AvsnfARzyvMc3rJeyl+EYx6T2XZ+2PbfsbfnpZzG9hOnv57bz1vC9lvYmcy9MJSdWnLsm4CHZ7YfCnwjnMdXscRkygFYsvRubIW4lD8O+zoXuEdm+z2xf+YWtopblz4/N2x7Ljs5OPT3qoKylPR8vw38aEH5Rii/uWF/wFap22Lng4X/Cnvb4R1s21cbdmGJ8CsdbQG+jF04ijgX6/NPZ7bVaVNX3icJ/qe0oZt2bXWDeLRLkG5eEqajXRNdZHPztrnY5oemdfogQTbnJWEcm1uybjBf7e4T9nNph32sgqXHd33GdtA+vvPEo13jO8846Ysm/euK5nxdf4dGNiebG5ol65Ygmyuizs98G/bWjqzf0mROnqteTUhQLOghQWPUQ4J085KgXKGHBNmcl4T52NxUdIPl2lyCdCujzv+NMQ+79Bxs9thLuccm9hzsnK8RdSRM5/oqH9KYus11ITYfsmmdGJDN+ZizbqvIoS55jCZoXvCSoNyChwTZnIcE6eYlYT5jtWmdPkiQzXlJmI/NTUU3mLfNVfm/TfKne0KdHayFgvWKxueFOmdmtj0gbLuuoP5jQtkbKvbZpl7K1ZiY6RsafzS0v6ik/vUUJzavDdtPL2jzklCWf3rz6rA9vyIc2JOoW8AzC8reF8ruF/5+FJYcLetzmsR8dsc+/0nY9siC+keHsktK+gDb5/tLJeVXhvInNezPT4dtby3Z35tC+WkVfSrjhNB2w9H2+0LbvykpT/9/D8hsq9OmrrxPErpdNL3aeXSDeLRLkG5eEqajXRNdZHPztrnY5oemdfogQTbnJWF4m1u6bjBv7W7HbqiIiaXHd33FduCL7zzxaJf4zjtO+qJJ/NkVzfm6/g6NbE42NzRL1i1BNldGmZ+5G/Ox3pbb3mROnrNedSQoFvSQoDHqIUG6eUlQrtBDgmzOS8J8bG4qusFybS5BulVRl2eNLQ+79Bxseuw+8rBTuMdmCjnYuV8jqkiYzvVVPuQ8bK4LsfmQTevEgGzOx9x16zuHuuQxmqB5wUuCcgseEmRzHhKkm5eE+YzVpnX6IEE25yVhPjY3Fd1g/jZXlSOty5/uCX3ZwVooWC9p+DDgX4DPYk90ZrkNS+zdO7f9cOCfw37fBfw8xW+ua1oPbPWQ72ErtaX8eWh7Skmbfwj9OzSz7ZCwn+vZd7W2lAvCPh9ccOwvlRzn69gb+u5ZUHYd8K3M3+lqcRdgmuc/b2en0Xj6fClwF3BgQf0fC/WLVs6D7fO9me3Ec5bdmE18saR9UX/eFbY9oqRNmnB/Qkl5Fc8Kbc9ytH1caPuSkvKvYTqk1GlTV57yfOwp5zuw16rmE/ZFbIa+Nv1sNNinV7u2ukF37R4NvBu4MRz7Zxv2dZNl63Y28CnsOnQrpuHxDfu7yXS1azIW6+qciT1c/a3w+XuKv5DKs8l0dYP+rnNgb7TdwlYdbcIm/WoX0/zQpM46O8+xyZfjmwXtlmZz98cWL7gVc4qvwr50r2OTOGxuLN02KT7P19b0t6zdUmxuF/D7bPtyXw7Hv0dBXbBVrPeWnsHwLD2+6zO2g/bxnSe2g27xnWeceOKVMur6B+2vK+fn2k9tzgfFhLHHNnkbyzM1m1Ns47O5dXz+OcRjc2PGhB4ffZNl21xffubPhL4/I7e9yZw8Jb3A77NssuxYUHnWYX2RTZatm9cPgXjGKowXN6S0ybNusmybW2c4H26jZn8x2Zx8OGMVNrdJ8Xkqx1pvc33nWWPKwy49B5s99lLusfGMI69fXoTivW02me71tQ8fUvGecs9jxC1D+JB1usG0bG4Tnw9Z1XZom5u6D+nNoU5xjOo7yZ2saj4F5RaUQx02h7rJsnXLovtU90X3qeo+1SrGmB828fm/Ze2WYnNtfV+ozpHW5U/3YOcINQfJ89pQ/zexE8pyFfbA4onAhzPbbwN+AjgHeDLwlND2fcCLsCRlm3oAP4wJmd32BOyf8+mSvt8f+AK2GknKiWE/l5ARJMPJwDexf0j+2O8rqH8UlsB9B5Y0znII8IPAx3J9BviFkj6n3NChz/fAnr79PHBnQf3bw8/dBWWwfb7vxZLJReX3oHwllyINH4M5Qp8paXP/8POGkvIqTg4/L3O0TR+ALOrXg7GnkbNPIjfRpqocLOn/aiyY+yjwn4CLgeOoPv9XAYfltp0IPBWbvDZzZZdX7CvFq11b3aC7dgcDVwBvBC5s0del67YGvA4L5PYDfg94P2ZvX6/p75S1azIW6+p8BfgdbB7ZD3MyLgp9qXqV8pR1g36uc2ArpD6Hdq+d7lu7mOaHpnWuZd+VN/K+VxFLt7nDMF/ro1iy5dZwrFsa9DcWmxtrrJ7Cvjc7HI/5u/lV/PIs3eZ+G0v2PQuLhx6OBat3YgFenoPY9r9jYOnxXZ+xXdpnaB7feeLRrvFd23HijVfKqOsf2GI2d7TY5425v6c25ysmHD626WpjeaZmc4pt/P6Sxz+HeGxuLN28PvrSba4vP/MM7CbZ9+a2N5mTp6RXF59l6bGg8qzD+iJL183rh0A8YxXGyxVC+zzr0m0O/D7cnG1OPtzqbE45Vr/N9Z1njSkPu/QcbPbYS7nHxjOOvH55EYr3tpny9bUPH1LxnnLPMKzNeX3IvnWDadmc14eEeGxu6j6kN4c6tTGq7ySHvU8VlFtYQznUIXOoS9ctRfep7kT3qeo+1SrGmB+UQx3m/gGozpG686drWIJtvaDsmTR7GvQ3KvZ/APZE51tD3dsoXhGsrt4LwvYk/L07/F32zz42lP9lbvtZYfuvFrQ5FLtgfjC3PT32cwraPD2UnV1Qdmooe1Wuzx8q6XMZbfv88FD/z0r294BQ/tGS8vR8n1dS/vxQ/isN+3NoqH9Fyf52YcHtzRSvLlfHR9i5yl7Kv8UG1k1YsP1Z4PGZ8reEvh1Z0DbV4aUF28q0qSsH+CTwv3LbvgCcW9GmjCQcb83RFvzatdUtu72LdilbdFvhL2GZusH2ao5PaVg/T8I0tGuiS1vtwALfNvVTEqahW3Z7F5u7D5bkfCy24kHTlWOKSPBrF9P80KTOeuhHHyQsx+Zeys4HgLqQMLzNxTA/gPms/wefP5awHJt7DztXw3lT2J5nf+y8yt4MPjSK7/qL7bJ9bhPfeeLRrvFd23HSZ7zSpH99UHUNgfj8TMWExQwd23Rhan5mEYpthvXPuzIl3aBfHz1hOTbXh5+5C/O7Li5o02ROnpJeffssCcuMBZVnbV6ep4svkrBc3cDvh8DqxirEFzdAf3nWhOXY3Dr9+nAJ88hPy4cb7jqnHGuzcug3zxpTHlY52H2PvZR7bDzjKEtXv1zxXjUJ07i+9p0HVLzXvDyPcs/N6/T9XX0XpmxzXXzIrkxNt7FzqFMbo/pOspgh51NYVm4hj3KozcqL0He587lPFeKLG9bRfapFLOE+VYhjflAOtVl5G98XqnOkTfKne0J/gGZvOLw38ApsRbE/zzbOcBRwGnBSxX7uwlYoeD/2D/oJ4Ah2rpBRVy99mjRdfW1v+Ny35Lj/NVc/Je3rpwvanISJmX/yND120ROpVU+rpsdK+5AOisML6lbRts8nhp9lieKvYk9VH1NSnp5v0fGy5WXnnO/PnZiBlv2vEmylvJdTbGdV7Ic93ZtfZQ/gB7BXS/8D5vTfgtlVtt6xmI1fn2t7INuDN2tDTbUpKz8As5lX5rZfAvy7kjaroot2bXWD7trFwtR1OxQbo3WrsK2CIbVroksb7XYBP4cFwR9vUL9Ppmhzfwq8HfgA8OKKeqsktvmhaZ0HY6+tvgv4BPZl62ZF/VUwNZt7GpYIvgD7Yv1G4H9jKwi39Su6MDXd8hyA3QhxHsPqBtPT7uPYTQrHAtdgK7E9luKk/DHY+TVZPWfVKL7b97hdYzvwxXeeeLRrfNdmnKwiXqnrX1eqriEQn5+pmHD6sc1U/cwUxTbNy0H+eRFNdJOPXswQfuajgfsBf13Qpm5OnpJemk+nn2OF6Ws3RV9kbN1i9UMgvrghZep5Vvlw8cQN8uGGuc4px9q8HPrNs8aSh1UOdpul3WPjGUd9onhvdUz5HocxmZLNFTHFeA/GsTn5kMXIhzRi9SG9OdSpjdGY5tQp2loMxGZzHt2mOKeObXOx5lBj1S32/CnEGzdMPYeq+1SnOz/I/21e3sb3heocaaf86Rr2z1rPbX9V2P7yiraPDHWuymw7CfjBgroPAb6Bibx/i3opl2OvcMw+LPm5cPzTcvs4k+2V4fJlV2APwB1QcOzfDG2ekdt+eUWbvw1tihKcG6Hs+Nzxt4B/X1AfzDB35ba17XO67VklxwCbYLcwvfNUnS/YYLgTuGfD/oBNRlvYBTbLacB3gC9jKw1k2WDfFfeKOCbUeUtB2cXYk+T7F5SlXBraH53ZdjDwF2zbUNZO67SpK09XwXl0bvuLsdcTtyXB/5R2F+3a6gbdtcsy5kpsU9YN4K+wL0ry15mmJExDuya6NKlzAnaN2ovNTU+uqFtFwjR0g+429xzsC7q0fA/jrBwT2/zQpM6TsFVdT8Dmy7/DkpnfV7HPMhKWY3N3hM+5mH/5bGzcnllSv46E4W0uhvnh57Br3QMa1s+TsByb2x+zt7uxIHILeElJ3WeH8hcUlG1Q7+v2ieK77eP2Fdulx28T33ni0a7xXZtx0jZe2aCZHVf1rytV1xCIz89UTBhPbONlin4mKLYpYkj/vAtT0w369dETlmNzffiZf4StFnxESbuqOXlKevU9n8JyY0HlWZuX5+niiyQsS7e+/BBYzViF+OIG6DfPmrAcm+vbh0uYR35aPtww84NyrM3Lob88a135BsPlYZWD3ffYS7rHxjOOslT55Rt0z8FO6RoRU7wH07zHIUXxXvPyPMo9N6/T93f1XqZsc119yC5MUbcxc6hTG6P6TjKO+1RhWbmFPMqhNisHfZdbxBzuU4U44wbdp1rMEu5ThfHnB+VQm5e38X2hOkdal1+Flm84PAEz/n8Efq+i3mcx5/uhwEFYsvIsLAl3KZasvAV4EPAzoc0vYSfdtB7Yk5/HYaLuzRz/XGxluPdijslNWDLxh7CnOI9h3ydFDwx9vQp7IjtP0Upq6bGvLGlzMpa4va2k7Hbg6sy2/xL6eyG20tyVmDH8m3D8ewIP7Njn9Pc/wBKx/xfT+G2ZOhdik8UTsVeSZo9Xdb4HhH1egRluk/4AvAhLHL8HeCv2VPcPh+Nfj71u9Ju5NumA3Us5+VX5Uo4Efgp4FNt2VMTfAqcAH8JWrTkES6BfiU2e9wK+FOrWaVNXnmUr9/d+BdtWTRft2ugG/Wo3NlPW7RVYEuFU7No9NENp10SXptpdi61oeRh2zXwz5qT09TrxJkzJ5o7BXgl9akn5kMQ0PzStc3Hm96uw1TK+hDnI51X0tW+mZHOwvQrt2eHvyzCn/kxs9ZihmJpueX4Zs8EbG9bvk6lp97PAf8RuRvgcNk+8GltA4w25uk/A5tx3Fuynia/bF4rv9j1uX7EdtIvvPLFd9m9PfAftxwk0j1ea2nFV/7pSdg2BuP1MxYT7MqXYZop+Jii2ySP/fLVjVT76TobyM5+GrXR4c8kxqubkqekFmk/zTCnHCtPWbqq+yFi6xeyHQJxxw1zyrPLh4ogb5MMNNz8ox9qsPKWvPGtd+VB5WOVg922zpHtswJeDbUofOdgpXiNiiPdgmvc4xMAUbS5lqvEejGNz8iF3Ih8yfh/Sm0Od4hiFOObUKdpaDMRkcx7dpjqnKoc6Dd2mkD+FeOOGOeRQdZ/qdOcH+b/NyqGd7wvVOdK6/Gola5gDuZ7Z9uGw7YwG7a8OdR8V/n4alqC8FvgWJsAm9srQH8q0a1oPLOG3BfxJwfF/DfhiaH8DNmCPBL7GTmNI9/P6inP5Nvs+mVrV5shQdmFB2YHYA3mfKCg7BVsJ5aZQ5zbsgv16dq4W5+kz2NOn12JPUm9hidEsB4Tjf7Ll8bz9AXgM9tri74TPVcDvYq/NLuIyzDbuW1IO5pRusfPNiU/FEq5VT2gD7MYG3o3Ad7FXkj4Pc9ruxp7UTemqDZjue4H/kNv+WuwC1JYE/1PaXbRroxv0o12WrivHJCxPt/+BJUqOa1C3ioT4tWuiS1ubS3k/Nk+1JSF+3aC7zSWhfG/msxWOtRebG9uS7nOtZbuY5oemdYr4IPA/W7aB5dgc2JfS+XH5i9iXwR4Shre5sefVI7Fg4qkN6xeRsByb+0fg13Pb/hsWk2S5D3ZjxEUl+2ni6/aF4rv6+t7YDprHd11iKW98B+3GSdt4pakdV/WvK2XXEIjTz1RMWMwYsY2XufiZim2G9c+7MEXd+vTRE5Zjc139zFPCMX6r4hhVc/KU9Op7PoVlxoKgPGvT8ix9+CIJy9Mti9cPgf7HKsQZNyT0m2dN97fWsh3Mw+a6+HAJ089Py4cbxuaUY21entJXnjWWPKxysM3azPEeG2g/jvJU+eV95GCndI2IKd6Dad/joHivWXkW5Z73ZYzv6r1M1eb68CG7MEXdxsyhTm2M6jvJYoaeT2FZuYUU5VCbl5eh73KHzZ9m97nWst3Ucs9leHOoCcuwOZjHfaow7vygHGrzcmju+0J1jrQuf5qyJ/RpB2uhYL1mB1PiF7Fz+s9jdyRyzsZ0OmnsjhRwGHZB+UNn+ydjg/fevfWoPz4J/Glu23XYaoIxELN2KX0Ecn0Ts26vYdykaB0xa5fnA8D5Y3ciEKNuh2GrfmY/nwIuCL/vN17X/j8x6lbGbmzlixeP3ZFArNpdAHwkt+33gc+P0JciYtUtyzpma3VvQR+aWLX7GjtfLX82OwO5X8N8llML9tHV110aiu+a0Vd81zReaWvHY8SfsV5HFBP6iTm2iVm3IhTbtEf+eXPko/vo6me+NGx/UM1x2s7Jseql+bQflGdth3yRfojJD4E4tVOetT/kwzVHPlw31lGOtS195FnrypWHbY5ysM1YVY6zzC9fdQ421mtE7PEexKtdFsV77VC850M+ZDfWkQ/ZlhhzqDHrFfucGrN2KZpP26M5tR9iyqHGqNsU8qcQp3ZFKIfaDPm+3VlH/m8bmvq+UJ0jrcuvpuyh5oHD9HNNzY5iYRfwrwu2Pw57e94N2KsqRTm7saet3z12Rwp4CrZy3Pc7298P+Dp2cX8YcCzwHMyRGJufx1YL/BXgocCrMJs9csxOZYhVu0OwV8GeyPaXHScCDxyzUxli1e112IqHj8XGU/qJ6foYq3YvwybXo4ATsGTL3cCTRuxTllh1y7MH+OOxO5EhZt1eib0N+EHYqrbvwcav5odqTsFWs30R8BBsdbZvYq+qj4FYdUvZH/MHXzZ2RwqIVbsN4CvA6dgccQZwK7ZCW8pB2Go2by/ZR1dfd44ovutOX/Fd03ilrR2PEX/Geh1RTOgj9tgmVt1AsY0X+ed+5KP72KCbn3k1cHmD47Sdk2PVS/OpH+VZfcgX8RG7HwLxapdnD8qzNkE+nB/5cH6UY/WxQfc8q/Kw7VAOtjt95jib+OWrzsHGeo2IPd6DeLVTvOdD8Z4f+ZB+5EP62CC+HGrMesU+p8aqneZTP5pTfcSeQ41Vtzx7iCt/CvFqpxyqD/m+3ZD/254N6n1fqPZ/6/Knh7Pvs4SFDxwehT0tmn7yT0HGyglYYu2dwHnAH2FPDW9hQsb41r4YeTRwDnDw2B1ZAT8OfAwLkr4BXEJxAn0Mng9sAncCn8H+DzERo3Zr7LygbWEX01iIUbcizWJ8q22M2m1gzs2dwC3Ya+qfOGaHCohRtzx7iC+Qi1W3t2CO3V3APwEXEt9qT7FqdzpwBeYbXgecRTwrFUG8ugE8AZsXjh67IyXEqN2hWCL+euxV81/CVkLcnanzUGyuPWrgvk0ZxXf90Fd8t6p4ZYz4M8brCCgm9DCF2CZG3UCxjRf5592Qj96eIf3MtnNyjHqB5lMvayjP6kG+iI8N4vdDIE7t8uxBedYmyIfrhnw4H8qx+ujD/1Ueth3KwfZDXznONVbjlyveG44YtVtD8Z4HxXvdkA/pQz6kj1hzqLHqBfHPqTFqt4bmUy+aU31sEH8ONUbd8uwhvvwpxKmdcqh+5Pv6kf/bnia+L1T7v3W+8b3Y91nCdXdvI+QY4B3Yhe4O4LvA54CXA0eM2C8hhBBCCCGEEO1QfCeEEEIIIYQQQgghhBCrQzlYIYQQQgghhBBCCCGEEEIIIYQQQghXkac+AAAAaElEQVQhhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCFEE/4fdeEOf6BVi+8AAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\left [ Assignment(rho, f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8), \\quad Assignment(u_0, (-f_curr_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_curr_8)/(f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8)), \\quad Assignment(u_1, (f_curr_0 + f_curr_1 + f_curr_2 - f_curr_6 - f_curr_7 - f_curr_8)/(f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8))\\right ]$$" + ], + "text/plain": [ + "⎡ \n", + "⎢ρ := f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr\n", + "⎣ \n", + "\n", + " -f_curr_0 + f_curr_2 - f_curr_\n", + "_6 + f_curr_7 + f_curr_8, u₀ := ──────────────────────────────────────────────\n", + " f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_\n", + "\n", + "3 + f_curr_5 - f_curr_6 + f_curr_8 f_cu\n", + "──────────────────────────────────────────────────, u₁ := ────────────────────\n", + "curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8 f_curr_0 + f_curr_1 \n", + "\n", + "rr_0 + f_curr_1 + f_curr_2 - f_curr_6 - f_curr_7 - f_curr_8 ⎤\n", + "────────────────────────────────────────────────────────────────────────────⎥\n", + "+ f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8⎦" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "for i, u_i in enumerate(u):\n", + " moments.append(Assignment(u_i, sum([ (c_j*f_curr[j])[i] for j, c_j in enumerate(c) ]) / sum(f_curr)))\n", + "\n", + "moments" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAASCAYAAACq/vK0AAAABHNCSVQICAgIfAhkiAAADZxJREFUeJztnHu8VUUVx7/3CgmoyEOStMRHIfjGUEQzMMwPftJ8hWYiHt/l9YUWlOYrxAQNUfNRqaFlL0208p3vSJRAAtNEhQv4AOFeBUTkefvjN/PZc+fMnLPP2fuKfT7n9/mczz5nZq21Z6+ZWbNmrdkHaqihhhpqqKFCNAItkc+iFPybA28Z+vvL0D7uyd8ArABeB/4AHAHUtQFvW2Mc8ASwEFgFNAMvAZcB3VPKqESPAIOB36H++xhYjnRxPzASqM9Ib/V9UOT+2wFjgKnoedea61Skj90DPG4f7h2R+ytTf1QJ3nMjvAC3O3S3l2lD7PPjjG1ekuIe7ue7JZ4nb3wDeAyNtVXAXOAeYGAJnqy63xRYA6wENinBP9vw9vXKnzTlB5bgbWt0B04DJgNvIN0tA/4BnErx/HERGnPNwDTDG7NbOwJjkS1pBlaj+TsJ2DNAv7mRPS8iby/n/kMiNLNM/e4A7QIEy4CJgfIPIwJdXA5sS+nJZLG3oRtjrnVAZ6APcCRwHBoYR6KFIC/etsZIYAYaFO8BmwH7Id2cYb4vLCPjctLpcRPgFuB01D8PA/PRQPkiMBToD1xXJb2F1feMQBsuQIO4Axpcf0KDuQuwD/ADYBQwHLjbk2lxbER2f3P9V6A969D43SPABzAAOBlYj57bl+HKGRuRAXBfhjZvBtzk1bcDLkYG86cB/kdKtCVPjEP90oSchKVoDBwBHAOMAH4b4Muq+z2B9sALpj6EzdAisQJ4zSmvM/ffgIzmxsIwNI/eBZ4CFgBbA0cDtwGHGpqWAK9t/5Wmvh7p/VuGtxdwqUNfhxyWi4HPAM+gRWklMvgjgBOAM4E7HL4PkdO2ZeQZRjnfQzRD0SLxEFq4i9BoPtVgFzQBHgZmIkVsFaHdydT/N1LfE3VCCzI+efGmRcHwD66Ct0OkfKyReXMZ/kr0eKmpf4Bwh3cDDs5AD4m+XysmZzxJXxwQaWNf9CyDAjKnAW8Dbwb4rAe6ONKeKWi3+2KAtx4ZqMVod9OCFq6QnJDBDyFLm13saeSEFq9KUKD6MdoTGepFwGe9uoOM3LkBvjx0f5Ypm1CifQcamqe88t6m/JUSvJWgQHU6/BpwOMU7iJ5o4WhBC64Pq7//BOq+Y+rcsVMH/Jpk3PUO8A1Bi8I6oJ9Xt8SU++hleOYb2acEaJ4wdYMCdUC2xeJJ04C+wF3mRl+P0B5n6u+O1AP0QKtnC7BDTrxpUaD6iRiDNRKPl6FLq8c65A22oEFaDpXSW8T0fTLJxO2R4t7uLtbKvBm4wXzv7/EMMOUPRtpzA1qEVlI8ab9naE5CE2YNMuQhOb8o0/Y82uzi1ArvG0OB6seobecDkfrlhHfkeej+DlN/fIn2XWBorvHKjzflvynBWwkK5D/PLzIybwzUWf3dEaizkQTXuP+IxLGIOaGQjMc7vfI5pnzzCP0p5jrSq7eRhKluYSi2tikKGVwEnIc8jVKxRdCqeBCaSK+SbFtiIZRYeMHFEuB5832/nHg3Jg4311klaCrRY2eSHMjaFPevlN4ipO+tgJ8h7/R4pO9S8CeBK/Me8/3YFPd1y6ejXUEn5LG5bRuL+v8Z83s2ivGG5Ewr0/Y82uziyylo2hqvIyO+L8W71q8CWwB/D/DloXv7/KX0vk+EJo1+Nzbs3Ap59Lb9oR3ZzuY631x3AK5AOcVh5hrDo+bq55o+MFc3itAdOSxPo3wVKGTswoaoxruFoZxFT4pX7nnIk3wmQL8F8gCaUKwdEoNYbrGYHqm3aDJX92Gy8H6S+D5a0bdEbf4K0svVEfpK9bgMJSY/j3YrN6FB81ZEfqX0FqEJOgLoiozmv8vwl5P5MvAOmhCjIjQx3pXm+x7IAIL0uyXQQGnDbOXsj3QSwjXOPbK02YVtU7nx25ZoBkajUNArKGfRhAz/N9H4ODPAl1X3HVGY9X3CYTyL/9fFoh2aGxDOPcUWiy4ohwSKJoDsR3vkOMaS1BY2B+qHlt93yt82389Gi/zVyCb4fNuj/MkcyhyuuQzF47Y2AncDbkUJmY8IZ90nIM+xwSnbxpTNCdDXoRVvPTKQpfCIkWM9uCy8laBA9u3pIlqfeHgY6TWGSvUIMnRvePeZj/rMj19WQ+/q293KTqP1iR8XfdBi537OC8j8iMRZsdvifR06exrmcwFee5JmR0NzhakfgMaqzQtdZepPjzyXfyrF/TQF6Ktps4t26OTMapSszIIC2cfokWjhcJ/7dbTD9ZGH7geasseIo6uhWeqV16Pw2DqUAM8DBfINQ11LPAzpjrkr0by4Ei0OTab8byiyU4eckVCoMwSbZ/JzIX805XbH0RFFAWyerg7N7UkOjx3Xp6W4bxBWCZO98l3RtutlikNVS9Hg6eyV72xkvZrivnMNrV2ksvDG0Ehpo+F/JqW4t4utkVF9DQ2A0C6hGj1a1KPQwVh0SmKtaecGio1kpfRW3/4gtPmgUO5jNMU6eygg83mnzCY0rzW/OyGj8DatYXmnOGXvo3FpE6tLUKIeZJRaKF4IrZznAu0PIUubXVSb3G4k/zE6yrR3AjL8ndDYfNTIGO/R56H7c0xZqdNnhxga3zPvY8pfLvtkYTTStvP8XBLb1C1Qb/XnftagE1UPAt8mOTZrncS1hKM/Pmxuw48I3WrKh5rfDeb3cQ7NByR2vRs6RfUuxTm+VA2xN70QGRkXNxkZIyk+Bjcbrdh7Ac865WnDSF9Ccbt3SEIdWXhjmEhxqGovdITwTooT/jPLyPOxGHXGDLRDuAvt2FxUo0eLDabc1nVDz3QiMmS3oQFSDX1o298FGRbQ0WAf40i21A3Az2kdTgjJnIIG6DB01LYf4eOuId6Z6IjfGSj0cTrylkHGbzXFBsbKSRtCy9LmkJxKQ1B5j9HBqI8mo2SyxQzk2MxB8/1WklNReejeOkqlTqDZHKMfqskagmrLed4AXI9CekNIdODCtn8ixQllHzaPtIJw7sNFHTo6C62PeUMShuqCxuaFKPx3r0OzjEQvDWjXNobiHF9qdEYGxE2yDCfdCn2+J+s6U15OYfalnoty4q0EBfLdnlq8ZOS6ScVq9VgKX3D40rwIGKO3+j7HKdvUoS0XCpxk6A4LyCx4tDea8v1Q2KoFuMSjsbwjnLIJaAFsRgbGemc7GNoXAu2yckI7rxCytNnFzRXetxQKVD9GbaTgnEj9fabePf6Zh+5fJLzTc2GP2/qO6cQyba4GBbLP8/ONjNkUH0N2EdJfDNuR7PY7laE9wdC+QvGBJbvLPxPtXOx3F7PQ4t0BObbLiLybkXZnYWNe1svojJJ/a9HWpyXAsz1aZf2Bkca7Go6OdTXS+ghaFt5PA7YxV7t7yKLHUrAx82bCXk5a+pA3txoNsD3QkV7fm3Fh4/nldhagZPnZKMfUI0IT6v8ZyEh1MfxWh2mS21l2Fmnb7OLTkNyGJMQQO/Jsy9c4ZXnovr25xozqASj3MReFSF18Gk6R+RiNEsUz0Vzw8ywurP7SvEy4AIUztzVyY0eceyMHZB3KMWzw6t0E9xkojzrJo1mO5v9JqF+uIUl8R7Er4VhbL5T0cj11u8qPC9Bb9CdZcS3q0dYqlqDubhprvRX3DdEsvJWiQHUeRx/Ccfx6kpfy3JhvtXo8Eb08F/prgB1JXua7pEp62+YVaCHr6PHYdwXmobdPQ+iH+mq+U2Zlhv7qoR6FDReiuG8LrY2Z5f2Q1h5UN5SoHezJu9rIODlwHzuO/OcKIUubXeSZ3IZsXvGxhncRMkguDkVzaBXJLjMv3dsY+uMUvzcwkOQlsaO9unpz73Wk67O0KFC9Di8hyT+F7KYLq7+PSe+gX2jkv0niZLo4jOSlu0JEhu3nF8x1dIDmQWQv56CxGboX0Lrhw4Aforcm56GH2wn9f0wHlKS8FsUoG9AE+UlMMIpVrkcvl3VEg68vOlWzAikDpMiuaLEaaO41BRm4eY68LLyfFIaiBetZ1MlNKME9CBnlRSQhiCx6PAUN8AXmXguRbnqjRWET5HFcZWRUSg+JvmeZe7q4HcWfz0Lb38fMtQUNtn6oT9bT+i8jrMx/Upyb2YB2KfY02AJav79heafQ2oNqJnzEL+aJWjnLCU8ei1+ihSBLm13sisbndFp77BsD96L3KA5Gi9xkNDb7IiNUh2yBPRGWl+7HGPkHo/nxHNJFb7QLtQbN3632RbH0lcSjBUtNmz8JnITm7Hr0DKH/yGok8eKt/qZTPgdhcR3KqQxH/5DwABpfPdAObBfUZ0cBf43IsDuLfdFu4ZYAzTJkQ7uilwXfSdO4QcDvTcM+QB7lEuQFjCDxSp9FnRo6NunDelsDzO8RFMfiPzINnIJioPtHZGXhrRQFqvM4dkPJ6plo8K5DnTENHZVzPZAsehyAJt5zyBv7GOliDupw/+WcSukh0XfoD/gsDkWG5l00XlaiEMJfUCzXf3/Byrw+Im8QSd/+uUJeH03oGf3dQGgc+R9395qlzS7s27JZ39y2KJAt3t4e9dFUkiOp76Hjm4d4tHnpHhTqmEjyYuAq8906ICGk6bNq/lOrQHU6vDxFe54OtP+2Ktp4DHLUl6A5thjZjpHET0ha2KhEC+H/IQMtIDY/0qeK9tVQQw011FBDDTXUUEMNNdRQQw011FBDnvgfsfRdXaS9UTkAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$5 ADD + 3 ASSIGNMENT + 8 MUL + 2 POW$$" + ], + "text/plain": [ + "5⋅ADD + 3⋅ASSIGNMENT + 8⋅MUL + 2⋅POW" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "count_ops(moments, visual=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAACmoAAAAXCAYAAAABMoIMAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJztnXm4HUWVwH8QCIGERWFEGJgAIgEFTVhERsALCKiIgjD6yejQKjgOCKijI+IID2UVxmERRhiVyCCibG58KCA+FjcWWaIEgsALS9hBiApE4M0fp/q7ffv23l1d1e+e3/fd7yW19K06t07VqVNLg6IoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoiqIoitIq5wKPATPN/3vAZORzZ4m8o8xqwGnABLAMkd3nXBbIA7ogk62Qcn3UdUEyUD0TutCe2qYrMlE96w5daVNtojIZRmUyjMokGZWL0iRqTyi28b2Nafvqo+PLMF2QSZ6OrcWgP3KypXLl0aO4n9QHPo2U8wOuCzJFqNP3hr/Ffo2WqHv4Pr76gupus6jdNHWw/Vt2wYbyFZWd4hN59obaI8VQe6RZ1B5phzbkrGNeNVRuim9k2QN1/ZalbI2tgZeRgTekZx4wDowBnyiRd5T5MSK3y4BjENlt5rJAHtAVmVwKPAzMcl2QBFTP+nSlPbVJl2SietYNutSm2kJlMozKZBiVSTIqF6Vp1J7wm32B04HrgGcR/T+vYN61gZcQ551LfG1j2r4G0fFlmK7IJEvHVkHKPYY48X3bqDlOtp80ZBbwoMnzA4vlSuN8892bOPjurpG3kbJu33ueef6mFfP7yJrAAYgu/xF4DngGuB5ZjFg+JZ+v42sU1d3uYFt3lT4nAj8HHkD0/SngFuAopD8owxgwv2SeNn7LrthQPqKyc8t6wLeAJcALiP18CvAKh2WyyYfob4o4ICVNnr2h9kg+ao8UR+2RPk3ZC2P4aSuAjnlVUbm5pY6vvqvUtRfq+i0L2xpXAH8CVo6E9cxDxyrkHVU2RWT2U9cF8YguyeRNSFmPcF2QBFTPhC61p7bomkxUz/yna22qDVQmw6hMhlGZJKNyUWyg9oTf3Ir8PkuBhZRz/nzMpN/JTtEK42sb0/bVR8eXYbokk6I6No5/GzXHCqY/2aR/GbjfTpEy+QekTSzn4Lu7RriRck5KfN2+9w7gz6RvXuwiH0dktgT4DnA8sknkTyb8IpLbnq/jaxTV3e5gW3eVPsuA3yB6fgKy0HsjIv+HgPVz8s+O/HuMwc0X65PfP9r+LbtkQ/mGys4trwEepb+Z7wTgavq3v5fdSO076yN9wVKyN17k2Rtqj+Sj9khx1B7pU8de8N1WAB3zqqJyc08dX30XacJeqOu3LJR/E2SgPzsW3iPfAZmWd1Q5CJHZga4L0hABUp9ejWd0TSYLgcXANNcFiaB61qdr7akIAfX0rIsyUT3zmy62qTwCRk/P8ghQmcQJUJnECRg9WzCPgPoy8ZGA7tVL7Ql/2Ql4LeLQ71HO+XM58AR+/K6+tTFtX4NMtfGlCbomkyI6Nk43N2q+Dlmkupy+Q3wtWwVTarMQceAnLf7V7XtnIjdFX18xv6/sDOzJsMxejWwmmAT2Scnr2/gaRXW3W9jUXWWQGSnhxyJ6cmZG3g2R3+kUYHX6my9mAF9EbvXZPiN/G79l12wonxhF2QX44z/5GVKWQ2LhXzXhX6/43AB/6hiyHHAVcA9wEtkbLyDf3lB7RGkKtUf6VLUXumArwGiOeU0wqnIL8GcsreOrzyLAnzqGNGkv1PVb5uY/wWTeJRbeI98BmZY3ZCcTfzKwJXKi5ynkVSiXIs4bEKPjfOAxE/cT5LRGV9iH/tWp8U+XX2sTUF25uiqTo5Ay7u66IBFUz7rbnooQUE3PuiyTrunZKOgYdLtN5REwenqWR4DKJE6AyiROwOjZgnkE1Jt8+jqmBFSvl6s6qT3RDXoUd/6sjryq7ZxYuLYxQeeFwlQdX+rQVZkU0bFxurlR82rgb8irvM41+XbNSL8DcAniSH4e0c8bgOMqpg37g6/E8s4ADkdueHwecRgfgTiNlwK3x9Lvap5zPLA50nc8itwO+Stg24Ty2eiL9qG/kX8ZcHek3HHKlPlE0nXngyZN0b73JOTGhB+a+k4CrwfeYv59mqnzecAjKeWJshxin10DPIm8OvAWYP+U9EVoa5w4wnzP6Snxvo2vUWzpbtF0o6q7ZcvbhO6G7XTvhLjZJu6SlLxdoQ2df6P5jitz0q0HfBO5Tes64HfAHxEZ59kqeb9lHeraUFXaUd64kRfvC6MsuwA//EIbmefcx/AGsVWRPvQvyKGRsgT4UccohyEbsXZE7PFJsjde5Nkbao+oPeLDXMK1PdLW/KCIveCrrQCjPebVoQlfVVnZ+SS3AP/GUvBno6at+jVpLxSxFcZJ91sO5Y8bbG9DTvb+JuML0sjLu6X5uwlycvglpJO9H9gLufp4T8QQmAV8G1gE7IEYH13hYeBo4HHgRfPvo5Ef/253xXJKV2XyS/M3y+htG9Wz7rYnm3RZJl3Ts1HQMeh2m7KFymQYlckwKpNkVC7JTMUxxVWd1J6YeuwBTGfYOaltTNB5oaDjyzBdlYlvOtYU+yEO5zOR0/MLTPiWKemPAK4FtgJ+jtyA9GNgJeDtFdOG33VLJGwm8Atk0fGvwKnm/0cit5XMiqUHmGf+boK8si7sO64BtgN+hGwEiNJkXzQN+C7yGu2NgQsRub6M3MryLYYpU+abTRzIwuvRkc81Jrxo37s5sqD4MnAWspBwJ/JbgSwi3YQcSjg3pTwhKyOvhTsHWMOU8RxkoWI+1V/T2dY48Tfz98WUeF9135bultHxUdXdsuVtUndvTojb2vyNy7VrtKHze5q/8c05cR4EPorc+rcd8pufBrwX6SuzqLN+mUddG6pKO8obN/LifUFlV52mdHNn8/cKpK5RliLj7SrAmxuvQT5N9z+bIRuxTkXG1CLk2Rtqj6g94tNcwpU90tb8oIi94KutADrmVaUJX1VZ2U0FuYVMdT+vjfo1bS/UtRUy889EFGNBQlyP7JPiWXlDvmOe8Qjwhkj4qsDTJv/DSIcbMh35AV4m/YpkH5mGGCd5k9IuEVBvp3cXZbI6UucbXBfEoHrWp4vtqQgB1fWsqzLpmp6Nio5Bd9tUHgGjp2d5BKhM4gSoTOIEjJ4tmEdAPZn4OqYEVK+XqzqpPdENehQ/pXshcktBvK7axnReGGcqji916aJMiujYOOkn09umR7afFETfHkJua3mFCdvd5PteQvq1Ed28DtHHOGtVTBv2B3MiYeeZsC8itzWG7EL/VovDYs+8wIQ/xvDi8MUmbodYeJN90dfo32izQiR8RcThPYncrlCnzB8zYR9jmDJ971KSN0LMN/GPFiwPyO0RkwxvyPw75HbN5+m3rzK0MU6sgMgr69YJn8bXEFu6W0ZvYXR1t0p56+rufciicRLHm2e/KyN/F7Ch859BxsH/Rtr1JHAb0j9lsR7wDQZvyboH+X2zblQq8lvWpY4NVaUd5Y0befE+MaqyC/DDLxS+zvPfU74n7Iv/rUIZA/yoI8g4chNwF3KYBYrdkJVnb6g9MpwuZFTtkSpl7ro9Ymt+UMVe8NlWgNEd8+pS11dVVnY+yS3An7E0Sg8/btRsun427IW6fsvM/JuYyCsS4npkOyCz8oYsNGnipzFAdjdP0r/+OcqVJm7NjGf7xhZImec7LkeTBNTrQLoqk+eQTsEHVM/6dLU95RFQXc+6LJMu6dmo6Bh0u01lETCaepZFgMokToDKJE7AaNqCWQTUk4mvY0pA9Xq5rJPaE/7To5jzZwbiJLswIU7bmM4L40zF8aUuXZVJno6Nk+7wbJse2X5SkJtpJoGDI2HrmrBFCenfauK+WeD7y6RdiPSp4VuV3mzy/iAl/WKSF0rvMuF7JOQ5xsTFbwZoqi/aFnG+p5U5XBT9cM0yf92EbZ2Qvkzf+5GU+NtN/DsKluddJuz7Kc/7tomv8lq/NsaJk03ay3LS+TK+htjS3TJ6C6Oru1XKW0d3X2nif5oSH9Zp3ZT4rmBD5x+hvyFnEnmV7No55dgQadenIouVY4i9MgO5ie1ZYPuUvEX64bpUtaGqtqO8cSMv3idGVXYBfviFziZ788GxJv7zFcoY4EcdAb6E3LIV3agxRv7GC8i3N9QeSWZU7ZEqZe66PWJrflDWXvDdVoDRHfPqUsdXVUV2PsktwJ+xNEqPYr76IgRUr2PT9bNlL9T1W6bm385kTDqJ0SPbAZmVF2QH+0vAvSnxTyGngFdMiFuEdLpRDkJ2TT+PXG8bH/Bdsz8ij0Mdl+NgxAn4rPn8mmSDIs4Eg4Nm3md+gWf6IhOAdRBn5uOIQixAjNMkHiL99Txt06ae7Yhcdb/EfOe+FcprE1/aUx05TdCsnvkiE5AJ/41Im3ockdHmGem7omdNj2Vl5dQ2PrWpquP+BFNXz0KOQMr0tYLpJ5iaMhljuOxFnWsTCXmngkwgvW5nVMw3FWzBacCX6fcp9yHOtRWyMtG8THyZH03QXL2arlNZO2tU7Ymqcz4X9Cjm/Hm3SbdfLLxJ2VWxw3xpY237X0LK2hxt4cv4AtX64gnK9cNFnKe+yKTsmJunY+NIvXygh5RlLCX+9cgrn3+PyCHKE8hC4Wqx8LWAP5nn/gh4P+k3JRZNG/YH10XC/s/k2ybl2b8z5Yu+enCWec5iBm/NCTnfPHOjhO9uoi8Kb+05H5F5/HORiY8usFQp8w3AMuS1j3GK9r2P0l/IjjIDaRP3pORPKs+PTNhWiTn6C/i7pcSnYWuciHKoKdtCZFEtC1/GV7Cru2V0fFR1t0p5oZ7uvs3EH5MS/ySi11Fsr01N0KxtYFvn1wb2Rja0LCH9lbwhsyP/HmNwjrk+yX0o5P+WE9SXW1Ubqko7yhs38uKhum95An/sTxeyG2O4jkV8iBMJ+XzxC+Vt1DzOxB+eU6YJ/K3jmxB74SuxdGNk1z0kz95Qe2SYUbVHqpa5bXukyXX8NuYHZewFW7YCjKa9UGZfTJQJmrUXqsoNysuuCblNkFzPttfWbOpnj+K2XpQJmqtj0/WzaS/U9VsO5I86R58zf6u84iov7xuRRn5lQtwGiBFwCWJ4RJkFvIb+O9tBjIZTkQnx9chV7ZcjV1LfX77oVggHl1ta+K4VGZZbGP4gYnjfjRgS+yMnRrYi+1rhU4A1YmFzgfcgHflELO7WAuX0RSYzkfZ0PbKA+ThiSD2W8qyV6bdv17SpZzORq8fPQa4S94022xOkt6nVqS6npvXMF5msiAzsZyKOouWQkwtXIf30Uwl5uqJnTeoYlJdT2/jQplYE3kv1cX8q69nfkBOsB1LuVQFTWSZ3MXha66WCz5zKMtmGQYfe5kgflnRLXZSpbAt+FtnYtj/ilHgDMlF8AdlMkkbTMvFlftRkvZoeJ8vao6NqT1Sd8/nM3ohzO34DV5Oy61HeDvOljbU5LwypYnO0hS/jS1Wb9R5ko0dRlhRI44tMyo65vuhYE5yB+Fw/xbBNugDpg+YC10bCn0BuBjkKeCewp8l7JfAFZNGzbNqwP4jm3Q1xct+UUvZ1kDFlaSRsrnnOFSQ7nbcEnkE2LoU02ReFGxE/kFLmkKiulS3zCsgNH3cgbTRO0b73MmRxOil+BdJvd0mS4VuRRYebU/KsY/6W9YfbGCeiHIz0h3cgt33m+Td80n2bultGx0dVd6uUt67uhhuhk/RsI2SjcfTmnjbWppq2DWzr/KPApUh7XQScS/ZmwcUZcQ9kxOX9lk3IraoNVbYdQbFxIyseqvuWfbI/Xcmuig/RZ7/QM+bv6inftVosXRq+1nEFZIPgIuR121XIszfUHlF7xOVcokpf2OQ6vm1bAcrZC7ZsBRg9e2ENyu2LidK0vVDHV1VWdk3YCr6srbWhn2XxdU3Jtr1Q11ZIzR9emX19QlzPxI2lPDQrL8AnTPyBCXH7mLikK9d3MHGnRMJ+C/xvLN3dwPEp3+2C6xg+MRKyHuKYfgTp3H6PXI29K7IYFd0NvD5S/9dE8k4iN4tciwz2/5IRnsRTwL9WqFNgvqNXIS/4I5PjKN4hLW/KnHbavW3a1LModU/i2KBKewK7etaEnAKq65mPMgkJT5/tmRDXJT2zpWMhWXJygYs2lRTW9LgfMDX0bHVEb3ZGTunUud0qoPsyudc8vykCui+TpP74FOCPJJ8EziNgatiCP2H4BN+3TXhZAqrLxOf5UUC1etkcJ/PsLLUnBqk657NNj/xTutOQhYPLE+Jsyi7PDvOpjbU9L2zS5rCBL+OLT74qX2RSZswtomPjZJ9Mb5Me6X7SD1LsRoFPZjx/OnJjxPdN2idIvpklL23YHwTm/zPM/9MWRjY18d+NhYe3I348Ic+qSP/5i1h4U31RWOZrUsqcRtkyv8Gk/1bK84r2vWnj70Gk396QVJ5VTfrbUp43DVkIe5Tytr3N8fSTJs0C4FUFyuLT+Nqm7ualG1XdrVLeurp7gYmfnRAXyuK4SJhP431R2ph/hISvJlyrZBmLkPdbNkFVf07ZdhQNTxs38uKTcOlb7pLsxmjOhxjgh1/oABN2Vsp3/czE71KhnAHu67gGxcbotD4rz95Qe0TtkThtzyWq9IVRJqm3Pt2mrQD27IU2bAXo1phXZl+MbarKDcrLzoad5WptzaZ+9kwal68+b7J+Nu2Fun7LofzRGzUfRnZRz8l4eBp5ecMd0km7nLN2QM8zf8PTGtNN+pNj6a4A/rFQSe2zHLLzN35iBMTJ/WukPvsiu9W3N+l2QE5aRHcDzzVx90b+D/AZ5LVjdwNP0697PDzKNOB9yGTtV1UrVxGfZHI4suh3PmKALgG+gZxQiivOHFP2IrdEtUFbeuY7VdsTSHuxqWeu8F0mqyIDUNJJ3i7pmW0dy5JT27hoU29JCFuKTG59GPd907OzkddyXA0cWbdyFfFJJochTsmHkE0Gv0GM94laNSyPTzKJ98fTEWdf+IrENvHJFlwPcbJtCtyJ3DSxM+0vrE3F+ZFLW1TtCcHlnK8pdgTWRE7Zx7Epuzw7zKc21va80AebIw1fxpepYLO6HnN90rE6rAachMgxfC1gnA2QRfJ5CXEhy5Absa5CFjO2R14Vl3RbW1basD8I9fpF80l7veF/xNKHhGVNujlnHtJ/xvuVpvqicBGk7CJe2TKH7Ttt4blo35t2u1CWPJLK8wLixE/7rQLkxqITKW/b2xpPPwecgOjxrshifx6+6H7bupuXblR1t0p56+rupsjvHr+5aSX6C7c+zr3K0OY8bV3zt+gbTspQZ/2yCHX8OWXaUUjRcSMtPglXvuUuym4j3PsQm9TNcMPYbvQ3AoSsivjdn0Pq2iZN1fEF4JsZ3zEP2Zh1F9Le4uTZG2qPqD2S9py25hJV+sImaduna8tesG0rQPfGvL0ovi/GJnXkBuVl17Sd5XJtbarv/2myfjbthbq2Qm7+i5DGtXEsvGfCxzIenpYX84UvII04TniSJ2mQnG/iwquPw53wO8bSHYkINC1/kFHuppljvvOChLjLkV3ryyfEfYfhmwaOZHDX/xeAvyKTCAqEg1y1/WfEeHoaucq8CgHVd3r7JJPnzed4RBk/jMjn4ITv/7Ap9ycS4ubTftuCdvQsTt5JnDB/kJGmSaq2J7CnZ1D/xBJU1zNfZRLyPWSCMi0hrkt6ZkvHQtLkFOYPcvI3iYs2lRRWdtwvQkD39exAxDAN2+I4bm7U9Ekm70BOV22BTDh/jkzMX5lThzQCui+TOO9D7MF1M9JkETA1bMHlETvwZWRyPwkck1eBFAKqy8TV/KgIAdXqZXOczLOzRt2eKDLnC/MHeQW3SI/8U7qnI47atRPibLax |