From 9cf82f641b8982df526107b99d883545ee02fe20 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Wed, 12 Jun 2019 21:01:42 +0200 Subject: Restructuring --- notebook/inspect_opencl_layout.ipynb | 612 ++++++++++++++++++ notebook/lbm_codegen.ipynb | 828 ++++++++++++++++++++++++ notebook/lbm_codegen_d3q19.ipynb | 1160 ++++++++++++++++++++++++++++++++++ 3 files changed, 2600 insertions(+) create mode 100644 notebook/inspect_opencl_layout.ipynb create mode 100644 notebook/lbm_codegen.ipynb create mode 100644 notebook/lbm_codegen_d3q19.ipynb (limited to 'notebook') diff --git a/notebook/inspect_opencl_layout.ipynb b/notebook/inspect_opencl_layout.ipynb new file mode 100644 index 0000000..521f93a --- /dev/null +++ b/notebook/inspect_opencl_layout.ipynb @@ -0,0 +1,612 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Thread layouting in OpenCL\n", + "## Setup PyOpenCL" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pyopencl as cl\n", + "import numpy" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "platform = cl.get_platforms()[0]\n", + "context = cl.Context(properties=[(cl.context_properties.PLATFORM, platform)])\n", + "queue = cl.CommandQueue(context)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "platform" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Initialize array for storing group and local IDs" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "nX = 16\n", + "nY = 16" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "np_data = numpy.ndarray(shape=(nX*nY, 4), dtype=numpy.int32)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "np_data[:,:] = 0" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 0],\n", + " [0, 0, 0, 0],\n", + " [0, 0, 0, 0],\n", + " ...,\n", + " [0, 0, 0, 0],\n", + " [0, 0, 0, 0],\n", + " [0, 0, 0, 0]], dtype=int32)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np_data" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "cl_data = cl.Buffer(context, cl.mem_flags.READ_WRITE | cl.mem_flags.USE_HOST_PTR, hostbuf=np_data)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cl_data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define OpenCL kernel\n", + "\n", + "Simply writes out group and local indices to the given data array." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "program = cl.Program(context, \"\"\"\n", + "__kernel void test(__global int* data)\n", + "{\n", + " const unsigned int gid = 4*(get_global_id(1)*get_global_size(1) + get_global_id(0));\n", + " data[gid + 0] = get_group_id(0);\n", + " data[gid + 1] = get_group_id(1);\n", + " data[gid + 2] = get_local_id(0);\n", + " data[gid + 3] = get_local_id(1);\n", + "}\"\"\").build()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Test output for a 1x1 work group size:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 0, 0, 0],\n", + " [ 1, 0, 0, 0],\n", + " [ 2, 0, 0, 0],\n", + " ...,\n", + " [13, 15, 0, 0],\n", + " [14, 15, 0, 0],\n", + " [15, 15, 0, 0]], dtype=int32)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "program.test(queue, (nX,nY), (1,1), cl_data)\n", + "queue.finish()\n", + "cl.enqueue_copy(queue, np_data, cl_data).wait();\n", + "np_data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test various work group sizes" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "def computeWorkgroup(i, j):\n", + " program.test(queue, (nX,nY), (i,j), cl_data)\n", + " queue.finish()\n", + " cl.enqueue_copy(queue, np_data, cl_data).wait();" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3],\n", + " [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3]], dtype=int32)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "computeWorkgroup(4,4)\n", + "np_data.reshape((nX, nY, 4))[:,:,0] # x-index of work group" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],\n", + " [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],\n", + " [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],\n", + " [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],\n", + " [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],\n", + " [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],\n", + " [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],\n", + " [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]], dtype=int32)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np_data.reshape((nX, nY, 4))[:,:,1] # y-index of work group" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "computeWorkgroup(8,2)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1]], dtype=int32)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np_data.reshape((nX, nY, 4))[:,:,0]" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],\n", + " [2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2],\n", + " [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],\n", + " [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3],\n", + " [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4],\n", + " [4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4],\n", + " [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5],\n", + " [5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5],\n", + " [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6],\n", + " [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6],\n", + " [7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7],\n", + " [7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7]], dtype=int32)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np_data.reshape((nX, nY, 4))[:,:,1]" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7],\n", + " [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7]], dtype=int32)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np_data.reshape((nX, nY, 4))[:,:,2]" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n", + " [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n", + " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]], dtype=int32)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np_data.reshape((nX, nY, 4))[:,:,3]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generate some visuals" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "def plotWorkgroupConfig(i,j):\n", + " computeWorkgroup(i,j)\n", + " data = np_data.reshape((nX, nY, 4))\n", + " for k in range(0,4):\n", + " plt.subplot(1,4,k+1)\n", + " plt.imshow(np_data.reshape((nX, nY, 4))[:,:,k])" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABrCAYAAABnlHmpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAACAxJREFUeJzt3U1oXNcZxvHn1YwlTZPKkj+jJnGahbNI7EKKiAktbSAtuN04m5Y4Gy0MLoFsuigVdBG6ateF0tYLY0NJQmkLcamhTQRNCi3GLoU2KSRWHCdWrVpRrUBiog9bbxeaiPlSR5q59545Z/4/CHfOyVHum0fXrw73zljm7gIAxG8gdAEAgGzQ0AEgETR0AEgEDR0AEkFDB4BE0NABIBE0dABIBA0dABLRVUM3s6Nm9paZzZjZVFZFYR355ods80O24VinnxQ1s5KktyV9XdKspIuSjrv7vzb7mkEb8mHdVf/f2VGuX1RuGEtaK9f/3PGyNa8pqWFN8/m9cU2pxf97i7lyea1uPDhwu2nNcKl+rjKw0rSmYqsN4+afp3/7x/KCu+/dbr6tskW9j7TYUbYS+bazpFta8WUj23x8eu22W9ei7W3ZY5Jm3P2KJJnZS5KOSdr0Gzesu3TEnqwvYM/+uvHavrGmr1vZU//NXtrdXPbSWH1zXB5rbvoro/XNenX0TtOa8s7mRrx79OO68YGRxaY1D909Xzc+VJltWnN46Hrd+JHBStOa0vjMe9WX28q3Vbao96r/uqNsJfJt54JPf/qSbHNQc+3+X93ccrlX0rWa8Wx1ro6ZnTSzS2Z2aVXLXZyu77TNl2w7xrWbH7INqJuG3rwFlpruV7j7KXefcPeJHRrq4nR9p22+ZNsxrt38kG1A3TT0WUn314zvk3R9k7XYPvLND9nmh2wD6uYe+kVJB83sQUn/lvS0pGcyqQrSNvO1HeWm5xFoMLfxims3P2QbUMcN3d1vm9lzkv4gqSTptLu/mVllfY5880O2+SHbsLrZocvdz0s6n1EtaEC++SHb/JBtOHxSFAAS0dUOHT2kXG75Hn7UmGu/BIgZO3QASAQNHQASQUMHgETQ0AEgETwUTcRaeaDpLzED0F/YoQNAImjoAJAIGjoAJIJ76InwsrX8xR8A+gc7dABIBA0dABJBQweARNDQASARPEVLxFpJWhrj5zPQz+gAAJAIGjoAJIKGDgCJ4B56IrwsLY9Z6DIABMQOHQASQUMHgETQ0AEgETR0AEgED0UT4SVpZdRDlwEgIHboAJAIGjoAJIKGDgCJ4B56IrzkWh29E7oMAAGxQweARNDQASARNHQASAQNHQASwUPRVJRc5Z0roasAEBA7dABIBA0dABLRtqGb2WkzmzezN2rmdpnZK2Z2uXocy7fMdJ347g3dc+hdfeGJ9zfmyDcb8z//ra5+50e69r2fbMyRbTbe9Et6zX+nv/ofN+bINryt7NDPSDraMDcladrdD0qaro7Rgclvj+j8C+ON0+Sbgc9+9VGNT002TpNtBj6nB/Sovtw4TbaBtX0o6u6vm9nnG6aPSXqi+vqspD9J+n6GdfWNrzxe0dVrq43T2863XF7T7tGPsy4vbo/v1fKND3WztJ7NlfVZrt0MjNlefeK3GqfJNrBO76Hvd/c5Saoe92VXEkS+eSLb/JBtYLm/bdHMTko6KUnD+kzep+srtdkO7hsJXE16uHbzQ7b56HSHfsPMxiWpepzfbKG7n3L3CXef2KGhDk/Xd7aUb2225Z38odgirt38kG1gne7Qz0malPTj6vHlzCqC1EG+gwO3dWBkMe+6onPr1ke6OnBHB0YWdXF9ims3P2QbWNuGbmYvav1Bxx4zm5X0vNa/Yb8ysxOS3pf0rTyLTNkzz/5Hr/3lEy3cvKMDX3xXkvaIfDNx4flXtfD3OS1/uKTzT/1SItvM/NMvaFEfaFXL+rP/XpJJZBvcVt7lcnyTf/VkxrX0pRd+dk/duDQ+s+Du/xX5du3ID79WN/7Nl35Bthk5bEfqxhd8Wkt+i2wD45OiAJAIGjoAJIK/bTERw6XbeujuTd9UAKAPsEMHgETQ0AEgETR0AEgE99ATURlY0aHKbOgyAATEDh0AEkFDB4BE0NABIBE0dABIBA9FE1GxVR0euh66DAABsUMHgETQ0AEgETR0AEgE99ATUbEBPTJYCV0GgIDYoQNAImjoAJAIGjoAJIKGDgCJMHcv7mRmH0h6T+u/fX2hsBNnp4i6H3D3vdv9IrLdko6ylaLPl2zz1TP5FtrQN05qdsndJwo/cZdiqDuGGluJpe5Y6qwVS82x1Nmol+rmlgsAJIKGDgCJCNXQTwU6b7diqDuGGluJpe5Y6qwVS82x1NmoZ+oOcg8dAJA9brkAQCIKb+hmdtTM3jKzGTObKvr8W2Fmp81s3szeqJnbZWavmNnl6nEsZI2txJCtFGe+ZJuvGPKNIdtCG7qZlST9VNI3JD0s6biZPVxkDVt0RtLRhrkpSdPuflDSdHXcMyLKVoosX7LNV0T5nlGPZ1v0Dv0xSTPufsXdVyS9JOlYwTW05e6vS7rZMH1M0tnq67OSniq0qPaiyFaKMl+yzVcU+caQbdEN/V5J12rGs9W5GOx39zlJqh73Ba6nUczZSr2dL9nmK+Z8eyrbohu6tZjjbTbZINv8kG2+yDcjRTf0WUn314zvkxTLbza+YWbjklQ9zgeup1HM2Uq9nS/Z5ivmfHsq26Ib+kVJB83sQTMblPS0pHMF19Cpc5Imq68nJb0csJZWYs5W6u18yTZfMefbW9m6e6H/SPqmpLclvSPpB0Wff4s1vihpTtKq1ncPJyTt1vpT7MvV467QdcaYbaz5ki35xpAtnxQFgETwSVEASAQNHQASQUMHgETQ0AEgETR0AEgEDR0AEkFDB4BE0NABIBH/A/6TJLQTZrhjAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plotWorkgroupConfig(1,1)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABrCAYAAABnlHmpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAB6NJREFUeJzt3U+IXeUZx/Hf00mMSWhL4vhn6p9WJIKiRcqQ4kYEkcRuRiiVRApZBETBTVcOdOHGRVYuBBFDCclGpbhJoAGxgdZFKWRKS0gLJqOJmakxcUz8A5V0Rh4Xcw03Jzfek3POe859n/v9QLj3nNx5z8NvXh4O58y5r7m7AAD5+0HXBQAAmkFDB4AgaOgAEAQNHQCCoKEDQBA0dAAIgoYOAEHQ0AEgiFoN3cy2m9n7ZjZvZrNNFYVV5JsO2aZDtt2xqk+KmtmEpBOSHpe0KOmopJ3u/p9r/cwNts5v1MZKx+u3Mll/DEmyTSuNjCNJ962/2Mg4/zh2acndb77efJvKNrKvdLFSttLgfIvzcNB8Ks6LE8c2XPWZVOMMGmvQPC2OVWWc0wvLWrrwjXXZFyL7bu4O+9yaGsfYKmne3T+UJDN7S9KMpGv+4m7URv3SHqtxyFVLv3649hiStPbJTxsZR5L+/tDbjYwzMTX/Ue/tdeXbVLaR/dnfrpStNDjf4jwcNJ+K82LbTx666jOpxhk01qB5Whyryjhbty1cfquO+kJkfXP3e9W55HK7pIW+7cXeviuY2TNmNmdmc8u6VONwY2dovmRbGXM3HbLtUJ2GbgP2XXX9xt33uvu0u0+v1boahxs7Q/Ml28qYu+mQbYfqNPRFSXf2bd8h6eN65aAP+aZDtumQbYfqXEM/KmmLmd0t6b+Sdkh6upGqIJFvSmSbznVnuzK5sbH7YmG9Xu4eXeWG7u4rZva8pHckTUja5+7/rjoerkS+6ZBtOmTbrTpn6HL3w5ION1QLCsg3HbJNh2y7w5OiABBErTN0AKjLNq00+kxISK+X+xhn6AAQBA0dAIKgoQNAEDR0AAiCm6IAOnXf+ouNfbldVBMlP8cZOgAEQUMHgCBo6AAQBNfQg+ALjkoo+QVHaNeJYxsGLtiBfvOlPsUZOgAEQUMHgCBo6AAQBA0dAILgpiiATnFDv4SSN/Q5QweAIGjoABAEDR0AguAaOoBOsWJRCaxYBADjhYYOAEHQ0AEgCBo6AATBTVEAnWLFouFYsQgAxgwNHQCCoKEDQBBcQw+ChzNKKPlwBtrFikVlsGIRAIwVGjoABEFDB4AgaOgAEAQ3RQF0ihWLSmDFIgAYLzR0AAhiaEM3s31mdt7Mjvft22xm75rZyd7rprRlxrX7d+d02wOn9PNHz1zeR77NOPXyn/TPHa/o+LN/uLyPbJvBvB1NZc7Q90vaXtg3K+mIu2+RdKS3jQp2PfUjHX5jqribfBsw+fiDuvelp4q7ybYBzNvRNPSmqLu/Z2Y/K+yekfRo7/0BSX+R9EKDdY2NRx5er9MLy8Xd5NuAHz54ly6d+7y4m2wb0OS85SnnEhIvQXeru5+VpN7rLRXHwWDkmw7ZpkO2HUt+U9TMnjGzOTObW9al1IcbK/3Zrnzxv67LCYe5mw5zN42qDf2cmU1JUu/1/LU+6O573X3a3afXal3Fw42dUvn2Z7vmxxtaLTBjzN10KmXL3G1O1QeLDknaJWlP7/VgYxVBIt+UyDadStmyYtFwZVcsGtrQzexNrd7omDSzRUkvavUX9kcz2y3pjKTfVKxz7D393Cf669++1tKFb3TXL05J0qTItxEf7Dmor46d0cqXX+tfv31VItvGFOftmgmTyLZzZf7KZec1/uuxhmsZS2+8dtsV2xNT80vu/pnIt7Z7Zmeu2D76xB6ybUhx3m7dtqDTC8tk2zGeFAWAIGjoABAE37YYBDeWhit7YwntYgm6MliCDgDGCg0dAIKgoQNAEFxDB9ApViwqgRWLAGC80NABIAgaOgAEQUMHgCC4KQqgU6xYVELiFYsAACOGhg4AQdDQASAIrqED6BRfLDdc2S+W4wwdAIKgoQNAEDR0AAiChg4AQZi7t3cws08lfaTV1deXWjtwc9qo+6fufvP1/hDZllIpWyn7fMk2rZHJt9WGfvmgZnPuPt36gWvKoe4cahwkl7pzqbNfLjXnUmfRKNXNJRcACIKGDgBBdNXQ93Z03LpyqDuHGgfJpe5c6uyXS8251Fk0MnV3cg0dANA8LrkAQBCtN3Qz225m75vZvJnNtn38Msxsn5mdN7Pjffs2m9m7Znay97qpyxoHySFbKc98yTatHPLNIdtWG7qZTUh6VdITku6XtNPM7m+zhpL2S9pe2Dcr6Yi7b5F0pLc9MjLKVsosX7JNK6N892vEs237DH2rpHl3/9Dd/y/pLUkzLdcwlLu/J+lCYfeMpAO99wckPdlqUcNlka2UZb5km1YW+eaQbdsN/XZJC33bi719ObjV3c9KUu/1lo7rKco5W2m08yXbtHLOd6Sybbuh24B9/JlNM8g2HbJNi3wb0nZDX5R0Z9/2HZI+brmGqs6Z2ZQk9V7Pd1xPUc7ZSqOdL9mmlXO+I5Vt2w39qKQtZna3md0gaYekQy3XUNUhSbt673dJOthhLYPknK002vmSbVo55zta2bp7q/8k/UrSCUkfSPp928cvWeObks5KWtbq2cNuSTdp9S72yd7r5q7rzDHbXPMlW/LNIVueFAWAIHhSFACCoKEDQBA0dAAIgoYOAEHQ0AEgCBo6AARBQweAIGjoABDEt1uHTloE0FJaAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plotWorkgroupConfig(4,4)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABrCAYAAABnlHmpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAB9JJREFUeJzt3c9rHPcdxvHn41WyqxgEkexEyo/WwTgEUUwowqWXEgjFbi/upSbOwSYYAoVeeqqhh177D5SCD8a+JKH0Eh8cRGpoc+hFKgkmLcRWEztS7dpR3ITg1K4kPj1oa3ZnR9Fqdma+8/36/QIzO+OR5uHx+sPszK5k7i4AQPx2hQ4AACgHAx0AEsFAB4BEMNABIBEMdABIBAMdABLBQAeARDDQASARIw10MztiZh+Z2ZKZnS4rFDbRb3Xotjp0G44V/aSombUkXZH0Q0krkhYkHXf3v2/1NY9a2zvaXeh4VXj+4NehIwz46+X7q+6+d6f9Nq3bJvpK/y7UrZTfr413+tbXx1sDX7fR7l/39uD/t057rW99cuzuwD5TrfW+9SuXHxvYJ5snL1M2T16mbJ68TNk815bXtHpnw1KYC030/+fudvuNjXCMQ5KW3P1jSTKztyQdlbTlP1xHu/U9e3mEQ5Zrfv6D0BEGtGaWrncf7qjfpnXbRH/0PxTqVsrvd9cLs33rX8xODHzdl/v7XwTf239/YJ/ZfTf61o9NLwzsc2JitW/98FMvDuyTzZOXKZsnL1M2T16mbJ5Dh5cfPFTkc6GJep6732iUSy5PS1ruWV/pbutjZq+b2aKZLa5p8MmMLW3bL90WxnO3OnQb0CgD3XK2DbyedPcz7j7n7nOPKOf1Hraybb90WxjP3erQbUCjDPQVSc/2rD8jafC1Goqi3+rQbXXoNqBRrqEvSDpgZs9J+qekVyS9WkoqSPRbJbqtzo67tfFO7vV/9Hh/uN0KD3R3Xzezn0ual9SSdNbd/1b0+6Ef/VaHbqtDt2GNcoYud78o6WJJWZBBv9Wh2+rQbTh8UhQAEjHSGToAjGp9vJX7Hn70GPIaOmfoAJAIBjoAJIKBDgCJYKADQCK4KQogqI12/g8Nw87RIgAkgoEOAIlgoANAIriGDiAob3vuL/7AznGGDgCJYKADQCIY6ACQCAY6ACSCm6IAguq01zS7j99S902uD7kfZ+gAkAgGOgAkgoEOAIngGjqAoCbH7urY9ELoGI32zpD7cYYOAIlgoANAIhjoAJAIBjoAJIKbogCCmmqt68TEaugYjfbakPtxhg4AiWCgA0AiGOgAkAiuoSfi+YNfa37+g9AxGq01EzoB8ly5/JgOP/Vi6BgNtzTUXpyhA0AiGOgAkAgGOgAkgoEOAIngpiiAoGy8o10vzIaO0WzvD7cbZ+gAkAgGOgAkYtuBbmZnzey2mX3Ys23SzN41s6vd5ePVxkzXqV/c0vR3PtHBlz59sI1+y0G31aHbZhrmDP2cpCOZbaclXXL3A5IudddRwMljE7r4xsAnXui3BHRbHbptpm1virr7e2a2L7P5qKSXuo/PS/qTpF+WmOuh8YPvj+va8lp2M/2WgG6rU2a36+MtfTE7UWa89FR8U/RJd78pSd3lEwW/D/LRb3Xotjp0G1jlN0XN7HUzWzSzxTXdr/pwD5Xebj/7fCN0nOTw3K1Ob7fr9+6GjpOMogP9lpnNSFJ3eXurHd39jLvPufvcI2oXPNxDZ6h+e7vdO9WqNWDEeO5Wp1C3Y53dtQVMXdEPFl2QdFLSb7rLt0tLBIl+q0S31SnU7UZb+nI/76Auw7YD3cze1OaNjj1mtiLp19r8B/u9mZ2S9Kmkn1YZMmWv/uxf+vNf/qPVOxv61nc/kaQ9ot9S0G11st2OtUyi2+CGeZfL8S3+6uWSszyU3vjddN96a2Zp1d0/F/2OjG6rk+320OFlXVteo9vAeJ0DAIlgoANAIvhpiwCC8rbr3n7eFloGztABIBEMdABIBAMdABLBNXQAQXXaa5rddyN0jEa7PuR+nKEDQCIY6ACQCAY6ACSCgQ4AieCmKICgJsfu6tj0QugYjfbOkPtxhg4AiWCgA0AiGOgAkAiuoQMIaqq1rhMTq6FjNNprQ+7HGToAJIKBDgCJYKADQCIY6ACQCHP3+g5m9pk2f3DYHkkx3gWpI/e33X3vTr+IbodSqFsp+n7ptlqN6bfWgf7goGaL7j5X+4FHFEPuGDLmiSV3LDl7xZI5lpxZTcrNJRcASAQDHQASEWqgnwl03FHFkDuGjHliyR1Lzl6xZI4lZ1Zjcge5hg4AKB+XXAAgEbUPdDM7YmYfmdmSmZ2u+/jDMLOzZnbbzD7s2TZpZu+a2dXu8vGQGfPE0K0UZ790W60Y+o2h21oHupm1JP1W0o8kzUo6bmazdWYY0jlJRzLbTku65O4HJF3qrjdGRN1KkfVLt9WKqN9zani3dZ+hH5K05O4fu/t/Jb0l6WjNGbbl7u9JupPZfFTS+e7j85J+Umuo7UXRrRRlv3RbrSj6jaHbugf605KWe9ZXutti8KS735Sk7vKJwHmyYu5Wana/dFutmPttVLd1D3TL2cbbbMpBt9Wh22rRb0nqHugrkp7tWX9G0o2aMxR1y8xmJKm7vB04T1bM3UrN7pduqxVzv43qtu6BviDpgJk9Z2aPSnpF0oWaMxR1QdLJ7uOTkt4OmCVPzN1Kze6XbqsVc7/N6tbda/0j6ceSrkj6h6Rf1X38ITO+KemmpDVtnj2ckjSlzbvYV7vLydA5Y+w21n7pln5j6JZPigJAIvikKAAkgoEOAIlgoANAIhjoAJAIBjoAJIKBDgCJYKADQCIY6ACQiP8BpLw4oSzcg78AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plotWorkgroupConfig(8,8)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABrCAYAAABnlHmpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAACCdJREFUeJzt3V+IXGcZx/Hfk5nsn7bGbP6nf1NKiqxtqbqkCEULRRK9iRdaml4YSiAg9EYQDHpRBC+8FESEICG5aYvkphFTYxuxFRRJxKCt0GRtk+5202zXxKptk90kjxc7LnNmZp3ZM+ff+873A+HM++bsnie/Pfvk5Jx5s+buAgCEb1XZBQAAskFDB4BI0NABIBI0dACIBA0dACJBQweASNDQASASNHQAiERfDd3MdpnZm2Y2aWYHsioKi8g3P2SbH7Itj6VdKWpmNUlnJX1J0rSkU5L2uPvflvuYIRv2Ed2a6niD4t+6MufuG1eaL9l2lzZbqXO+trqe3KneMpZ0s568ZvK6te9TU8s+7cf31n1qHb5vO8zV6zcT46FV19v2Gakl50ZXzbftM2oLLePkn+v81ILmLt8w+kI+/nfudtuvw6nTsx2SJt39LUkysxck7Za07BduRLfqEXu8j0PG7xU/eqHxckX5km13abOVOudb37A5Mb65aazt4+Y3JBvV1fXt33JXx5LN8dpYe9OfX5ts1gtrb7TtU/9keyNev/Y/ifHda6607XP/bbOJ8QOj0237PDg8kxh/emg0Md6xc2rppegLmWs6d/+vfm653CFpqmk83ZhLMLP9ZnbazE4v6Fofhxs4XfMl29Q4d/NDtiXqp6G3X0ZIbf/mc/eD7j7h7hOrNdzH4QZO13zJNjXO3fyQbYn6aejTku5qGt8paWaZfbFy5Jsfss0P2Zaon3vopyRtN7N7Jb0r6UlJT2VSFSTyzRPZ5mfF2d7/0Ec6ceJMEbUFq7a1t/1SN3R3v25mz0g6Iakm6ZC7v5H28yGJfPNDtvkh23L1c4Uudz8u6XhGtaAF+eaHbPNDtuVhpSgARKKvK3RUh42OaNWnxssuo9r+XHYB6OTsX27RztsfLruMipvsaS+u0AEgEjR0AIgEDR0AIkFDB4BI8FAUQKlYWNRdrwuLuEIHgEjQ0AEgEjR0AIgE99AjcX20pn+Orym7jGpjYVElsbCoFywsAoCBQkMHgEjQ0AEgEjR0AIgED0UBlIqFRd2xsAgABgwNHQAiQUMHgEhwDz0SN4alD+7j72eEh4VFvWBhEQAMFBo6AESChg4AkaChA0AkeCgKoFQsLOqOhUUAMGBo6AAQCRo6AESCe+iR8GHX1fuulV0GsGIsLOoFC4sAYKDQ0AEgEjR0AIgEDR0AIsFDUQClYmFRdywsAoABQ0MHgEh0behmdsjMZs3s9aa5dWb2spmda2zH8i0zXm/4ab3qv9Af/NdLc+SbjbmfHdXUMz/QzHd/tDRHttnY961L2vLA23rosXeW5si2fL1coR+WtKtl7oCkk+6+XdLJxhgp3K579Bk92jpNvhm47dHPadO3n26dJtsM7H1ijY4/13Zjl2xL1vWhqLu/ZmbbWqZ3S3qs8fqIpN9K+k6GdQ2MMduoj/3D1ukV5zsyvKDxbTNZlxe2bcP66OK/9KehxWwuLM5y7mbgC58f1fmphdbpVNmyUrQX+a4U3ezuFyWpsd2U8vOgM/LND9nmh2xLlvvbFs1sv6T9kjSiW/I+3EBJZLv5EyVXEx/O3fyQbT7SXqFfMrOtktTYzi63o7sfdPcJd59YreGUhxs4PeXbnO3Q2tFCCwwY525+yLZkaa/Qj0naK+mHje2LmVUEiXzzRLb5SZUtC4u663VhUdeGbmbPa/FBxwYzm5b0rBa/YD83s32S3pH09bSFDrq/+h91Re9rQdf0O/+lJG0Q+WbizPd/pctn3tX8B1f1m68dksg2M0998z29+vuPNXf5hu7+7Nuq10wi29L18i6XPcv81uMZ1zKQHrRHEuNX/Oicu/9D5Nu3h59Nvtv2pS/+mGwz8txPtyTGO3ZO6fzUAtmWjJWiABAJGjoARIL/bTES6+of6oktp8ouo9JeKrsAdMTCol7wI+gAYKDQ0AEgEjR0AIgE99ABlIqFRd3xE4sAYMDQ0AEgEjR0AIgEDR0AIsFD0Uisr13XN9bMlV1GpbX9MDpUAguLesHCIgAYKDR0AIgEDR0AIsE9dAClYmFRdywsAoABQ0MHgEjQ0AEgEjR0AIiEuXtxBzN7X9IFLf709RBXwRRR9z3uvnGlH0S2PUmVrRR8vmSbr8rkW2hDXzqo2Wl3nyj8wH0Koe4QauwklLpDqbNZKDWHUmerKtXNLRcAiAQNHQAiUVZDP1jScfsVQt0h1NhJKHWHUmezUGoOpc5Wlam7lHvoAIDsccsFACJReEM3s11m9qaZTZrZgaKP3wszO2Rms2b2etPcOjN72czONbZjZdbYSQjZSmHmS7b5CiHfELIttKGbWU3STyR9WdK4pD1mNl5kDT06LGlXy9wBSSfdfbukk41xZQSUrRRYvmSbr4DyPayKZ1v0FfoOSZPu/pa7z0t6QdLugmvoyt1fk3S5ZXq3pCON10ckfbXQoroLIlspyHzJNl9B5BtCtkU39DskTTWNpxtzIdjs7hclqbHdVHI9rULOVqp2vmSbr5DzrVS2RTd06zDH22yyQbb5Idt8kW9Gim7o05LuahrfKWmm4BrSumRmWyWpsZ0tuZ5WIWcrVTtfss1XyPlWKtuiG/opSdvN7F4zG5L0pKRjBdeQ1jFJexuv90p6scRaOgk5W6na+ZJtvkLOt1rZunuhvyR9RdJZSX+X9L2ij99jjc9LuihpQYtXD/skrdfiU+xzje26susMMdtQ8yVb8g0hW1aKAkAkWCkKAJGgoQNAJGjoABAJGjoARIKGDgCRoKEDQCRo6AAQCRo6AETiv+qQRIYl3vsZAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plotWorkgroupConfig(16,2)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAABrCAYAAABnlHmpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAB8RJREFUeJzt3c9rHOcdx/HPV7vWj6Qoln8rv0NwKGoS2iIcCqUEQrHbi3tpiHOoDwZDoZeeauih1/4DpeCDsS9JKLnEpQY3MTQptBS5UNKkEFt1nEiNY0e1W5I0tlb2twdtzGq98oxmZ+bZ59n3C8zsPH525+uPhq8ez+xK5u4CAMRvJHQBAIBy0NABIBE0dABIBA0dABJBQweARNDQASARNHQASAQNHQAS0VdDN7N9Zvaemc2b2ZGyisIq8q0O2VaHbMOxop8UNbOGpHOSvitpUdKcpAPu/o/1njNqYz6ue+/+uhPjmcdemWhkzrk5ljlFPpbv3z4+1sqcs6X5eeacrY2VzDl/ffvGkrtv32i+ebIddp/qWqFspd75PvH0/9bsn3v7njueF3JOr3lVzbm40NLS1ZtWVV8Ydl+eu1nzmn0cY4+keXe/IElm9oqk/ZLW/cKN6149Y8/d9UVHvjqTeeD/zExmzvnv49n/+bj++I3MOZI08+hHmXOe3zWXOedHk0uZcxrT8x+0H24o3zzZDrs3/NVC2Uq98z19+m9r9vfe//U7nhdyTq95Vc3Zs3fh9kNV0BeGXce5e1f9XHJ5QNJCx/5ie2wNMztsZmfN7GxL+RooJOXIl2wL49ytDtkG1E9Dtx5jd1zDcPej7j7r7rOblOM6CL6UmS/ZFsa5Wx2yDaifhr4o6aGO/QclZV+bQF7kWx2yrQ7ZBtTPNfQ5SbvN7DFJ/5L0gqQXS6kKEvlWiWyrs+FsbVNTzW0766gtXpfyTSvc0N19xcx+Ium0pIakY+7+btHXw1rkWx2yrQ7ZhtXPCl3ufkrSqZJqQRfyrQ7ZVodsw+GTogCQiL5W6ADQt2ZTt3ZMha5isOW8hs4KHQASQUMHgETQ0AEgETR0AEgEN0UBBHWrOaLlbfy0xTKwQgeARNDQASARNHQASATX0AEE5U3T9a20ojKwQgeARNDQASARNHQASAQNHQASwZ0IAEHdakjXp1hbloEUASARNHQASAQNHQASwTV0AEF5U7oxZaHLSAIrdABIBA0dABJBQweARNDQASAR3BQFEJQ3pOXNHrqMJLBCB4BE0NABIBE0dABIBNfQAQTlDVdr883QZSSBFToAJIKGDgCJoKEDQCJo6ACQCG6KAgir4Wretxy6iiSwQgeARNDQASARmQ3dzI6Z2RUze6djbIuZvW5m59vbqWrLTNehn17Wriff19PPfnh7jHzL8a6f1Zv+W/3Zf397jGzLwXk7mPKs0I9L2tc1dkTSGXffLelMex8FHHx+Uqdemu4eJt8S3K9H9A19u3uYbEvAeTuYMm+KuvtbZvZo1/B+Sc+2H5+Q9AdJPyuxrqHxnW9N6OJCq3uYfEswZdv1hX/ePUy2JSjzvG02b2nr5s/KLC85F3LOK3oNfae7X5Kk9nZHwddBb+RbHbKtDtkGVvnbFs3ssKTDkjSue6o+3FAh22qRb3U6sx3dMRm4mnQUXaFfNrNpSWpvr6w30d2Puvusu89u0ljBww2dXPmSbSGcu9UplG3zPr5ZlqXoCv2kpIOSftnevlZaRZDIt0pkW51C2Y6OrOjhyWtV1hW9uZzzMhu6mb2s1Rsd28xsUdIvtPoF+42ZHZL0oaQfFqxz6L3444/15p++0NLVm3r4m+9L0jaRbyn+7n/RNX2ilm7oj/47iWxL033eNhsmkW1wed7lcmCdv3qu5FqG0ku/3rVmvzE9v+Tu/xb59u0pe2bN/hv+KtmWpPu83bN3QRcXWmQbGJ8UBYBE0NABIBH8tEUAQY03VvTEV9Z9Qww2gBU6ACSChg4AiaChA0AiuIYOIKiJkWU9ObEYuowksEIHgETQ0AEgETR0AEgEDR0AEsFNUQBBTVhLT419FLqMJLBCB4BE0NABIBE0dABIBNfQAQQ1YSP62uhE6DKSwAodABJBQweARNDQASARNHQASIS5e30HM/tE0gda/e3rS7UduDx11P2Iu2/f6JPINpdC2UrR50u21RqYfGtt6LcPanbW3WdrP3CfYqg7hhp7iaXuWOrsFEvNsdTZbZDq5pILACSChg4AiQjV0I8GOm6/Yqg7hhp7iaXuWOrsFEvNsdTZbWDqDnINHQBQPi65AEAiam/oZrbPzN4zs3kzO1L38fMws2NmdsXM3ukY22Jmr5vZ+fZ2KmSNvcSQrRRnvmRbrRjyjSHbWhu6mTUk/UrS9yTNSDpgZjN11pDTcUn7usaOSDrj7rslnWnvD4yIspUiy5dsqxVRvsc14NnWvULfI2ne3S+4+7KkVyTtr7mGTO7+lqSrXcP7JZ1oPz4h6Qe1FpUtimylKPMl22pFkW8M2dbd0B+QtNCxv9gei8FOd78kSe3tjsD1dIs5W2mw8yXbasWc70BlW3dDtx5jvM2mHGRbHbKtFvmWpO6GvijpoY79ByXF8tthL5vZtCS1t1cC19Mt5mylwc6XbKsVc74DlW3dDX1O0m4ze8zMRiW9IOlkzTUUdVLSwfbjg5JeC1hLLzFnKw12vmRbrZjzHaxs3b3WP5K+L+mcpH9K+nndx89Z48uSLklqaXX1cEjSVq3exT7f3m4JXWeM2caaL9mSbwzZ8klRAEgEnxQFgETQ0AEgETR0AEgEDR0AEkFDB4BE0NABIBE0dABIBA0dABLxf0GZN9fdsbfGAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plotWorkgroupConfig(2,16)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebook/lbm_codegen.ipynb b/notebook/lbm_codegen.ipynb new file mode 100644 index 0000000..e1593c0 --- /dev/null +++ b/notebook/lbm_codegen.ipynb @@ -0,0 +1,828 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import *\n", + "init_printing()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Characteristic constants" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "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": { + "text/latex": [ + "$$\\left [ \\left[\\begin{matrix}-1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\-1\\end{matrix}\\right]\\right ]$$" + ], + "text/plain": [ + "⎡⎡-1⎤ ⎡0⎤ ⎡1⎤ ⎡-1⎤ ⎡0⎤ ⎡1⎤ ⎡-1⎤ ⎡0 ⎤ ⎡1 ⎤⎤\n", + "⎢⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥⎥\n", + "⎣⎣1 ⎦ ⎣1⎦ ⎣1⎦ ⎣0 ⎦ ⎣0⎦ ⎣0⎦ ⎣-1⎦ ⎣-1⎦ ⎣-1⎦⎦" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "c" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "w = [Rational(*x) for x in [(1,36), (1,9), (1,36), (1,9), (4,9), (1,9), (1,36), (1,9), (1,36)]]" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAAVCAYAAABi32uwAAAABHNCSVQICAgIfAhkiAAABzlJREFUeJztnWuoFVUUx3++wqKHWZEEkhZC4Y1I0jJMt0VShqFSXwIhpPoSFIhF9KBjlFlRaFpEVJr5IeghfQgrE063J0pZUWkvtRdpXtOyNM1HH9Ye7pxhzjkzc/acM/vM+sFh7t2zZ85ea/9n1r77dUFRFEVRFEVxggGOhj6bO1oaRVEURVGU4nAqte2ko8GJgTGZ3wUWAMsi6dcCS4H3gL/sTVYlLMDpwGHgiVDaw8A64GdgP/AHsBG4Dzilwb0uBV4FfgMO2OPbwPSEZWkHLn01AJgLfAzsBfYhfroVGOSuyE5waXdWfUDxNeL6WfJJI2Hm0P9CurFJXp/tLmt9a8xIjsaMdHa71Ack18g+pH20APix3s2MNaRS5/xn9vxeYBPpKvtmm39qKO0gUtnPA4sQh26w+X4FRsbc5x57fiewHFgIPGOveyRhWdqBS1+ttL/vAJ4FlgBf2bRXkAenKLi0O4s+wA+NuH6WfNJIwEhgD+KDJI0pn+0ua31rzEiOxox0drvSB2TXSJVQz1QYQ+PG1FRgDFIRQd6klb0G6KO2VTy0Tt4H7b2fiqRfZ9PXAifEXDckYVnagStfzbTXbkG6FwOGAKvtuRtcFNgRLjWSVh/gj0Zc+sk3jYDY/Q7wA/AoyRpTPttd1vrWmJEcjRnp7HahD2hNI1UyNqbi8iYx+iSk22x5grwA59NvXMBARBz/AKclvE9RMGT3VfAXxi0xeXvsuU9aL2IuGPLRSJw+wF+NGFrzk48auQ04AkxG3jfNGlPdYjeUs75BY0YaDBozGuFCH9C6RqqEGlODM9wgLVcDxwCvJcw/wx6/CKVdAoxGuih323v2AP8C64GPnJS080R9NcIet8TkDdLGAcOQIRNfSaOROH1AOTQS5yffNHIu0gW/BOgFLktwTTfYnQW1OxkaM8oVM1zoAxxrpB2NqVlIyy/aKgyYDxyPtDYvBCYhRi8K5RlvjzuAT4HzIvfoRSay7XRT5I4R9VWfPY6OyXtW6OdzkLFiX2mkkST6gHJoJM5PPmlkMPAi8BNwV4rrfLc7K2p3PBoz+iljzHChD8hRIwb3w3xDkYllLzfIs53aZYZrkJn6YR6y5w4B3wGXI84aC7xpz1UTlLsTGLL76np77ffA8FD6YGTlQeCzqxyV1SUGNxpJog/wVyOG1vzkk0buR1bgTAylVWg8zNcNdocxlKe+wxg0ZiTFoDGjHq70Aa1rpEpomC9uawSXTEMKt7pBnhHI5LMRwGyk9bwR6YoMCCaZDUBaiuuAv5EVCrOAX4Ap1L6kfSPOVy8hQjgb+BpZYbAYWQExHREASIDylWYaSaIP6H6N1POTLxqZgPRGPUa67nPf7c6K2l0fjRlCGWOGK32AY43k3ZiahSxXfCNB3h2Ig6Yhe0KsDJ3bbY9bgM8j1+0H3rI/T8hc0s4T56sjwDVIt+V2ZF+euUglTwJ22Xy/t6+YzkmqkUb6gO7XSD0/+aCRYHjvW+DelNf6bHcrqN3N0ZhRvpjhSh+Qo0YMbof5BiHjt2vSFgRpRR6lf2nnbPv7hjr5g+XVd2b4rrwx5OOrY5EK30exlvgGGPLTSFQf4K9GDPn5qSgaGUZtt3ujz+LQdb7bHYeh++s7DoPGjKQYNGbE4VIf0LpGqrRpNd9kpDXYqDuuHmfYY9AV2YuMa45BZvEfjOTvscdtGb6rCGTx1Rxk/PgF4L88CtUGsmokqg/obo1k9VNRNHIAeK7OuXHABcD7wDfUDgH6bndW1O70aMxojuqjTTHD4LZnailS8LiJX+fQv4QzzED6N9j6IHJulU1/IJJ+BdK1uQf5CzjMCjq/SZmhNV+dGJM2Htkqfy+1KzQCVuC33Vn0AX5qxNCaPsBfjUDjCejdaLeh/fW9gs7XtUFjRlIMGjPicK0PyKaRgCoZe6Zm2g/0F3oiUgkg3W/zI/k/RMYto1yJdKH1Irsg70IcNAWp6O3ATZFr5gEXAXcjLdT1wJnIGOphmz+6b0YwJ+xQc/Oc4tJXa5Gu2S+RB2EsMpHwANJNGbefiO92Z9EH+KMRl/oAvzSShm6xu9P17fv7ADRmaMxwrw/IppGmGBr3TFVoPNdhWyjveJs2r869eoAnkRUGfUgF/omMXVaoXdIZZjjwOLAV6ZLbBbwOXFwn/0bknyeeXOd8XlRw56vbkR1r9yAPw1bgaWBUg+/33e6s+gA/NFLBnT7AL41EqRDfM9VNdlfobH37/j4AjRnbQnnLGDPy0gek10hAlVDPVBhD8mG+Ziy094rbOKxdDENalkX6Z5ZxuPZVWe3Ogg++ysNPandx0fdB5yirr3ywuwj6iFKlSWMq+Gxu4Us2IS3ETjID2RY+bhy1SLj2VVntzoIPvsrDT2p3cdH3Qecoq698sLsI+gBZDRjtPQNks6qAUdROPusDluVfNkVRFEVRlMJzHHBHJK3SgXIoiqIoiqIoiqIoiqIoiuV/+NxET1RKZ7IAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\left [ \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{4}{9}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}, \\quad \\frac{1}{9}, \\quad \\frac{1}{36}\\right ]$$" + ], + "text/plain": [ + "[1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "w" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAHZJREFUKJFjYKACCGFgYJjMwMBwmIGB4RMDA8N/BgaGJYQ0XYAq/MzAwHCdWE2ODAwMqgwMDIwMDAwOuDSxoPH3EzKVgYGBgYkYRaOaBlwTeuQGQDEDAwODBJS2ZGBgWABlv2FgYChBN6SBAZJ0cOEH5LiMzgAA6XoX52TB9a4AAAAASUVORK5CYII=\n", + "text/latex": [ + "$$1$$" + ], + "text/plain": [ + "1" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sum(w)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAADEAAAAdCAYAAAAD8oRRAAAABHNCSVQICAgIfAhkiAAAAi9JREFUWIXt1t2LTVEYx/GPGTTFNJkhk5SoKeRCopAkaoSkhD/BtXArLrwkF8pbjQvuKcqUuZIkbpTxFlFKTSF5yUthzMTF2jNzzmnvffY5Z585Ls63dmu113p+63n2evazFk3+D6ZUOL8d3xq4fi7swuZGLJwn5zC90U6U0lLB3NboGa6TL1VTSRBrcb9ejkwWJ9DVaCfiqGQnuvCpXo7UQtYgevAyYewkbmEIP/EZgzgs+87NxSjOVKObtU4fQD9exYwN4yGe4wNmYDVW4m3UHyqjvxd92IjbOeoW0Zcy1pbw/hj+4kIG/QF8FKpfxbqF6dSCo9hSYtQp/V/4lfD+StT2pNhCh7AD/UJKVaxbGEQv1mB3idFW3CzjSBzbo/ZJmXnbhAP0Wl66y4X8KwzuvGwF4CCO4DTuClv+GHPK2F3FD8npU5XuG6yL+tMUV4w03kcLjD0DQtVJow3fhUCq1o37wjewI+pvwJ0yjozRLVS7buzEIqEkrkix6cVMXM9Z1yYTpfRUtEg1LMBvPEuZczma05Gzrqn4giWylcc0BoUUmB0z1iqU1YFadePSaUSoRofwtIoFCpkXtaMxY+uFkzctlarRHWePEOn8MmKLhVwtpcXEoXQvwfZs5ETcz1+L7jjtsl279+GPcMe5KNx0L+F1tNA7LE2wHRJKZt66RczKMGeZcI48EvJ7BF/xQKjtnQl2qyJn9uesO6kcF4JY2GhHauGF8JWbNGlSB/4BwBaSc+fGZlwAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\frac{\\sqrt{3}}{3}$$" + ], + "text/plain": [ + "√3\n", + "──\n", + "3 " + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "c_s = sqrt(Rational(1,3))\n", + "c_s" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Moments" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "rho, tau = symbols('rho tau')" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([f_next_0, f_next_1, f_next_2, f_next_3, f_next_4, f_next_5,\n", + " f_next_6, f_next_7, f_next_8], dtype=object)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f_next = symarray('f_next', q)\n", + "f_next" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([f_curr_0, f_curr_1, f_curr_2, f_curr_3, f_curr_4, f_curr_5,\n", + " f_curr_6, f_curr_7, f_curr_8], dtype=object)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f_curr = symarray('f_curr', q)\n", + "f_curr" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$$\\left[\\begin{matrix}u_{0}\\\\u_{1}\\end{matrix}\\right]$$" + ], + "text/plain": [ + "⎡u₀⎤\n", + "⎢ ⎥\n", + "⎣u₁⎦" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u = Matrix(symarray('u', d))\n", + "u" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy.codegen.ast import Assignment" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "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 [ 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/jPLyPOxGHXGDLRDuAvt2FxUo0eLDabc1nVDz3QiMmS3oQFSDX1o298FGRbQ0WAf40i21A3Az2kdTgjJnIIG6DB01LYf4eOuId6Z6IjfGSj0cTrylkHGb