02_iteration-domains.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polyhedral Representation\n",
    "\n",
    "将一个命令式语言的statement作为identifer, 加上loop的约束, 将其转化为Union Set \n",
    "\n",
    "| Imperative Program | | Iteration Domain |\n",
    "|-|-|-|\n",
    "| <pre><code class=\"c language-c\">const int N = 100; <br>double X[N], Y[N], Z[2*N]; <br>  for (int i = 0; i &lt;= 2*N; ++i) <br>S:  Z[i] = 0.; <br>    for (int i = 0; i &lt;= N; ++i) <br>      for (int j = 0; j &lt;= N; ++j) <br>T:      Z[i + j] += A[i] * B[j]; <br> </code></pre> | $ \\to $ | $( \\{ T[i,j] \\mid 0 < i \\le j < 10; S[i,0] \\mid 0 < i \\le 20 \\})$  |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statement instances\n",
    "\n",
    "考虑以下计算多项式乘积的代码片段\n",
    "每个多项式都由其系数数组表示\n",
    "\n",
    "```c\n",
    "const int N = 100;\n",
    "double X[N], Y[N], Z[2*N];\n",
    "\n",
    "for (int i = 0; i <= 2*N; ++i)\n",
    "S:  Z[i] = 0.;\n",
    "for (int i = 0; i <= N; ++i)\n",
    "    for (int j = 0; j <= N; ++j)\n",
    "T:      Z[i + j] += A[i] * B[j];\n",
    "```\n",
    "\n",
    "为其中一些statement添加label作为identifier.\n",
    "\n",
    "Statement `S` 初始化数组`Z`的元素,statement `T`计算它们.\n",
    "Statement `S` 被包含在循环中,他将会按以下顺序被执行 $((2\\mathtt{N} + 1))$ 次:\n",
    "\n",
    "* `Z[0] = 0.; /* i = 0 */`\n",
    "* `Z[1] = 0.; /* i = 1 */`\n",
    "* ...\n",
    "* `Z[2*N] = 0. /* i = 2*N */;`\n",
    "\n",
    "把循环中每个单独执行的statement称为 **statement instances**\n",
    "这样每个instance可以通过语句label和封闭循环迭代器的值来标识, 例如：\n",
    "\n",
    "* $(\\mathtt{S}(0))$\n",
    "* $(\\mathtt{S}(1))$\n",
    "* ...\n",
    "* $(\\mathtt{S}(2 \\mathtt{N}))$\n",
    "\n",
    "如果一个statement包含在多个循环中，其instance由所有迭代器的值按循环的顺序来标识, 比如statement `T` 会包含以下这些例子:\n",
    "\n",
    "* $(\\mathtt{T}(0,0))$ for `Z[0] += A[0] * B[0] /* i = 0, j = 0 */`,\n",
    "* $(\\mathtt{T}(0,1))$ for `Z[1] += A[0] * B[1] /* i = 0, j = 1 */`,\n",
    "* ...\n",
    "* $(\\mathtt{T}(\\mathtt{N},\\mathtt{N}))$ for `Z[2*N] += A[N] * B[N] /* i = N, j = N */`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iteration domain\n",
    "\n",
    "一个statement的所有instance的集合被称为 *(iteration) domain*.\n",
    "<!-- The set of all instances of a statement is referred to as *(iteration) domain* of a statement. -->\n",
    "\n",
    "<!-- Iteration domains can be expressed using the set-builder notation. -->\n",
    "迭代域可以使用 set-builder 表示法来表示:\n",
    "\n",
    "比如, $( \\mathcal{D}_\\mathtt{S} = \\{ \\mathtt{S}[i] : 0 \\leq i \\leq \\mathtt{N} \\}. )$\n",
    "表达式 $( 0 \\leq i \\leq \\mathtt{N} )$ 被循环起始(`i=0`)到结束(`i<=N`)所约束.\n",
    "<!-- The loop iterates for $(i)$ ranging from 0 to $(\\mathtt{N})$ inclusive. -->\n",
    "\n",
    "<!-- The set we defined is *named*, i.e. the statement name $( \\mathtt{S} )$ is mentioned before the list of loop iterator names. -->\n",
    "\n",
    "<!-- Note that $( \\mathtt{N} )$ here is treated as a symbolic constant. -->\n",
    "注意, 这里 $( \\mathtt{N} )$ 被看作是符号常量. 在多面体模型中,这些符号常量通常被称为 *(structure) parameters*:\n",
    "<!-- In the context of the polyhedal model, such symbolic constants are referred to as *(structure) parameters*.\n",
    "Parameters are explicitly mentioned in the set description, -->\n",
    "$( \\mathcal{D}\\_\\mathtt{S} = [N] \\rightarrow \\{ \\mathtt{S}[i] : 0 \\leq i \\leq N \\} )$,\n",
    "本质上是将其转换为从参数值到domain set的具体instance的映射。\n",
    "这样的参数集可以在 *isl* 中如下定义:\n",
    "<!-- essentially turning it into a map from the parameter value to a concrete instance of a domain set.\n",
    "Such parametric sets can be defined in *isl* as follows. -->\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { S[i] : 0 <= i <= N }\n"
     ]
    }
   ],
   "source": [
    "import isl\n",
    "D_S = isl.set(\"[N] -> {S[i] : 0 <= i <= N}\")\n",
    "print(D_S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " \n",
    "<!-- Similarly, we can define the iteration domain of the statement `T`, -->\n",
    "同样, 我们可以为statement `T`定义iteration domain,\n",
    "$( \\mathcal{D}_\\mathcal{T} = [N] \\rightarrow \\{ \\mathtt{T}(i,j) : 0 \\leq i,j \\leq N \\} )$.\n",
    "\n",
    "<!-- This domain is defined as a set of two-dimensional vectors.\n",
    "Each component of a vector corresponds to an enclosing loop, in nesting order. -->\n",
    "这个domain被定义为一组二维向量, 向量的每个分量都以嵌套顺序对应于一个封闭循环。\n",
    "\n",
    "### Question\n",
    "\n",
    "<!-- Define the variable `D_T` so that it contains the iteration domain of `T` using *isl* notation, then print it. -->\n",
    "使用 *isl* 表示定义变量 `D_T` 使其包含 `T` 的迭代域，然后打印它。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { T[i, j] : 0 <= i <= N and 0 <= j <= N }\n"
     ]
    }
   ],
   "source": [
    "D_T = isl.set(\"[N] -> { T[i,j] : 0 <= i, j <= N }\")\n",
    "print(D_T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- As you may see, *isl* output uses conjunction (logical *and*) to combine inequalities around different iterators.\n",
    "This can be handy if the loop bounds differ.  Use `and` in *isl* syntax to express conjunction. -->\n",
    "如上所示，*isl* 输出使用连词(logical *and*)来组合不同迭代器周围的不等式,\n",
    "这可以很方便循环边界不同的情况.\n",
    "\n",
    "### Question\n",
    "使用 `and` operator 分离 `i` 和 `j`的bounds来重新定义`D_T` ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { T[i, j] : 0 <= i <= N and 0 <= j <= N }\n"
     ]
    }
   ],
   "source": [
    "D_T = isl.set(\" [N] -> { T[i, j] : 0 <= i <= N and 0 <= j <= N }\")\n",
    "print(D_T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- \n",
    "Note that the output does not necessarily reproduce the textual form of the input.\n",
    "Rather, it represents the same set after simplification.\n",
    "In particular, redundant inequalities are eliminated and components that appear in first positions are expressed using simpler equations.\n",
    "Here is an example of the simplification result. -->\n",
    "\n",
    "注意，print的输出不一定再现输入的文本形式\n",
    "相反，它表示简化后的同一集合, 比如消除了多余的不等式，并使用更简单的方程来表示出现在第一位的分量。\n",
    "比如下面这个例子:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ [i, j] : i >= 0 and j > 0 and j >= -i }\n"
     ]
    }
   ],
   "source": [
    "print(isl.set(\"{[i,j]: i+j >= 0 and i >= 0 and j > 0 and j >= 1}\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling Non-Unit Strides\n",
    "\n",
    "<!-- Consider now the following code snippet that negates every odd element in the array. -->\n",
    "考虑以下代码片段，它取数组中的每个奇数元素的负数。\n",
    "```c\n",
    "const int N;\n",
    "double A[2*N];\n",
    "\n",
    "for (int i = 1; i < 2*N; i += 2)\n",
    "R:  A[i] = -A[i];\n",
    "```\n",
    "\n",
    "如果数组`A`存储的是复数的实部和虚部，那么上面的代码操作计算复共轭.\n",
    "\n",
    "`R` 的迭代域现在应该限制为 `i` 的奇数值, 这可以使用模运算符来实现: \n",
    "<!-- If `A` stores real and imaginary parts of complex numbers, one after other, this operation computes a complex conjugate.\n",
    "\n",
    "The iteration domain of `R` should now be restricted to odd values of `i`.\n",
    "This can be achieved using a modulo operator -->\n",
    "$( \\mathcal{D}_\\mathtt{R}: [N] \\rightarrow \\{ [i] : 0 \\leq i < N \\wedge i\\mod 2 = 1 \\} $)\n",
    "<!-- which can be translated directly into *isl* syntax. -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { R[i] : (1 + i) mod 2 = 0 and 0 <= i < N }\n"
     ]
    }
   ],
   "source": [
    "D_R = isl.set(\"[N] -> {R[i]: (i mod 2) = 1 and 0 <= i < N}\")\n",
    "print(D_R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- As you may see, *isl* transforms the modulo operation into division with floor rounding.\n",
    "This transformation is a combination of two properties of the modulo operation -->\n",
    "*isl* 将模运算转换为带底舍入的除法, 这种转换是模运算的两个属性的组合\n",
    "\n",
    "$( a \\mod b = c \\Leftrightarrow (a + c) \\mod b = 0 )$,\n",
    "\n",
    "$( a \\mod b \\equiv a - b \\lfloor a/b \\rfloor )$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "<!-- Define the set representing the iteration domain of `Q` in the following code, then print it -->\n",
    "在下面的代码中定义代表`Q`的迭代域的集合，然后打印出来\n",
    "```c\n",
    "const int N;\n",
    "double A[2*N];\n",
    "\n",
    "for (int i = 1; i < 2*N; i += 2)\n",
    "Q:  A[i] = -A[i];\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { Q[i] : (1 + i) mod 2 = 0 and 0 <= i < N }\n"
     ]
    }
   ],
   "source": [
    "D_Q = isl.set(\"[N] -> {Q[i]: (i mod 2) = 1 and 0 <= i < N}\")\n",
    "print(D_Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Handling Conditions\n",
    "\n",
    "<!-- Conditional constructs inside the loop also limit the iteration domain of the statements they enclose.\n",
    "The complex conjugate computation can also be re-written using a branch inside the unit-stride loop. -->\n",
    "循环内的条件构造也限制了它们所包含的语句的迭代域,\n",
    "复共轭计算也可以使用单位步长循环内的分支语句来重写。\n",
    "```c\n",
    "const int N;\n",
    "double A[2*N];\n",
    "\n",
    "for (int i = 1; i < 2*N; ++i)\n",
    "    if (i % 2 == 1)\n",
    "P:    A[i] = -A[i];\n",
    "```\n",
    "<!-- The definition of the iteration domain should also include the constraint imposed by the branches surrounding the statement. -->\n",
    "迭代域的定义还应该包括语句周围的分支所施加的约束。\n",
    "\n",
    "### Question\n",
    "<!-- Define the set representing the iteration domain of `P` and print it. -->\n",
    "定义代表`P`的迭代域的集合并打印它。\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { P[i] : (1 + i) mod 2 = 0 and 0 <= i < N }\n"
     ]
    }
   ],
   "source": [
    "D_P = isl.set(\"[N] -> {P[i]: (i mod 2) = 1 and 0 <= i < N}\")\n",
    "print(D_P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- Even though the set of statement instances are identical for `P` and `R`, these domains would be considered different *because of different statement names*. -->\n",
    "即使 `P` 和 `R` 的语句实例集是相同的，这些域也会被认为是不同的因为是不同的*statement name*。\n",
    "\n",
    "### Question\n",
    "如何修改`D_P`, 让他等价于`D_R`?\n",
    "<!-- What do you need to change in the definition of `D_P` to make it equal to `D_R`? \n",
    "Change it and check using `D_P.is_equal(D_R)`. -->\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "D_P = isl.set(\"[N] -> {R[i]: i mod 2 = 1 and 0 <= i < N}\")\n",
    "print(D_P.is_equal(D_R))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iteration Domains as Presburger Sets\n",
    "\n",
    "<!-- Because *isl* operates on Presburger Sets, it can encode any iteration domain that can be expressed using Presburger formulas.\n",
    "This typically involves the statements surrounded by loops and branches with the so called *static control flow*.\n",
    "That is, loop bounds and branch conditions are Presburger formulas of outer loop bounds and parameters.\n",
    "The values of parameters are unknown but must be constant throughout the execution.\n",
    "As a corollary, the control flow cannot depend on *values* being computed.\n",
    "Hence the program parts amenable to the polyhedral modeling are referred to as *static control parts* or *SCoPs*. -->\n",
    "\n",
    "由于*isl*在Presburger Sets上运行，因此它可以编码任何可以使用Presburger公式表示的迭代域,\n",
    "这通常涉及由具有所谓的*static control flow*的循环和分支包围的语句。\n",
    "\n",
    "也就是说，loop bounds和分支条件是外部边界和参数的Presburger公式,\n",
    "其中参数的值未知，但在整个执行过程中必须保持不变。\n",
    "\n",
    "作为推论，控制流不能依赖于被计算的*value*, \n",
    "因此，适合多面体建模的程序部件被称为*static control parts*或*SCoPs*。\n",
    "\n",
    "### Question\n",
    "\n",
    "<!-- Define the iteration domain of the only statement, enclosed in two loops and a branch with a disjunctive constraint.\n",
    "\n",
    "Hint: use the operator `or` and parentheses to ensure the precedence if necessary. -->\n",
    "定义包含在两个循环和一个具有析取约束的分支中的statement的迭代域。\n",
    "\n",
    "Hint: 如果有必要，请使用运算符`or`和括号来确保优先级\n",
    "```c\n",
    "for (int i = 0; i < 10; ++i)\n",
    "  for (int j = 0; j < 10; ++j)\n",
    "    if (i < j - 1 || i > j + 1)\n",
    "      Z[i][j] = 0.;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ Z[i, j] : (i >= 0 and 2 + i <= j <= 9) or (i <= 9 and 0 <= j <= -2 + i) }\n"
     ]
    }
   ],
   "source": [
    "D_Z = isl.set(\" { Z[i,j] : 0 <= i,j < 10 and (i < j - 1 or i > j + 1) } \")\n",
    "print(D_Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### \n",
    "<!-- It is also possible to express certain common mathematical operations using the Presburger formulas -->\n",
    "也可以使用Presburger公式表达某些常见的数学运算:\n",
    "\n",
    "* `i >= max(a,b)` $( \\Leftrightarrow i \\geq a \\wedge i \\geq b )$ (lower bound only)\n",
    "* `i <= min(a,b)` $( \\Leftrightarrow i \\leq a \\wedge i \\leq b )$ (upper bound only)\n",
    "* `a = ceil(b/c)` $( \\Leftrightarrow a = \\lfloor (b - 1)/c \\rfloor + 1 )$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Putting Domains Together\n",
    "<!-- In summary, an iteration domain of a statement is a set of multidimensional vectors constrained by the affine expressions that appear in loop bounds and branch conditions surrounding the statement.\n",
    "\n",
    "Iteration domains of multiple statements live in different *spaces* due to different statement names, even if they are enclosed by the same loops.\n",
    "They can be managed together by putting them into a union set.  For example, the combined domain of -->\n",
    "\n",
    "总之，statement的迭代域是一组受仿射表达式约束的多维向量，仿射表达式出现在statement周围的循环边界和分支条件下。\n",
    "\n",
    "由于statement名称不同，多个语句的迭代域存在于不同的*Spaces*中，即使它们被相同的循环包围,\n",
    "通过将它们放入*Union set*中可以统一操作, 比如组合domain:\n",
    "```c\n",
    "for (int i = 0; i < 10; ++i)\n",
    "  for (int j = 0; j < 10; ++j) {\n",
    "    if (i < j - 1)\n",
    "S1:   Z[i][j] = 0.;\n",
    "    if (i > j + 1)\n",
    "S2:   Z[i][j] = 0.;\n",
    "  }\n",
    "```\n",
    "<!-- is defined as -->\n",
    "被定义为:\n",
    "$( \\mathcal{D} = \n",
    "    \\{\\mathtt{S1}[i,j]: 0 \\leq i,j < 10 \\wedge i < j - 1 \\} \\cup \n",
    "    \\{\\mathtt{S2}[i,j]: 0 \\leq i,j < 10 \\wedge i > j + 1 \\},\n",
    "$)\n",
    "<!-- which translates to *isl* notation as -->\n",
    "*isl*表示如下:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ S2[i, j] : 0 <= i <= 10 and 0 <= j <= 10 and j <= -2 + i; S1[i, j] : 0 <= i <= 10 and j >= 2 + i and 0 <= j <= 10 }\n"
     ]
    }
   ],
   "source": [
    "D = isl.union_set(\"{S1[i,j]: 0 <= i,j <= 10 and i < j - 1; S2[i,j]: 0 <= i,j <= 10 and i > j + 1}\")\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting Iteration Domains\n",
    "可以绘制1D/2D的非参数迭代域集合:\n",
    "<!-- Polyhedral playground supports 1D and 2D plots of non-parametric sets. -->\n",
    "<!-- If the iteration domain has these properties, you can plot it immediately using -->\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALjRJREFUeJzt3X9sXfV9//HXuQ44xva9lt0QJ7NBztXd+NWWkHRTQ/gxtc3U0krVfnRllFEhsSJCfhStgpRKdO0aF1gRLhnJUlVVMkabP2JWJg1nWctCGaqIAmkZ7SBNUHHKoiyO77k2ToOIP98/zjceN7Gvfe89933u5/r5kI6SXs7JfXL1Cefd43PvDZxzTgAAAEZSSQcAAID5heEDAACYYvgAAACmGD4AAIAphg8AAGCK4QMAAJhi+AAAAKYYPgAAgKkFSQeca3JyUm+99Zba29sVBEHSOQAAYA6ccxobG9PSpUuVSpW+tlF3w8dbb72l3t7epDMAAEAFhoeH1dPTU3Kfuhs+2tvbJUXx6XQ64RoAADAXhUJBvb29U+fxUupu+Dj7o5Z0Os3wAQCAZ+ZyywQ3nAIAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAU3X3bpdacU4aGZHGx6W2NqmrS6r3zzCj2QbNNnxslvzsphn1ruGvfOTz0sCAlMtJixZJfX3Rr7lc9Hg+n3Th+Wi2QbMNH5slP7tphjdcnQnD0ElyYRhW/WcNDTnX2upcEERbNFtH29nHWluj/eoFzTZotuFjs3N+dtOMpJVz/m7Y4WNoyLmmJudSqeIFfe6WSkX71cPipplmmpPnYzfNqAc1HT727dvnPvnJT7olS5Y4Se6pp54q+ueTk5PugQcecEuWLHELFy50N9xwg/uv//qvmsTPZHQ0mpZnW9TvXdytrdFxSaGZZpqTbfa1m2bUi3LO32Xf8/H222/rgx/8oLZs2TLtP3/ooYf0yCOPaMuWLdq/f7+6u7v1sY99TGNjY5X/bKhMO3ZIExPS5OTc9p+cjPbfubO2XaXQbINmGz42S3520wwvVTPl6JwrH5OTk667u9t985vfnHrst7/9rctkMm7btm1z+jOrvfIxOelcNnv+zw9n24IgOm5ysqKnrQrNNNOcbLOv3TSjntT0ykcpb7zxho4dO6Y1a9ZMPdbc3KwbbrhBL7zwwrTHnD59WoVCoWirxsiIdPhwtFzL4Vx03MmTVT19RWi2QbMNH5slP7tphq9iHT6OHTsmSVq8eHHR44sXL576Z+fq7+9XJpOZ2np7e6tqGB+v6nAZ/nRoCs02aLbhY7PkZzfN8FVNPufj3K/Tdc7N+BW7mzZtUhiGU9vw8HBVz93WVtXham+v7vhK0GyDZhs+Nkt+dtMMX8U6fHR3d0vSeVc5jh8/ft7VkLOam5uVTqeLtmp0dUnZbPmfjBcE0XGdnVU9fUVotkGzDR+bJT+7aYavYh0++vr61N3drb1790499s4772jfvn1atWpVnE81oyCQ1q2r7Nj165P5OF+abdBsw8dmyc9umuGrwLnybvsZHx/Xr371K0nS8uXL9cgjj+gP//AP1dnZqUsuuUQPPvig+vv79b3vfU+5XE6bN2/Wf/zHf+i1115T+xyulxUKBWUyGYVhWPFVkHxe6umRTp2a21u5UimppUU6elTq6KjoKatGsw2abfjYLPnZTTPqRVnn73LfSvPss886Sedtt912m3Pu/z5krLu72zU3N7vrr7/evfLKKzV5q04p5X563p49VT1dLGimmebk+dhNM+oBH6/+/831ewPqaVHTbINmGz42O+dnN81IGsPHe4yOOjcwEH04zXsXdjYbPZ7Px/I0saLZBs02fGx2zs9umpGkcs7fZd/zUWtx3PMxHeeiD6cZG4veqtXZWf83LtFsg2YbPjZLfnbTjCSUc/6eN8MHAAConXLO3zX5kDEAAICZMHwAAABTDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADC1IOkAK85JIyPS+LjU1iZ1dUlBkHRVaTTboNmGj82Sn90+NmN+afgrH/m8NDAg5XLSokVSX1/0ay4XPZ7PJ114Pppt0GzDx2bJz24fmzFPuToThqGT5MIwrPrPGhpyrrXVuSCItuj/D0Tb2cdaW6P96gXNNmi24WOzc352+9iMxlLO+bthh4+hIeeampxLpYr/Ep67pVLRfvXwF5JmmmlOno/dPjaj8ZRz/g6ccy7Zay/FCoWCMpmMwjBUOp2u6M/I56WeHunUKWlycvb9UymppUU6elTq6KjoKatGsw2abfjYLPnZ7WMzGlM55++GvOdjxw5pYmJufxGlaL+JCWnnztp2lUKzDZpt+Ngs+dntYzPQcFc+nIturjpyJPr9XAWBtGyZdOiQ/V3hNNug2YaPzZKf3T42o3GVc/5uuOHjxIno7u5KnTgRvS3NEs02aLbhY/PZ5/Wt28dmNK55/WOX8fHqjh8bi6ejHDTboNmGj82Sn90+NgNSAw4fbW3VHd/eHk9HOWi2QbMNH5slP7t9bAakBhw+urqkbLb8n2MGQXRcZ2dtukqh2QbNNnxslvzs9rEZkBpw+AgCad26yo5dvz6Zm69otkGzDR+bJT+7fWwGpAa84VTy833vNNug2YaPzZKf3T42ozHN6xtOpegv1O7d0VSfmuXfMJWK9hscTPYvIs02aLbhY7PkZ7ePzUDDfry6c3P/roM9e2IIjwnNNmi24WOzc352+9iMxsJ3u7zH6KhzAwPOZbPFfxmz2ejxfD6Wp4kVzTZotuFjs3N+dvvYjMYx77/bZTrOSSdPRu9rb2+P7vKu95utaLZBsw0fmyU/u31shv/m9SecAgAAe/P+hlMAAFC/GD4AAIAphg8AAGCK4QMAAJhi+AAAAKYYPgAAgCmGDwAAYIrhAwAAmGL4AAAAphg+AACAKYYPAABgiuEDAACYYvgAAACmGD4AAIAphg8AAGCK4QMAAJhakHSAFeekkRFpfFxqa5O6uqQgSLqqNJpt0GzDx2bJ325gOvWynhv+ykc+Lw0MSLmctGiR1NcX/ZrLRY/n80kXno9mGzTb8LFZ8rcbmE7drWdXZ8IwdJJcGIZV/1lDQ861tjoXBNEWzXzRdvax1tZov3pBsw2abfjY7Jy/3cB0rNZzOefvhh0+hoaca2pyLpUqfqHP3VKpaL96+I8IzTTTnDxfu4HpWK7ncs7fgXPOxXkl5d1339VXv/pV/dM//ZOOHTumJUuW6POf/7y+8pWvKJWa/ac8hUJBmUxGYRgqnU5X1JDPSz090qlT0uTk7PunUlJLi3T0qNTRUdFTVo1mGzTb8LFZ8rcbmI71ei7n/B37PR8PPvigtm3bpi1btuiXv/ylHnroIT388MN67LHH4n6qGe3YIU1MzO3FlqL9JiaknTtr21UKzTZotuFjs+RvNzCdel7PsV/5+OQnP6nFixfru9/97tRjf/Inf6KLLrpI//iP/zjr8dVe+XAuuoHmyJHo93MVBNKyZdKhQ/Z3/tJsg2YbPjZL/nYD00liPSd65WP16tX60Y9+pNdff12S9LOf/UzPP/+8PvGJT0y7/+nTp1UoFIq2aoyMSIcPl/diS9H+hw9LJ09W9fQVodkGzTZ8bJb87QamU+/rOfbh495779XNN9+syy67TBdccIGWL1+ujRs36uabb552//7+fmUymamtt7e3qucfH6/qcI2NVXd8JWi2QbMNH5slf7uB6dT7eo59+Ni1a5eeeOIJPfnkk3rppZe0Y8cO/d3f/Z127Ngx7f6bNm1SGIZT2/DwcFXP39ZW1eFqb6/u+ErQbINmGz42S/52A9Op9/Uc+/DxpS99Sffdd58++9nP6v3vf79uvfVWffGLX1R/f/+0+zc3NyudThdt1ejqkrLZ8n9WFQTRcZ2dVT19RWi2QbMNH5slf7uB6dT7eo59+JiYmDjvLbVNTU2anOvttlUKAmndusqOXb8+mRvGaLZBsw0fmyV/u4Hp1Pt6jv3dLp///Of17//+7/qHf/gHXXnllXr55Zf1V3/1V7r99tv14IMPzno8n/NBcy3RbMPHZsnfbmA68+pzPh577DH96Z/+qe666y5dfvnl+uu//mt94Qtf0Ne//vW4n2pGHR3S7t3R5Dbb55qlUtF+g4PJ/seDZhs02/CxWfK3G5hOXa/nyj9ItTaS+G6XPXtiCI8JzTZotuFjs3P+dgPTsVrPfLfLe4yOOjcw4Fw2W/yCZ7PR4/l8LE8TK5pt0GzDx2bn/O0GpmOxnhP9bpdqxXHPx3Sciz40ZWwsegtRZ2f93yBGsw2abfjYLPnbDUynluu5nPP3vBk+AABA7SR6wykAAEApDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADC1IOkAK85JIyPS+LjU1iZ1dUlBkHRVaTTboNmGj81AKT6u6XppbvgrH/m8NDAg5XLSokVSX1/0ay4XPZ7PJ114Pppt0GzDx2agFB/XdN01uzoThqGT5MIwrPrPGhpyrrXVuSCItmjmi7azj7W2RvvVC5pt0GzDx2agFB/XtFVzOefvhh0+hoaca2pyLpUqfqHP3VKpaL96WCg000wzUL98XNOWzeWcvwPnnDO+2FJSoVBQJpNRGIZKp9MV/Rn5vNTTI506JU1Ozr5/KiW1tEhHj0odHRU9ZdVotkGzDR+bgVJ8XNPWzeWcvxvyno8dO6SJibm92FK038SEtHNnbbtKodkGzTZ8bAZK8XFN13Nzw135cC66gebIkej3cxUE0rJl0qFD9nf+0myDZhs+NgOl+Limk2gu5/zdcMPHiRPRHbyVOnEieuuRJZpt0GzDx2agFB/XdBLN8/rHLuPj1R0/NhZPRzlotkGzDR+bgVJ8XNP13txww0dbW3XHt7fH01EOmm3QbMPHZqAUH9d0vTc33PDR1SVls+X/rCoIouM6O2vTVQrNNmi24WMzUIqPa7remxtu+AgCad26yo5dvz6ZG91otkGzDR+bgVJ8XNP13txwN5xKvB/bCs02aAaS5+Oa5nM+jHV0SLt3R5NbapZ/w1Qq2m9wMNn/6NFsg2YbPjYDpfi4puu6ufIPUq2NJL7bZc+eGMJjQrMNmm342AyU4uOatmrmu13eY3TUuYEB57LZ4hc8m40ez+djeZpY0WyDZhs+NgOl+LimLZrn/Xe7TMc56eTJ6L3L7e3Rnbz1fmMbzTZotuFjM1CKj2u6ls3z+hNOAQCAvXl/wykAAKhfDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADDF8AEAAEwtSDrAinPSyIg0Pi61tUldXVIQJF1VGs02aAaS5+OaprlyDX/lI5+XBgakXE5atEjq64t+zeWix/P5pAvPR7MNmoHk+bimaY6BqzNhGDpJLgzDqv+soSHnWludC4Joi2a+aDv7WGtrtF+9oNkGzUDyfFzTNM+snPN3ww4fQ0PONTU5l0oVv9DnbqlUtF89LBSaaW6kZqAUH9c0zaWVc/4OnHMu7qspv/nNb3TvvffqmWee0alTp/S7v/u7+u53v6sVK1bMemyhUFAmk1EYhkqn0xU9fz4v9fRIp05Jk5Oz759KSS0t0tGjUkdHRU9ZNZpt0Awkz8c1TfPsyjl/x37Px+joqK699lpdcMEFeuaZZ/SLX/xC3/rWt9Rh+Orv2CFNTMztxZai/SYmpJ07a9tVCs02aAaS5+OapjlesV/5uO+++/Sf//mf+slPflLR8dVe+XAuuoHmyJHo93MVBNKyZdKhQ/Z3/tJsg2YgeT6uaZrnppzzd+zDxxVXXKE/+qM/0tGjR7Vv3z79zu/8ju666y7dcccd0+5/+vRpnT59uii+t7e34uHjxInoDt5KnTgRvfXIEs02aAaS5+OapnluEv2xy5EjR7R161blcjnt2bNHd955p9avX6+dM1zH6e/vVyaTmdp6e3urev7x8aoO19hYdcdXgmYbNAPJ83FN0xy/2IePyclJXXPNNdq8ebOWL1+uL3zhC7rjjju0devWaffftGmTwjCc2oaHh6t6/ra2qg5Xe3t1x1eCZhs0A8nzcU3THL/Yh48lS5boiiuuKHrs8ssv15tvvjnt/s3NzUqn00VbNbq6pGy2/J9VBUF0XGdnVU9fEZpt0Awkz8c1TXP8Yh8+rr32Wr322mtFj73++uu69NJL436qaQWBtG5dZceuX5/MzXk026AZSJ6Pa5rm+MV+w+n+/fu1atUq/c3f/I0+85nP6MUXX9Qdd9yh7du365Zbbpn1eD7ng+ZaohlIno9rmubZJXrD6Yc+9CE99dRT+v73v6+rrrpKX//61/Xoo4/OafCIS0eHtHt3NLmlZvk3TKWi/QYHk/0PNc02aAaS5+OapjlmlX+Qam0k8d0ue/bEEB4Tmm3QDCTPxzVN88z4bpf3GB11bmDAuWy2+AXPZqPH8/lYniZWNNugGUiej2ua5ukl/t0u1Yjjno/pOCedPBm9d7m9PbqTt95vxqPZBs1A8nxc0zQXS/QTTqtVq+EDAADUTqI3nAIAAJTC8AEAAEwxfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMLkg6w4pw0MiKNj0ttbVJXlxQESVeVRrMNH5uBUnxc0zTbqJfmhr/ykc9LAwNSLictWiT19UW/5nLR4/l80oXno9mGj81AKT6uaZpt1F2zqzNhGDpJLgzDqv+soSHnWludC4Joi2a+aDv7WGtrtF+9oNmGj81AKT6uaZptWDWXc/5u2OFjaMi5pibnUqniF/rcLZWK9quHhUIzzUAlfFzTNDdecznn78A554wvtpRUKBSUyWQUhqHS6XRFf0Y+L/X0SKdOSZOTs++fSkktLdLRo1JHR0VPWTWabfjYDJTi45qm2YZ1cznn74a852PHDmliYm4vthTtNzEh7dxZ265SaLbhYzNQio9rmmYb9dzccFc+nItuoDlyJPr9XAWBtGyZdOiQ/Z2/NNvwsRkoxcc1TbONJJrLOX833PBx4kR0B2+lTpyI3npkiWYbPjYDpfi4pmm2kUTzvP6xy/h4dcePjcXTUQ6abfjYDJTi45qm2Ua9Nzfc8NHWVt3x7e3xdJSDZhs+NgOl+LimabZR780NN3x0dUnZbPk/qwqC6LjOztp0lUKzDR+bgVJ8XNM026j35oYbPoJAWreusmPXr0/mhkKabfjYDJTi45qm2Ua9NzfcDacS78e2QjOQPB/XNM02+JwPYx0d0u7d0eSWmuXfMJWK9hscTPbkQrMNH5uBUnxc0zTbqOvmyj9ItTaS+G6XPXtiCI8JzTZ8bAZK8XFN02zDqpnvdnmP0VHnBgacy2aLX/BsNno8n4/laWJFsw0fm4FSfFzTNNuwaJ733+0yHeekkyej9y63t0d38tb7DYQ02/CxGSjFxzVNs41aNs/rTzgFAAD25v0NpwAAoH4xfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMMHwAAwBTDBwAAMLUg6QArzkkjI9L4uNTWJnV1SUGQdFVpNAPJ83FN02yD5so1/JWPfF4aGJByOWnRIqmvL/o1l4sez+eTLjwfzUDyfFzTNNugOQauzoRh6CS5MAyr/rOGhpxrbXUuCKItmvmi7exjra3RfvWCZiB5Pq5pmm3QPLNyzt8NO3wMDTnX1ORcKlX8Qp+7pVLRfvWwUGgGkufjmqaZ5npoLuf8HTjnXC2vrPT39+vLX/6yNmzYoEcffXTW/QuFgjKZjMIwVDqdrug583mpp0c6dUqanJx9/1RKammRjh6VOjoqesqq0Qwkz8c1TbMNmmdXzvm7pvd87N+/X9u3b9cHPvCBWj7NeXbskCYm5vZiS9F+ExPSzp217SqFZiB5Pq5pmm3QHK+aXfkYHx/XNddco8cff1x/+7d/q6uvvtrkyodz0Q00R45Ev5+rIJCWLZMOHbK/85dmIHk+rmmabdA8N3Vx5WPt2rW66aab9NGPfrTkfqdPn1ahUCjaqjEyIh0+XN6LLUX7Hz4snTxZ1dNXhGYgeT6uaZpt0By/mgwfP/jBD/TSSy+pv79/1n37+/uVyWSmtt7e3qqee3y8qsM1Nlbd8ZWgGUiej2uaZhs0xy/24WN4eFgbNmzQE088oYULF866/6ZNmxSG4dQ2PDxc1fO3tVV1uNrbqzu+EjQDyfNxTdNsg+b4xT58HDhwQMePH9eKFSu0YMECLViwQPv27dO3v/1tLViwQGfOnCnav7m5Wel0umirRleXlM2W/7OqIIiO6+ys6ukrQjOQPB/XNM02aI5f7MPHRz7yEb3yyis6ePDg1LZy5UrdcsstOnjwoJqamuJ+yiJBIK1bV9mx69cncxMkzUDyfFzTNNugOX41/5wPSbrxxhvN3u0i8X5sKz42A6X4uKZptkHz7Ori3S5J6uiQdu+OJrfULP+GqVS03+BgsidEmoHk+bimabZBc8wq/yDV2kjiu1327IkhPCY0A8nzcU3TbIPmmfHdLu8xOurcwIBz2WzxC57NRo/n87E8TaxoBpLn45qm2QbN06ur73YpVxz3fEzHuehDU8bGorcQdXbW/02PNAPJ83FN02yD5mLlnL/nzfABAABqZ97fcAoAAOoXwwcAADDF8AEAAEwxfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMLkg6w4pw0MiKNj0ttbVJXlxQESVeV5mMzMBNf17OP3TTboLlyDX/lI5+XBgakXE5atEjq64t+zeWix/P5pAvP52MzMBNf17OP3TTboDkGrs6EYegkuTAMq/6zhoaca211LgiiLZr5ou3sY62t0X71wsdmYCa+rmcfu2m2QfPMyjl/N+zwMTTkXFOTc6lU8Qt97pZKRfvVw0LxsRmYia/r2cdummmuh+Zyzt+Bc84ZX2wpqVAoKJPJKAxDpdPpiv6MfF7q6ZFOnZImJ2ffP5WSWlqko0eljo6KnrJqPjYDM/F1PfvYTbMNmmdXzvm7Ie/52LFDmpiY24stRftNTEg7d9a2qxQfm4GZ+Lqefeym2QbN8Wq4Kx/ORTfQHDkS/X6ugkBatkw6dMj+zl8fm4GZ+Lqefeym2QbNc1PO+bvhho8TJ6I7eCt14kT01iNLPjYDM/F1PfvYTbMNmudmXv/YZXy8uuPHxuLpKIePzcBMfF3PPnbTbIPm+DXc8NHWVt3x7e3xdJTDx2ZgJr6uZx+7abZBc/wabvjo6pKy2fJ/VhUE0XGdnbXpKsXHZmAmvq5nH7tptkFz/Bpu+AgCad26yo5dvz6ZG918bAZm4ut69rGbZhs0x6/hbjiVeD82kDRf17OP3TTboHl28/qGUyl60Xbvjia31Cz/hqlUtN/gYLL/0fOxGZiJr+vZx26abdAcs8o/SLU2kvhulz17YgiPiY/NwEx8Xc8+dtNsg+aZ8d0u7zE66tzAgHPZbPELns1Gj+fzsTxNrHxsBmbi63r2sZtmGzRPb95/t8t0nJNOnozeu9zeHt3JW+83avrYDMzE1/XsYzfNNmguNq8/4RQAANib9zecAgCA+sXwAQAATDF8AAAAUwwfAADAFMMHAAAwxfABAABMMXwAAABTDB8AAMAUwwcAADDF8AEAAEwxfAAAAFMMHwAAwBTDBwAAMMXwAQAATDF8AAAAUwuSDrDinDQyIo2PS21tUleXFARJVwGV8XE9+9gs+dlNsw2aK9fwVz7yeWlgQMrlpEWLpL6+6NdcLno8n0+6EJg7H9ezj82Sn90026A5Bq7OhGHoJLkwDKv+s4aGnGttdS4Ioi2a+aLt7GOtrdF+QL3zcT372Oycn90026B5ZuWcvxt2+Bgacq6pyblUqviFPndLpaL96mmhAOfycT372Oycn90001wPzeWcvwPnnIvzSkp/f78GBwf13//932ppadGqVav04IMP6vd+7/fmdHyhUFAmk1EYhkqn0xU15PNST4906pQ0OTn7/qmU1NIiHT0qdXRU9JRAzfi4nn1slvzsptkGzbMr5/wd+z0f+/bt09q1a/XTn/5Ue/fu1bvvvqs1a9bo7bffjvupZrRjhzQxMbcXW4r2m5iQdu6sbRdQCR/Xs4/Nkp/dNNugOV6xX/k41//+7//q4osv1r59+3T99dfPun+1Vz6ci26gOXIk+v1cBYG0bJl06FD9362M+cPH9exjs+RnN802aJ6bRK98nCsMQ0lSZ2fntP/89OnTKhQKRVs1Rkakw4fLe7GlaP/Dh6WTJ6t6eiBWPq5nH5slP7tptkFz/Go6fDjndM8992j16tW66qqrpt2nv79fmUxmauvt7a3qOcfHqzpcY2PVHQ/Eycf17GOz5Gc3zTZojl9Nh4+7775bP//5z/X9739/xn02bdqkMAyntuHh4aqes62tqsPV3l7d8UCcfFzPPjZLfnbTbIPm+NVs+Fi3bp2efvppPfvss+rp6Zlxv+bmZqXT6aKtGl1dUjZb/s+qgiA6boafDgGJ8HE9+9gs+dlNsw2a4xf78OGc0913363BwUH9+Mc/Vl9fX9xPUVIQSOvWVXbs+vXcbIr64uN69rFZ8rObZhs0xy/2d7vcddddevLJJ/XDH/6w6LM9MpmMWlpaZj2ez/kAivm4nn1slvzsptkGzbNL9N0uW7duVRiGuvHGG7VkyZKpbdeuXXE/1Yw6OqTdu6PJLTXLv2EqFe03OMjggfrk43r2sVnys5tmGzTHrPIPUq2NJL7bZc+eGMKBGvNxPfvY7Jyf3TTboHlmfLfLe4yOOjcw4Fw2W/yCZ7PR4/l8LE8DmPBxPfvY7Jyf3TTboHl6iX63S7XiuOdjOs5FH5oyNha9haizk5tL4S8f17OPzZKf3TTboLlYOefveTN8AACA2qmrj1cHAAB4L4YPAABgiuEDAACYYvgAAACmGD4AAIAphg8AAGCK4QMAAJhi+AAAAKYYPgAAgCmGDwAAYIrhAwAAmGL4AAAAphg+AACAKYYPAABgiuEDAACYYvgAAACmFiQdACTNOWlkRBofl9rapK4uKQiSriqNZjs+dtNsg+bKceUD81Y+Lw0MSLmctGiR1NcX/ZrLRY/n80kXno9mOz5202yD5hi4OhOGoZPkwjBMOgUNbGjIudZW54Ig2qL/PxBtZx9rbY32qxc02/Gxm2YbNM+snPM3wwfmnaEh55qanEuliv8SnrulUtF+9fAfEZrt+NhNM8310FzO+Ttwzjnjiy0lFQoFZTIZhWGodDqddA4aTD4v9fRIp05Jk5Oz759KSS0t0tGjUkdHreumR7MdH7tptkHz7Mo5f3PPB+aVHTukiYm5/UWUov0mJqSdO2vbVQrNdnzsptkGzfHiygfmDeeim6uOHIl+P1dBIC1bJh06ZH9XOM12fOym2QbNc1PO+ZvhA/PGiRPR3d3VHN/VFV/PXJ+TZhs+dtNsg+a54ccuwDTGx6s7fmwsno5y0GzHx26abdAcP4YPzBttbdUd394eT0c5aLbjYzfNNmiOH8MH5o2uLimbLf/nmEEQHdfZWZuuUmi242M3zTZojh/DB+aNIJDWravs2PXrk7kJkmY7PnbTbIPm+HHDKeYV3qtvw8dmyc9umm3QPDtuOAVm0NEh7d4dTfWpWVZ/KhXtNziY7AmRZjs+dtNsg+aYVf5BqrXBx6vDwly/62DPnqRL/w/NdnzsptkGzTPju12AORgddW5gwLlstvgvYzYbPZ7PJ114Pprt+NhNsw2ap8d3uwBlcE46eTJ6X3t7e3SXd1I3Pc4VzXZ87KbZBs3F+IRTAABgihtOAQBA3WL4AAAAphg+AACAKYYPAABgiuEDAACYYvgAAACmGD4AAIAphg8AAGCK4QMAAJhi+AAAAKYYPgAAgCmGDwAAYIrhAwAAmGL4AAAAphg+AACAqQVJB6CxOCeNjEjj41Jbm9TVJQVB0lWl0WzDx2bJz26abdBcOa58IBb5vDQwIOVy0qJFUl9f9GsuFz2ezyddeD6abfjYLPnZTbMNmmPg6kwYhk6SC8Mw6RTM0dCQc62tzgVBtEWzdbSdfay1NdqvXtBsw8dm5/zsptkGzTMr5/zN8IGqDA0519TkXCpVvKDP3VKpaL96+AtJM82l+NhNM8310FzO+TtwzrlaXFF5/PHH9fDDD+t//ud/dOWVV+rRRx/VddddN+txhUJBmUxGYRgqnU7XIg0xyeelnh7p1ClpcnL2/VMpqaVFOnpU6uiodd30aLbhY7PkZzfNNmieXTnn75rc87Fr1y5t3LhR999/v15++WVdd911+vjHP64333yzFk+HhOzYIU1MzG1RS9F+ExPSzp217SqFZhs+Nkt+dtNsg+Z41eTKxx/8wR/ommuu0datW6ceu/zyy/XpT39a/f39JY/lyocfnItuVDpyJPr9XAWBtGyZdOiQ/R3WNNvwsVnys5tmGzTPTaJXPt555x0dOHBAa9asKXp8zZo1euGFF87b//Tp0yoUCkUb6t/IiHT4cHmLWor2P3xYOnmyNl2l0GzDx2bJz26abdAcv9iHjxMnTujMmTNavHhx0eOLFy/WsWPHztu/v79fmUxmauvt7Y07CTUwPl7d8WNj8XSUg2YbPjZLfnbTbIPm+NXscz6Cc67XOOfOe0ySNm3apDAMp7bh4eFaJSFGbW3VHd/eHk9HOWi24WOz5Gc3zTZojl/sw8f73vc+NTU1nXeV4/jx4+ddDZGk5uZmpdPpog31r6tLymbL/5lgEETHdXbWpqsUmm342Cz52U2zDZrjF/vwceGFF2rFihXau3dv0eN79+7VqlWr4n46JCQIpHXrKjt2/fpkbiik2YaPzZKf3TTboDl+NXm3y65du3Trrbdq27Zt+vCHP6zt27frO9/5jl599VVdeumlJY/l3S7+4H3vNmi242M3zTZonl3in/Px53/+53r00Uf1ta99TVdffbWee+45/eu//uusgwf80tEh7d4dTcipWVZSKhXtNziY7MmFZhs+Nkt+dtNsg+aYVf5BqrXBx6v7Z67fG7BnT9Kl/4dmGz42O+dnN802aJ4Z3+0Cc6Ojzg0MOJfNFi/sbDZ6PJ9PuvB8NNvwsdk5P7tptkHz9Oriu10qxT0ffnMu+nCasbHorVqdncndQDhXNNvwsVnys5tmGzQXK+f8zfABAACqlvgNpwAAADNh+AAAAKYYPgAAgCmGDwAAYIrhAwAAmGL4AAAAphYkHXCus+/8LRQKCZcAAIC5OnvenssneNTd8DE2NiZJ6u3tTbgEAACUa2xsTJlMpuQ+dfchY5OTk3rrrbfU3t6uIOaPiisUCurt7dXw8DAfYFZDvM42eJ3t8Frb4HW2UavX2TmnsbExLV26VKlZvsmu7q58pFIp9fT01PQ50uk0C9sAr7MNXmc7vNY2eJ1t1OJ1nu2Kx1nccAoAAEwxfAAAAFPzavhobm7WAw88oObm5qRTGhqvsw1eZzu81jZ4nW3Uw+tcdzecAgCAxjavrnwAAIDkMXwAAABTDB8AAMAUwwcAADA1b4aPxx9/XH19fVq4cKFWrFihn/zkJ0knNZz+/n596EMfUnt7uy6++GJ9+tOf1muvvZZ0VsPr7+9XEATauHFj0ikN5ze/+Y0+97nPqaurSxdddJGuvvpqHThwIOmshvLuu+/qK1/5ivr6+tTS0qJly5bpa1/7miYnJ5NO895zzz2nT33qU1q6dKmCINA///M/F/1z55y++tWvaunSpWppadGNN96oV1991aRtXgwfu3bt0saNG3X//ffr5Zdf1nXXXaePf/zjevPNN5NOayj79u3T2rVr9dOf/lR79+7Vu+++qzVr1ujtt99OOq1h7d+/X9u3b9cHPvCBpFMazujoqK699lpdcMEFeuaZZ/SLX/xC3/rWt9TR0ZF0WkN58MEHtW3bNm3ZskW//OUv9dBDD+nhhx/WY489lnSa995++2198IMf1JYtW6b95w899JAeeeQRbdmyRfv371d3d7c+9rGPTX3HWk25eeD3f//33Z133ln02GWXXebuu+++hIrmh+PHjztJbt++fUmnNKSxsTGXy+Xc3r173Q033OA2bNiQdFJDuffee93q1auTzmh4N910k7v99tuLHvvjP/5j97nPfS6hosYkyT311FNT/3tyctJ1d3e7b37zm1OP/fa3v3WZTMZt27at5j0Nf+XjnXfe0YEDB7RmzZqix9esWaMXXnghoar5IQxDSVJnZ2fCJY1p7dq1uummm/TRj3406ZSG9PTTT2vlypX6sz/7M1188cVavny5vvOd7ySd1XBWr16tH/3oR3r99dclST/72c/0/PPP6xOf+ETCZY3tjTfe0LFjx4rOjc3NzbrhhhtMzo1198VycTtx4oTOnDmjxYsXFz2+ePFiHTt2LKGqxuec0z333KPVq1frqquuSjqn4fzgBz/QSy+9pP379yed0rCOHDmirVu36p577tGXv/xlvfjii1q/fr2am5v1l3/5l0nnNYx7771XYRjqsssuU1NTk86cOaNvfOMbuvnmm5NOa2hnz3/TnRt//etf1/z5G374OCsIgqL/7Zw77zHE5+6779bPf/5zPf/880mnNJzh4WFt2LBB//Zv/6aFCxcmndOwJicntXLlSm3evFmStHz5cr366qvaunUrw0eMdu3apSeeeEJPPvmkrrzySh08eFAbN27U0qVLddtttyWd1/CSOjc2/PDxvve9T01NTedd5Th+/Ph5Ex/isW7dOj399NN67rnn1NPTk3ROwzlw4ICOHz+uFStWTD125swZPffcc9qyZYtOnz6tpqamBAsbw5IlS3TFFVcUPXb55Zdr9+7dCRU1pi996Uu677779NnPflaS9P73v1+//vWv1d/fz/BRQ93d3ZKiKyBLliyZetzq3Njw93xceOGFWrFihfbu3Vv0+N69e7Vq1aqEqhqTc0533323BgcH9eMf/1h9fX1JJzWkj3zkI3rllVd08ODBqW3lypW65ZZbdPDgQQaPmFx77bXnvVX89ddf16WXXppQUWOamJhQKlV8KmpqauKttjXW19en7u7uonPjO++8o3379pmcGxv+yock3XPPPbr11lu1cuVKffjDH9b27dv15ptv6s4770w6raGsXbtWTz75pH74wx+qvb196mpTJpNRS0tLwnWNo729/bz7aFpbW9XV1cX9NTH64he/qFWrVmnz5s36zGc+oxdffFHbt2/X9u3bk05rKJ/61Kf0jW98Q5dccomuvPJKvfzyy3rkkUd0++23J53mvfHxcf3qV7+a+t9vvPGGDh48qM7OTl1yySXauHGjNm/erFwup1wup82bN+uiiy7SX/zFX9Q+rubvp6kTf//3f+8uvfRSd+GFF7prrrmGt3/WgKRpt+9973tJpzU83mpbG//yL//irrrqKtfc3Owuu+wyt3379qSTGk6hUHAbNmxwl1xyiVu4cKFbtmyZu//++93p06eTTvPes88+O+1/k2+77TbnXPR22wceeMB1d3e75uZmd/3117tXXnnFpC1wzrnajzgAAACRhr/nAwAA1BeGDwAAYIrhAwAAmGL4AAAAphg+AACAKYYPAABgiuEDAACYYvgAAACmGD4AAIAphg8AAGCK4QMAAJhi+AAAAKb+H/32jjLQKXAKAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "from plot.plotter import plot_set_points\n",
    "\n",
    "\n",
    "plot_set_points(D,color='blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- visualizations can be useful to see the size of the domain or whether individual domains intersect. -->\n",
    "可视化可以有效的检查domain的size或者独立的domain交集.\n",
    "<!-- In this example, domains are completely disjoint, which means they can be traversed by separate loop nests. -->\n",
    "这里例子中,domain完全不相交,意味他们可以被分离的循环所穿过.\n",
    "\n",
    "<!-- If the domain is parametric, you will have to first *fix* all parameter values to a constant by -->\n",
    "如果domain是参数化的,我们首先需要*fix*所有的参数到constant, 通过以下几种方法:\n",
    "\n",
    "* 创建一个参数集，其中域值是固定的;\n",
    "* 将domain与这个新set相交;\n",
    "* 映射所有参数.\n",
    "\n",
    "<!-- Here is how it can be done for `T` statement of our first example. -->\n",
    "下面是对第一个例子的`T` statement 的操作方法:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIvNJREFUeJzt3W9s1eX9//HX5xxCbdrT05xamaTFlOMZmwiKsDmIOv+SEGMk2w01hrDEO25QMMQbQxPnDEldSBbOYmTDGxvN4uAGRU3mDiNxgy2GBdBGwtxSyxbbCKIt5/SUsho51+8GaX/fKn/O+fQ673N6+nwkn+S7fnt6rjxtcl45PZwTOOecAAAAPIhU+gAAAKB2MCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeDPH+g4LhYI++eQTxWIxBUFgffcAACAE55zy+bzmz5+vSOTKz0uYD4tPPvlE7e3t1ncLAAA8GBgYUFtb2xX//+bDIhaLSbp0sKamJuu7BwAAIYyMjKi9vX3ycfxKzIfFxJ8/mpqaGBYAAMww13oZAy/eBAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADemP+rkHJwzmloaEijo6NqbGxUS0sLb75VBnS2Q2sbdLZBZxvV0nlGP2ORzWaVTqeVSqXU2tqqjo4Otba2KpVKKZ1OK5vNVvqINYHOdmhtg8426Gyj6jo7Y7lczklyuVxuWj8nk8m4hoYGFwSBC4LASZq8Jr7W0NDgMpmMp5PPTnS2Q2sbdLZBZxuWnYt9/J6RwyKTybhoNOoikciUiF+9IpGIi0aj/OKGRGc7tLZBZxt0tmHdudjH78A554p9duPFF1/Uz3/+8ylfmzdvns6cOVPsj9DIyIji8bhyuVyod97MZrNqa2vThQsXVCgUrvn9kUhE9fX1GhwcVHNzc8n3N1vR2Q6tbdDZBp1tVKJzsY/fJb/GYvHixTp9+vTkdeLEiVAHDGv37t0aGxsrKqR06dNUx8bG1N3dXeaT1RY626G1DTrboLONau5c8jMWb7zxhnp7e0Pf4XSesXDOKZVK6dSpUyrh2AqCQAsXLlRfXx+vRC4Cne3Q2gadbdDZRqU6l+0Zi76+Ps2fP18dHR16/PHHderUqat+//j4uEZGRqZcYQ0NDam/v7+kkNKl/wj9/f0aHh4Ofd+zCZ3t0NoGnW3Q2Ua1dy5pWNx5553q7u7WgQMH9Nprr+nMmTNatWqVhoaGrnibrq4uxePxyau9vT30YUdHR0PfVpLy+fy0bj9b0NkOrW3Q2QadbVR755KGxZo1a/TDH/5QS5Ys0YMPPqg//vGPki79redKtm7dqlwuN3kNDAyEPmxjY2Po20q65mfI4xI626G1DTrboLONau88rTfIamho0JIlS9TX13fF76mrq1NTU9OUK6yWlhYlk8mS/zYUBIGSyaQSiUTo+55N6GyH1jbobIPONqq987SGxfj4uD788EPdeOONvs5zVUEQqLOzM9RtN23axIuCikRnO7S2QWcbdLZR7Z1L+lchzz77rB555BEtWLBAZ8+e1bZt23To0CGdOHFCN910U1E/g/exmBnobIfWNuhsg842auZ9LAYHB/XEE09o0aJF+sEPfqC5c+fqyJEjRY8KH5qbm7Vv3z4FQaBI5OrHj0QiCoJAPT09/MKWiM52aG2DzjbobKOqO0/r/T1DsP6skAMHDng6+exEZzu0tkFnG3S2Ydm5pj8rZMK5c+dcOp12yWRySsxkMunS6bTLZrMeTgw626G1DTrboLMNq85l+awQH6b7GovLcc5peHhY+XxesVhMiUSCFwGVAZ3t0NoGnW3Q2Ua5Oxf7+F0TwwIAAJRX2d7SGwAA4EoYFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwJs5lT6AD845DQ0NaXR0VI2NjWppaVEQBJU+Vs2hsx1a26CzDTrbqJbOM/oZi2w2q3Q6rVQqpdbWVnV0dKi1tVWpVErpdFrZbLbSR6wJdLZDaxt0tkFnG1XX2RnL5XJOksvlctP6OZlMxjU0NLggCFwQBE7S5DXxtYaGBpfJZDydfHaisx1a26CzDTrbsOxc7OP3jBwWmUzGRaNRF4lEpkT86hWJRFw0GuUXNyQ626G1DTrboLMN687FPn4HzjkX9tmOrq4uPffcc9q8ebN27NhR1G1GRkYUj8eVy+XU1NRU8n1ms1m1tbXpwoULKhQK1/z+SCSi+vp6DQ4Oqrm5ueT7m63obIfWNuhsg842KtG52Mfv0K+xOHr0qHbt2qWlS5eG/RGh7N69W2NjY0WFlKRCoaCxsTF1d3eX+WS1hc52aG2DzjbobKOqO4d5OiSfz7tUKuUOHjzovv/977vNmzd7fyrlcgqFgksmk1/7O9K1riAIXDKZdIVCoeT7nI3obIfWNuhsg842KtW52MfvUM9YbNiwQQ8//LAefPDBa37v+Pi4RkZGplxhDQ0Nqb+/X67Ev94459Tf36/h4eHQ9z2b0NkOrW3Q2QadbVR755KHxZ49e/Tee++pq6urqO/v6upSPB6fvNrb20s+5ITR0dHQt5WkfD4/rdvPFnS2Q2sbdLZBZxvV3rmkYTEwMKDNmzfr97//va677rqibrN161blcrnJa2BgINRBJamxsTH0bSUpFotN6/azBZ3t0NoGnW3Q2Ua1dy5pWBw/flxnz57V8uXLNWfOHM2ZM0eHDh3Sr371K82ZM0cXL1782m3q6urU1NQ05QqrpaVFyWSy5HcSC4JAyWRSiUQi9H3PJnS2Q2sbdLZBZxvV3rmkYfHAAw/oxIkT6u3tnbxWrFihJ598Ur29vYpGo+U6p6RLUTo7O0PddtOmTbyFbJHobIfWNuhsg842qr3ztN7HQpLuvfde3X777byPRY2hsx1a26CzDTrbqMn3saiU5uZm7du3T0EQKBK5+vEjkYiCIFBPTw+/sCWisx1a26CzDTrbqOrOof4x6zRYf1bIgQMHPJ18dqKzHVrboLMNOtuw7FzTnxUy4dy5cy6dTrtkMjklZjKZdOl02mWzWQ8nBp3t0NoGnW3Q2YZVZ5PPCgljuq+xuBznnIaHh5XP5xWLxZRIJHgRUBnQ2Q6tbdDZBp1tlLtzsY/fNTEsAABAedXsizcBAED1YlgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8mVPpA/jgnNPQ0JBGR0fV2NiolpYWBUFQ6WPVHDrbobUNOtugs41q6Tyjn7HIZrNKp9NKpVJqbW1VR0eHWltblUqllE6nlc1mK33EmkBnO7S2QWcbdLZRdZ2dsVwu5yS5XC43rZ+TyWRcQ0ODC4LABUHgJE1eE19raGhwmUzG08lnJzrbobUNOtugsw3LzsU+fs/IYZHJZFw0GnWRSGRKxK9ekUjERaNRfnFDorMdWtugsw0627DuXOzjd+Ccc8U+u7Fz507t3LlT//3vfyVJixcv1gsvvKA1a9YU+yM0MjKieDyuXC6npqamom83IZvNqq2tTRcuXFChULjm90ciEdXX12twcFDNzc0l399sRWc7tLZBZxt0tlGJzsU+fpf0Gou2tja9/PLLOnbsmI4dO6b7779fjz76qE6ePBnqkGHs3r1bY2NjRYWUpEKhoLGxMXV3d5f5ZLWFznZobYPONuhso5o7l/SMxeUkEglt375dTz31VFHfP51nLJxzSqVSOnXqlEo5dhAEWrhwofr6+nglchHobIfWNuhsg842KtW5LM9Y/F8XL17Unj17dP78ea1cufKK3zc+Pq6RkZEpV1hDQ0Pq7+8vKaR06T9Cf3+/hoeHQ9/3bEJnO7S2QWcbdLZR7Z1LHhYnTpxQY2Oj6urq9PTTT2v//v265ZZbrvj9XV1disfjk1d7e3vow46Ojoa+rSTl8/lp3X62oLMdWtugsw0626j2ziUPi0WLFqm3t1dHjhzRj3/8Y61fv17//Oc/r/j9W7duVS6Xm7wGBgZCH7axsTH0bSUpFotN6/azBZ3t0NoGnW3Q2Ua1dy55WMydO1c333yzVqxYoa6uLt12221Kp9NX/P66ujo1NTVNucJqaWlRMpks+W9DQRAomUwqkUiEvu/ZhM52aG2DzjbobKPaO0/7nTedcxofH/dxlmsKgkCdnZ2hbrtp0yZeFFQkOtuhtQ0626CzjWrvXNK/Cnnuuee0Zs0atbe3K5/Pa8+ePXr55ZeVyWT00EMPFfUzeB+LmYHOdmhtg8426GyjZt7H4tNPP9W6deu0aNEiPfDAA/rHP/5R0qjwobm5Wfv27VMQBIpErn78SCSiIAjU09PDL2yJ6GyH1jbobIPONqq687Te3zME688KOXDggKeTz050tkNrG3S2QWcblp1r+rNCJpw7d86l02mXTCanxEwmky6dTrtsNuvhxKCzHVrboLMNOtuw6lyWzwrxYbqvsbgc55yGh4eVz+cVi8WUSCR4EVAZ0NkOrW3Q2QadbZS7c7GP3zUxLAAAQHmV/S29AQAAvophAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8IZhAQAAvGFYAAAAbxgWAADAG4YFAADwhmEBAAC8YVgAAABvGBYAAMAbhgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8GZOpQ/gg3NOQ0NDGh0dVWNjo1paWhQEQaWPVXPobIfWNuhsg842qqXzjH7GIpvNKp1OK5VKqbW1VR0dHWptbVUqlVI6nVY2m630EWsCne3Q2gadbdDZRtV1dsZyuZyT5HK53LR+TiaTcQ0NDS4IAhcEgZM0eU18raGhwWUyGU8nn53obIfWNuhsg842LDsX+/g9I4dFJpNx0WjURSKRKRG/ekUiEReNRvnFDYnOdmhtg8426GzDunOxj9+Bc84V++xGV1eXenp69K9//Uv19fVatWqVfvGLX2jRokXF/giNjIwoHo8rl8upqamp6NtNyGazamtr04ULF1QoFK75/ZFIRPX19RocHFRzc3PJ9zdb0dkOrW3Q2QadbVSic7GP3yW9xuLQoUPasGGDjhw5ooMHD+rLL7/U6tWrdf78+VCHDGP37t0aGxsrKqQkFQoFjY2Nqbu7u8wnqy10tkNrG3S2QWcb1dy5pGcsvuqzzz7TDTfcoEOHDumee+4p6jbTecbCOadUKqVTp06plGMHQaCFCxeqr6+PVyIXgc52aG2DzjbobKNSncvyjMVX5XI5SVIikbji94yPj2tkZGTKFdbQ0JD6+/tLCild+o/Q39+v4eHh0Pc9m9DZDq1t0NkGnW1Ue+fQw8I5py1btuiuu+7SrbfeesXv6+rqUjwen7za29vD3qVGR0dD31aS8vn8tG4/W9DZDq1t0NkGnW1Ue+fQw2Ljxo364IMP9Ic//OGq37d161blcrnJa2BgIOxdqrGxMfRtJSkWi03r9rMFne3Q2gadbdDZRrV3DjUsOjs79dZbb+kvf/mL2trarvq9dXV1ampqmnKF1dLSomQyWfLfhoIgUDKZvOqfbPD/0dkOrW3Q2QadbVR755KGhXNOGzduVE9Pj9555x11dHSU61yXFQSBOjs7Q91206ZNvCioSHS2Q2sbdLZBZxvV3rmkfxXyk5/8RK+//rrefPPNKe9dEY/HVV9fX9TP4H0sZgY626G1DTrboLONmnkfi507dyqXy+nee+/VjTfeOHnt3bs31CHDaG5u1r59+xQEgSKRqx8/EokoCAL19PTwC1siOtuhtQ0626CzjaruPK339wzB+rNCDhw44OnksxOd7dDaBp1t0NmGZeea/qyQCefOnXPpdNolk8kpMZPJpEun0y6bzXo4Mehsh9Y26GyDzjasOpfls0J8mO5rLC7HOafh4WHl83nFYjElEgleBFQGdLZDaxt0tkFnG+XuXOzjd00MCwAAUF4mb+kNAADwfzEsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3syp9AF8cM5paGhIo6OjamxsVEtLi4IgqPSxag6d7dDaBp1t0NlGtXSe0c9YZLNZpdNppVIptba2qqOjQ62trUqlUkqn08pms5U+Yk2gsx1a26CzDTrbqLrOzlgul3OSXC6Xm9bPyWQyrqGhwQVB4IIgcJImr4mvNTQ0uEwm4+nksxOd7dDaBp1t0NmGZediH79n5LDIZDIuGo26SCQyJeJXr0gk4qLRKL+4IdHZDq1t0NkGnW1Ydy728TtwzrlSnuE4fPiwtm/fruPHj+v06dPav3+/1q5dW/TtR0ZGFI/Hlcvl1NTUVMpdS7r0lE9bW5suXLigQqFwze+PRCKqr6/X4OCgmpubS76/2YrOdmhtg8426GyjEp2Lffwu+TUW58+f12233aZXXnkl1MGma/fu3RobGysqpCQVCgWNjY2pu7u7zCerLXS2Q2sbdLZBZxvV3LnkZyym3DgITJ+xcM4plUrp1KlTKuXYQRBo4cKF6uvr45XIRaCzHVrboLMNOtuoVOeyPWNRqvHxcY2MjEy5whoaGlJ/f39JIaVL/xH6+/s1PDwc+r5nEzrbobUNOtugs41q71z2YdHV1aV4PD55tbe3h/5Zo6Oj0zpLPp+f1u1nCzrbobUNOtugs41q71z2YbF161blcrnJa2BgIPTPamxsnNZZYrHYtG4/W9DZDq1t0NkGnW1Ue+eyD4u6ujo1NTVNucJqaWlRMpks+W9DQRAomUwqkUiEvu/ZhM52aG2DzjbobKPaO8+od94MgkCdnZ2hbrtp0yZeFFQkOtuhtQ0626CzjWrvXPK/ChkdHdVHH30kSVq2bJl++ctf6r777lMikdCCBQuueXvex2JmoLMdWtugsw0626ip97E4duyYli1bpmXLlkmStmzZomXLlumFF14IddBSNTc3a9++fQqCQJHI1Y8fiUQUBIF6enr4hS0Rne3Q2gadbdDZRlV3ntb7e4Zg/VkhBw4c8HTy2YnOdmhtg8426GzDsnNNf1bIhHPnzrl0Ou2SyeSUmMlk0qXTaZfNZj2cGHS2Q2sbdLZBZxtWncv2WSHTNd3XWFyOc07Dw8PK5/OKxWJKJBK8CKgM6GyH1jbobIPONsrdudjH75oYFgAAoLyq5i29AQDA7MGwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3jAsAACANwwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeDOn0gfwwTmnoaEhjY6OqrGxUS0tLQqCoNLHqjl0tkNrG3S2QWcb1dJ5Rj9jkc1mlU6nlUql1Nraqo6ODrW2tiqVSimdTiubzVb6iDWBznZobYPONuhso+o6O2O5XM5Jcrlcblo/J5PJuIaGBhcEgQuCwEmavCa+1tDQ4DKZjKeTz050tkNrG3S2QWcblp2LffyekcMik8m4aDTqIpHIlIhfvSKRiItGo/zihkRnO7S2QWcbdLZh3bnYx+/AOedKfZbj1Vdf1fbt23X69GktXrxYO3bs0N13313UbUdGRhSPx5XL5dTU1FTqXSubzaqtrU0XLlxQoVC45vdHIhHV19drcHBQzc3NJd/fbEVnO7S2QWcbdLZRic7FPn6X/BqLvXv36plnntHzzz+v999/X3fffbfWrFmjjz/+ONRBS7V7926NjY0VFVKSCoWCxsbG1N3dXeaT1RY626G1DTrboLONqu5c6lMh3/3ud93TTz895Wvf+ta33E9/+lOvT6VcTqFQcMlk8mt/R7rWFQSBSyaTrlAolHyfsxGd7dDaBp1t0NlGpToX+/hd0jMWX3zxhY4fP67Vq1dP+frq1av17rvvXvY24+PjGhkZmXKFNTQ0pP7+frkS/3rjnFN/f7+Gh4dD3/dsQmc7tLZBZxt0tlHtnUsaFp9//rkuXryoefPmTfn6vHnzdObMmcvepqurS/F4fPJqb28PfdjR0dHQt5WkfD4/rdvPFnS2Q2sbdLZBZxvV3jnU+1h89Q03nHNXfBOOrVu3KpfLTV4DAwNh7lKS1NjYGPq2khSLxaZ1+9mCznZobYPONuhso9o7lzQsrr/+ekWj0a89O3H27NmvPYsxoa6uTk1NTVOusFpaWpRMJkt+J7EgCJRMJpVIJELf92xCZzu0tkFnG3S2Ue2dSxoWc+fO1fLly3Xw4MEpXz948KBWrVrl9WCXEwSBOjs7Q91206ZNvIVskehsh9Y26GyDzjaqvXPJ72Oxd+9erVu3Tr/+9a+1cuVK7dq1S6+99ppOnjypm2666Zq3530sZgY626G1DTrboLONmnofi8cee0w7duzQSy+9pNtvv12HDx/W22+/XdSo8KG5uVn79u1TEASKRK5+/EgkoiAI1NPTwy9siehsh9Y26GyDzjaqunOof8w6DdafFXLgwAFPJ5+d6GyH1jbobIPONiw71/RnhUw4d+6cS6fTLplMTomZTCZdOp122WzWw4lBZzu0tkFnG3S2YdW5rJ8VMh3TfY3F5TjnNDw8rHw+r1gspkQiwYuAyoDOdmhtg8426Gyj3J2LffyuiWEBAADKq2wv3gQAALgShgUAAPCGYQEAALxhWAAAAG8YFgAAwBuGBQAA8GaO9R1O/OvWkZER67sGAAAhTTxuX+tdKsyHRT6flyS1t7db3zUAAJimfD6veDx+xf+/+RtkFQoFffLJJ4rFYt7fEay9vV0DAwO88VYZ0dkOrW3Q2QadbZSzs3NO+Xxe8+fPv+oHn5k/YxGJRNTW1la2n9/U1MQvrQE626G1DTrboLONcnW+2jMVE3jxJgAA8IZhAQAAvKmZYVFXV6ef/exnqqurq/RRahqd7dDaBp1t0NlGNXQ2f/EmAACoXTXzjAUAAKg8hgUAAPCGYQEAALxhWAAAAG9qZli8+uqr6ujo0HXXXafly5frb3/7W6WPVFMOHz6sRx55RPPnz1cQBHrjjTcqfaSa1NXVpe985zuKxWK64YYbtHbtWv373/+u9LFqzs6dO7V06dLJNxFauXKl/vSnP1X6WDWvq6tLQRDomWeeqfRRas6LL76oIAimXN/4xjcqcpaaGBZ79+7VM888o+eff17vv/++7r77bq1Zs0Yff/xxpY9WM86fP6/bbrtNr7zySqWPUtMOHTqkDRs26MiRIzp48KC+/PJLrV69WufPn6/00WpKW1ubXn75ZR07dkzHjh3T/fffr0cffVQnT56s9NFq1tGjR7Vr1y4tXbq00kepWYsXL9bp06cnrxMnTlTkHDXxz03vvPNO3XHHHdq5c+fk17797W9r7dq16urqquDJalMQBNq/f7/Wrl1b6aPUvM8++0w33HCDDh06pHvuuafSx6lpiURC27dv11NPPVXpo9Sc0dFR3XHHHXr11Ve1bds23X777dqxY0elj1VTXnzxRb3xxhvq7e2t9FFm/jMWX3zxhY4fP67Vq1dP+frq1av17rvvVuhUgB+5XE7SpQc9lMfFixe1Z88enT9/XitXrqz0cWrShg0b9PDDD+vBBx+s9FFqWl9fn+bPn6+Ojg49/vjjOnXqVEXOYf4hZL59/vnnunjxoubNmzfl6/PmzdOZM2cqdCpg+pxz2rJli+666y7deuutlT5OzTlx4oRWrlyp//3vf2psbNT+/ft1yy23VPpYNWfPnj167733dPTo0Uofpabdeeed6u7u1je/+U19+umn2rZtm1atWqWTJ0+qpaXF9CwzflhM+OpHsDvnvH4sO2Bt48aN+uCDD/T3v/+90kepSYsWLVJvb6+y2az27dun9evX69ChQ4wLjwYGBrR582b9+c9/1nXXXVfp49S0NWvWTP7fS5Ys0cqVK5VMJrV7925t2bLF9Cwzflhcf/31ikajX3t24uzZs197FgOYKTo7O/XWW2/p8OHDamtrq/RxatLcuXN18803S5JWrFiho0ePKp1O6ze/+U2FT1Y7jh8/rrNnz2r58uWTX7t48aIOHz6sV155RePj44pGoxU8Ye1qaGjQkiVL1NfXZ37fM/41FnPnztXy5ct18ODBKV8/ePCgVq1aVaFTAeE457Rx40b19PTonXfeUUdHR6WPNGs45zQ+Pl7pY9SUBx54QCdOnFBvb+/ktWLFCj355JPq7e1lVJTR+Pi4PvzwQ914443m9z3jn7GQpC1btmjdunVasWKFVq5cqV27dunjjz/W008/Xemj1YzR0VF99NFHk//7P//5j3p7e5VIJLRgwYIKnqy2bNiwQa+//rrefPNNxWKxyWfi4vG46uvrK3y62vHcc89pzZo1am9vVz6f1549e/TXv/5VmUym0kerKbFY7GuvD2poaFBLSwuvG/Ls2Wef1SOPPKIFCxbo7Nmz2rZtm0ZGRrR+/Xrzs9TEsHjsscc0NDSkl156SadPn9att96qt99+WzfddFOlj1Yzjh07pvvuu2/yf0/8zW79+vX63e9+V6FT1Z6JfzJ97733Tvn6b3/7W/3oRz+yP1CN+vTTT7Vu3TqdPn1a8XhcS5cuVSaT0UMPPVTpowGhDA4O6oknntDnn3+u1tZWfe9739ORI0cq8jhYE+9jAQAAqsOMf40FAACoHgwLAADgDcMCAAB4w7AAAADeMCwAAIA3DAsAAOANwwIAAHjDsAAAAN4wLAAAgDcMCwAA4A3DAgAAeMOwAAAA3vw/oHVRlPXUHpIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from plot.plotter import plot_set_points\n",
    "\n",
    "\n",
    "fixer = isl.set(\"[N] -> {T[i,j]: N = 5}\")\n",
    "D_T = D_T.intersect(fixer)\n",
    "D_T = D_T.project_out(isl.dim_type.PARAM, 0, 1)\n",
    "plot_set_points(D_T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- The `isl.dim_type.param` means that you want to project out parameters, the two following numbers are the\n",
    "position of the first parameter and the number of consecutive parameters that will be projected out.\n",
    "For the sake of plotting, project out all of them. Note that if you forget to fix the parameter size, the set will become unbounded and cannot be plotted. -->\n",
    "`isl.dim_type.param`表示你想映射的参数，后面两个数据分别是第一个参数的位置和将要映射的连续参数的数量.\n",
    "\n",
    "为了plot出来，因此把所有维度都映射出来。注意，如果忘记*fix*参数大小，该集将变得unbounds，无法plot. \n",
    "\n",
    "<!-- Similarly to parameters, you can project out domain dimensions until you get at most two before plotting.\n",
    "Use `isl.dim_type.set` to get those. -->\n",
    "与参数类似，在plot之前可以投影domain dimensions获得两个维度.\n",
    "使用'isl.dim_type.set'来获取这些\n",
    "\n",
    "### Question\n",
    "<!-- In the following loop nest performing LU decomposition, -->\n",
    "下面执行LU分解的代码:\n",
    "\n",
    "```c\n",
    "for (i = 0; i < N; i++) {\n",
    "  for (j = 0; j < i; j++) {\n",
    "    for (k = 0; k < j; k++) {\n",
    "Sa:    A[i][j] -= A[i][k] * A[k][j];\n",
    "    }\n",
    "Sb: A[i][j] /= A[j][j];\n",
    "  }\n",
    "  for (j = i; j < N; j++) {\n",
    "    for (k = 0; k < i; k++) {\n",
    "Sc:   A[i][j] -= A[i][k] * A[k][j];\n",
    "    }\n",
    "  }\n",
    "}\n",
    "```\n",
    "\n",
    "<!-- 1. Define the union set with all iteration domains.\n",
    "1. Check if iteration domains of `Sa` and `Sc` overlap.\n",
    "1. Plot projections on the domain on (i,j), (j,k). -->\n",
    "1.  为所有的iteration domain定义union set.\n",
    "2.  检查iteration domain `Sa` 和 `Sc` 是否有overlap\n",
    "3.  绘制domain在(i,j)和(j,k)\n",
    "\n",
    "<!-- Hint: it is impossible to project out domain dimensions *directly* from the union set\n",
    "because it may have sets that live in different spaces.\n",
    "However, a union set may be constructed from sets, and you can take a union of union sets. -->\n",
    "Hint: 不能*直接*从union set中映射出domain dimension,因为他们可能存在不同的space中, 但是可以从set中做union得到他.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[N] -> { Sc[i, j, k] : i <= j < N and 0 <= k < i; Sa[i, j, k] : i < N and j < i and 0 <= k < j; Sb[i, j] : i < N and 0 <= j < i }\n"
     ]
    }
   ],
   "source": [
    "D_Sa = isl.set(\"[N] -> { Sa[i,j,k] :0 <= k < j < i < N }\")\n",
    "D_Sb = isl.set(\"[N] -> { Sb[i,j] :0 <= j < i < N }\")\n",
    "D_Sc = isl.set(\"[N] -> { Sc[i,j,k] :0 <= k < i <= j < N }\")\n",
    "D = isl.union_set(D_Sa).union(D_Sb).union(D_Sc)\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE 这里需要去掉标识名再做交集.\n",
    "D_Sa_ = isl.set(\"[N] -> { [i,j,k] :0 <= k < j < i < N }\")\n",
    "D_Sc_ = isl.set(\"[N] -> { [i,j,k] :0 <= k < i <= j < N }\")\n",
    "D_Sa_.intersect(D_Sc_).is_empty()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_Sa = D_Sa.intersect(isl.set(\"[N] -> { Sa[i,j,k]: N = 8}\"))\n",
    "D_Sb = D_Sb.intersect(isl.set(\"[N] -> { Sb[i,j]: N = 8}\"))\n",
    "D_Sc = D_Sc.intersect(isl.set(\"[N] -> { Sc[i,j,k]: N = 8}\"))\n",
    "D_Sa = D_Sa.project_out(isl.dim_type.PARAM, 0, 1)\n",
    "D_Sb = D_Sb.project_out(isl.dim_type.PARAM, 0, 1)\n",
    "D_Sc = D_Sc.project_out(isl.dim_type.PARAM, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ [i, j] : (i <= 7 and 0 < j < i) or (i > 0 and i <= j <= 7); Sb[i, j] : i <= 7 and 0 <= j < i }\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAALeRJREFUeJzt3W1wXOV5//HfWS0IZ/UUqUK1x0KYjQZCAIfgNuUhkzZQZjzEplKHNh2aOsO0UzfGBrudSZy+IC9UKW/SkTppmDHJMPbQlL7AENKWFc6EOGk7NEDqwQMpUdZZIQoIVVqtJIsRI53r/yJ/a1Cx7H2495zV3t/PzJlJxN461+9i95yL1dmzgZmZAAAAHEjEXQAAAKgfDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGeSUe8wDEO9+eabam5uVhAEUe8eAACUwcw0Pz+vLVu2KJFY/32JyAeLN998U93d3VHvFgAAODAxMaGtW7eu+88jHyyam5sl/aqwlpaWqHcPAADKMDc3p+7u7tXz+HoiHyzO/fmjpaWFwQIAgA3mYpcxcPEmAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnIv9USDWYmaanp7WwsKCmpiZ1dHR4d/OtMAw1NjamqakpdXZ2qre394I3MKk35Pc7v0QPfM8v0YOaORdaxAqFgkmyQqFQ8e/K5/M2PDxs6XTaJK1u6XTahoeHLZ/PV15wjcvlctbX12fJZHJND5LJpPX19Vkul4u7xKoiv9/5zeiB7/nN6EFU58Jiz98bdrDIZDKWSqUsCAILgmBNM8/9LJVKWSaTcVR57RkYGFiTe71tYGAg7lKrgvx+5zejB77nN6MHUZ4L63qwyGQy1tDQYIlE4oJPpEQiYQ0NDXU5XBT7YqrXFxX5/c5vRg98z29GD6I+F1ZlsOjp6Tlv0V/84hedF7aefD5vqVTqoo18f0NTqVRd/Vkkl8uV9GI6t9XL24Hk9zu/GT3wPb8ZPYjjXFjs+bukq1peeOEFvfXWW6vbiRMnJEn33HNPKb+mIkePHtXi4qLCMCzq8WEYanFxUceOHatyZdE5ePBgWesOHTrkuJJ4kN/v/BI98D2/RA9q+VwYmJmVu/jBBx/UP//zP2tsbKzoK0/n5ubU2tqqQqFQ8neFmJl6e3t15swZlVJ2EAS66qqrSqqzVoVhqMbGRi0vL5e8NplMamlpaUNfJU1+v/NL9MD3/BI9iOtcWOz5u+zOvvfee3rsscd03333XbDApaUlzc3NrdnKNT09rWw2W1IjpV/9S8hms5qZmSl737VibGysrBeTJC0vLyubzTquKFrk9zu/RA98zy/Rg1o/F5Y9WDz11FOanZ3VF77whQs+bmhoSK2tratbd3d3ubvUwsJC2WslaX5+vqL1tWBqaqqi9ZOTk44qiQf5/c4v0QPf80v0oNbPhWUPFt/+9re1c+dObdmy5YKPO3z4sAqFwuo2MTFR7i7V1NRU9lpJF/0O+Y2gs7OzovVdXV2OKokH+f3OL9ED3/NL9KDWz4VlXWMxPj6uq666SsePH9fdd99d0lqusaiM739bJL/f+SV64Ht+iR7U5TUWjz76qC6//HLddddd5SwvWxAE2r9/f1lrDxw4sOGHCklKJBLatWtXWWt37969oV9MEvl9zy/RA9/zS/Sg5s+FpX6OdWVlxa644gr70pe+VOpSM+M+Fi74/vlt8vud34we+J7fjB7UzX0sJOn73/++Xn/9dd13332lLnWira1NTzzxhIIguOjUmUgkFASBjh8/rra2tmgKjEBPT48GBgZKWjM4OKienp4qVRQt8vudX6IHvueX6EFNnwvLHl3KFPV3hYyOjjqqvPYUezvbwcHBuEutCvL7nd+MHvie34weRHkurOvvCjknn8/byMjIeb/RbWRkxGZnZx1UXNtyuZz19/ef91v9+vv7bXx8PO4Sq4r8fuc3owe+5zejB1GdC4s9f1d0581yVPKpkPWYmWZmZjQ/P6/m5ma1t7fXxYWapQjDUNlsVpOTk+rq6lI6nd7wFyiVgvx+55foge/5JXpQ7XNhsefvuhgsAABAdVX9lt4AAAD/F4MFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGeScRfggplpenpaCwsLampqUkdHh4IgiLusSIVhqLGxMU1NTamzs1O9vb1KJPyZG8nvd36JHvieX+JcUDP5LWKFQsEkWaFQqPh35fN5Gx4etnQ6bZJWt3Q6bcPDw5bP5ysvuMblcjnr6+uzZDK5pgfJZNL6+vosl8vFXWJVkd/v/Gb0wPf8ZpwLospf7Pl7ww4WmUzGUqmUBUFgQRCsaea5n6VSKctkMo4qrz0DAwNrcq+3DQwMxF1qVZDf7/xm9MD3/GacC6LMX9eDRSaTsYaGBkskEhd8MSUSCWtoaKjLJ1SxB5R6PbCQ3+/8ZvTA9/xmnAuizl+1weKNN96we++919rb223Tpk22fft2e/HFF50Xtp58Pm+pVOqijXx/Q1OpVF29FZbL5Uo6oJzb6uUtUfL7nd+MHvie34xzQRz5iz1/l3RlTz6f16233qpLLrlEzzzzjF599VV9/etfV1tbWym/piJHjx7V4uKiwjAs6vFhGGpxcVHHjh2rcmXROXjwYFnrDh065LiSeJDf7/wSPfA9v8S5oJbzB2ZmxT74y1/+sv793/9dP/7xj8ve4dzcnFpbW1UoFNTS0lLSWjNTb2+vzpw5oxLKVhAEuuqqqzQ2NrbhrxAOw1CNjY1aXl4ueW0ymdTS0tKGvlKc/H7nl+iB7/klzgVx5S/2/F3Ss+vpp5/Wjh07dM899+jyyy/XjTfeqEceeeSCa5aWljQ3N7dmK9f09LSy2WxJjZR+9S8hm81qZmam7H3XirGxsbIOKJK0vLysbDbruKJokd/v/BI98D2/xLmg1vOXNFicOXNGDz/8sHp7ezU6Oqq9e/fqwIEDF3xrZWhoSK2tratbd3d32cUuLCyUvVaS5ufnK1pfC6ampipaPzk56aiSeJDf7/wSPfA9v8S5oNbzlzRYhGGoT3ziExocHNSNN96oP//zP9ef/dmf6eGHH153zeHDh1UoFFa3iYmJsottamoqe60kNTc3V7S+FnR2dla0vqury1El8SC/3/kleuB7folzQa3nL2mw2Lx5s6699to1P/voRz+q119/fd01jY2NamlpWbOVq6OjQ+l0uuS/DQVBoHQ6rfb29rL3XSt6e3uVTJZ3w9RkMql0Ou24omiR3+/8Ej3wPb/EuaDW85c0WNx666167bXX1vzs5z//uXp6epwWtZ4gCLR///6y1h44cGBDX6xzTiKR0K5du8pau3v37g1/0Rb5/c4v0QPf80ucC2o+fymfYf3JT35iyWTS/uZv/sbGxsbsH/7hH+xDH/qQPfbYY84/B7se3z+7bMZn2Mnvd34zeuB7fjPOBbV8H4uSb5D1ve99z6677jprbGy0a665xo4cOVKVwi6k1LuNjY6Olr2vWlXqXfcGBwfjLtkp8vud34we+J7fjHNB1Pnr+pbeZsXfH73enkjvV+yBpR4PKGbk9z2/GT3wPb8Z54Io89f9YGH2q7eCRkZGzvuNbiMjIzY7O+ug4tqWy+Wsv7//vN9s2N/fb+Pj43GXWFXk9zu/GT3wPb8Z54Ko8hd7/i7pzpsuVHLnzfWYmWZmZjQ/P6/m5ma1t7dv+ItzShWGobLZrCYnJ9XV1aV0Ol0XF2kVi/x+55foge/5Jc4F1c5f7Pm7LgYLAABQXVW5pTcAAMCFMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcCYZdwEumJmmp6e1sLCgpqYmdXR0KAiCuMuKVBiGGhsb09TUlDo7O9Xb26tEwp+5kfx+55foAcdBelAz+S1ihULBJFmhUKj4d+XzeRseHrZ0Om2SVrd0Om3Dw8OWz+crL7jG5XI56+vrs2QyuaYHyWTS+vr6LJfLxV1iVZHf7/xm9IDjID2IKn+x5+8NO1hkMhlLpVIWBIEFQbCmmed+lkqlLJPJOKq89gwMDKzJvd42MDAQd6lVQX6/85vRA46D9CDK/HU9WGQyGWtoaLBEInHBg0kikbCGhoa6fEIVe0Ct1wMr+f3Ob0YPOA7Sg6jzV2WweOihhz5QcFdXV1UKW08+n7dUKnXRRr6/oalUqq7eCsvlciUdUM9t9fKWMPn9zm9GDzgO0oM48hd7/i75yqaPfexjeuutt1a306dPl/orKnL06FEtLi4qDMOiHh+GoRYXF3Xs2LEqVxadgwcPlrXu0KFDjiuJB/n9zi/RA46D9KCW8wdmZsU++Ktf/aqeeuopnTp1quwdzs3NqbW1VYVCQS0tLSWtNTP19vbqzJkzKqFsBUGgq666SmNjYxv+CuEwDNXY2Kjl5eWS1yaTSS0tLW3oK+XJ73d+iR5wHKQHceUv9vxd8qtrbGxMW7Zs0bZt2/S5z31OZ86cueDjl5aWNDc3t2Yr1/T0tLLZbEmNlH71LyGbzWpmZqbsfdeKsbGxsg6okrS8vKxsNuu4omiR3+/8Ej3gOEgPaj1/SYPFJz/5SR07dkyjo6N65JFH9Pbbb+uWW27R9PT0umuGhobU2tq6unV3d5dd7MLCQtlrJWl+fr6i9bVgamqqovWTk5OOKokH+f3OL9EDjoP0oNbzlzRY7Ny5U7//+7+v66+/XnfccYf+5V/+RdKv/taznsOHD6tQKKxuExMTZRfb1NRU9lpJam5urmh9Lejs7KxofVdXl6NK4kF+v/NL9IDjID2o9fwV/aExlUrp+uuv19jY2LqPaWxsVEtLy5qtXB0dHUqn0yX/bSgIAqXTabW3t5e971rR29urZLK8G6Ymk0ml02nHFUWL/H7nl+gBx0F6UOv5KxoslpaW9LOf/UybN292Vc8FBUGg/fv3l7X2wIEDG/pinXMSiYR27dpV1trdu3dv6IvWJPL7nl+iBxwH6UHN5y/lM6x/+Zd/aT/84Q/tzJkz9vzzz9tnP/tZa25uLumz4dzHonK+f4af/H7nN6MHHAfpQd3cx+KNN97QH/3RH+nqq69Wf3+/Lr30Uj3//PPq6ekp5ddUpK2tTU888YSCILjof3kkEgkFQaDjx4+rra0tmgIj0NPTo4GBgZLWDA4ORvrvqZrI73d+iR5wHKQHNZ2/7NGlTFF/V8jo6KijymtPsbc0HhwcjLvUqiC/3/nN6AHHQXoQZf66/q6Qc/L5vI2MjJz3G91GRkZsdnbWQcW1LZfLWX9//3m/2bG/v9/Gx8fjLrGqyO93fjN6wHGQHkSVv9jzd0l33nShkjtvrsfMNDMzo/n5eTU3N6u9vX3DX5xTqjAMlc1mNTk5qa6uLqXT6Q1/kVopyO93fokecBykB9XOX+z5uy4GCwAAUF1Vu6U3AADAehgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADiTjLsAF8xM09PTWlhYUFNTkzo6OhQEQdxlRSoMQ42NjWlqakqdnZ3q7e1VIuHP3Eh+v/NLHAd8zy/Rg5rJbxErFAomyQqFQsW/K5/P2/DwsKXTaZO0uqXTaRseHrZ8Pl95wTUul8tZX1+fJZPJNT1IJpPW19dnuVwu7hKrivx+5zfjOOB7fjN6EFX+Ys/fG3awyGQylkqlLAgCC4JgTTPP/SyVSlkmk3FUee0ZGBhYk3u9bWBgIO5Sq4L8fuc34zjge34zehBl/roeLDKZjDU0NFgikbjgATWRSFhDQ0NdPqGKPanU68mF/H7nN+M44Ht+M3oQdf5IBovBwUGTZA888IDzwtaTz+ctlUpdtJHvb2gqlaqrt8JyuVxJJ5VzW728LU5+v/ObcRzwPb8ZPYgjf7Hn77Kv7nrhhRd05MgR3XDDDeX+irIcPXpUi4uLCsOwqMeHYajFxUUdO3asypVF5+DBg2WtO3TokONK4kF+v/NLHAd8zy/Rg5rOX87UMj8/b729vXbixAn79Kc/Hdk7FmEYWjqd/sDfkS62BUFg6XTawjAseZ+1ZmVl5QMX6hW7JZNJW1lZiTtCRcjvd34zjgO+5zejB3Hlr+o7Fvv27dNdd92lO+6446KPXVpa0tzc3JqtXNPT08pmszKzktaZmbLZrGZmZsred60YGxvT8vJyWWuXl5eVzWYdVxQt8vudX+I44Ht+iR7Uev6SB4vHH39cP/3pTzU0NFTU44eGhtTa2rq6dXd3l1zkOQsLC2WvlaT5+fmK1teCqampitZPTk46qiQe5Pc7v8RxwPf8Ej2o9fwlDRYTExN64IEH9Nhjj+myyy4ras3hw4dVKBRWt4mJibIKlaSmpqay10pSc3NzRetrQWdnZ0Xru7q6HFUSD/L7nV/iOOB7foke1Hr+wEp4L+Wpp55SX1+fGhoaVn+2srKiIAiUSCS0tLS05p+dz9zcnFpbW1UoFNTS0lJSsWam3t5enTlzpqS3gIIg0FVXXaWxsbENfxe2MAzV2NhY1tvhyWRSS0tLG/qOjOT3O7/EccD3/BI9iCt/sefvko4wt99+u06fPq1Tp06tbjt27NC9996rU6dOXXSoqFQQBNq/f39Zaw8cOLChn0jnJBIJ7dq1q6y1u3fv3vAnFfL7nV/iOOB7foke1Hr+kt6xOJ/f/u3f1sc//nENDw8X9fhK3rGQpNnZWW3dulXvvvtuUR+zSSQS2rRpk9544w21tbWVvL9aND4+riuvvLLkdblcTj09Pe4Lihj5/c4vcRzwPb9ED+LIX5V3LGpBW1ubnnjiidU/v1xIIpFQEAQ6fvx4XTyRzunp6dHAwEBJawYHB+vmpEJ+v/NLHAd8zy/Rg5rOX9aHWSsQ9XeFjI6OOqq89hR7W+fBwcG4S60K8vud34zjgO/5zehBlPnr+rtCzsnn8zYyMnLeb3QbGRmx2dlZBxXXtlwuZ/39/ef9dsv+/n4bHx+Pu8SqIr/f+c04Dvie34weRJW/2PN3xddYlKrSayzOx8w0MzOj+fl5NTc3q729fcNfnFOqMAyVzWY1OTmprq4updPpurhQr1jk9zu/xHHA9/wSPah2/mLP33UxWAAAgOqq24s3AQBA7WKwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOBMMu4CXDAzTU9Pa2FhQU1NTero6FAQBHGXFakwDDU2NqapqSl1dnaqt7dXiYQ/c6Pv+XkN0APf80v0oGbyW8QKhYJJskKhUPHvyufzNjw8bOl02iStbul02oaHhy2fz1decI3L5XLW19dnyWRyTQ+SyaT19fVZLpeLu8Sq8j0/rwF64Ht+M3oQVf5iz98bdrDIZDKWSqUsCAILgmBNM8/9LJVKWSaTcVR57RkYGFiTe71tYGAg7lKrwvf8vAboge/5zehBlPnrerDIZDLW0NBgiUTigieURCJhDQ0NdfmEKvakWq8nV9/z8xqgB77nN6MHUeevymDxzW9+066//nprbm625uZm+63f+i3713/916oUtp58Pm+pVOqijXx/Q1OpVF29FZbL5Uo6qZ7b6uXPAr7n5zVAD3zPb0YP4shf7Pm7pKvbtm7dqq997Wt68cUX9eKLL+ozn/mM7r77br3yyiul/JqKHD16VIuLiwrDsKjHh2GoxcVFHTt2rMqVRefgwYNlrTt06JDjSuLhe35eA/TA9/wSPajp/GWPLv/fhz/8YfvWt77lfOI5nzAMLZ1Of+DvSBfbgiCwdDptYRiWvM9as7Ky8oELFYvdksmkraysxB2hIr7n5zVAD3zPb0YP4spflXcs3m9lZUWPP/64zp49q5tvvnndxy0tLWlubm7NVq7p6Wlls1mZWUnrzEzZbFYzMzNl77tWjI2NaXl5uay1y8vLymazjiuKlu/5eQ3QA9/zS/Sg1vOXPFicPn1aTU1Namxs1N69e/Xkk0/q2muvXffxQ0NDam1tXd26u7vLLnZhYaHstZI0Pz9f0fpaMDU1VdH6yclJR5XEw/f8vAboge/5JXpQ6/lLHiyuvvpqnTp1Ss8//7z+4i/+Qnv27NGrr7667uMPHz6sQqGwuk1MTJRdbFNTU9lrJam5ubmi9bWgs7OzovVdXV2OKomH7/l5DdAD3/NL9KDW85c8WFx66aX6yEc+oh07dmhoaEjbt2/XyMjIuo9vbGxUS0vLmq1cHR0dSqfTJd9JLAgCpdNptbe3l73vWtHb26tksrwbpiaTSaXTaccVRcv3/LwG6IHv+SV6UOv5K77nsZlpaWnJRS0XFQSB9u/fX9baAwcO1MWtXROJhHbt2lXW2t27d2/421z7np/XAD3wPb9ED2o+fylXhB4+fNh+9KMf2S9/+Ut7+eWX7Stf+YolEgl79tlnnV9Vuh7fP7tsxn0cfM/Pa4Ae+J7fjB7UzX0sJicn9fnPf15XX321br/9dv3nf/6nMpmMfvd3f7eUX1ORtrY2PfHEEwqC4KL/9ZlIJBQEgY4fP662trZoCoxAT0+PBgYGSlozODionp6eKlUULd/z8xqgB77nl+hBTecve3QpU9TfFTI6Ouqo8tpT7G2tBwcH4y61KnzPz2uAHvie34weRJm/rr8r5Jx8Pm8jIyPn/Ua3kZERm52ddVBxbcvlctbf33/eb/fs7++38fHxuEusKt/z8xqgB77nN6MHUeUv9vwdmJV4h40Kzc3NqbW1VYVCoaJPiLyfmWlmZkbz8/Nqbm5We3v7hr84p1RhGCqbzWpyclJdXV1Kp9Mb/kLFUvien9cAPfA9v0QPqp2/2PN3XQwWAACguoo9f/vzn3QAAKDqGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOJOMuwAXzEzT09NaWFhQU1OTOjo6FARB3GVFKgxDjY2NaWpqSp2dnert7VUi4c/c6PtzwPf8Ej3wPb9ED2ol/4Y+88zOzmpkZES9vb3q7OzUtm3bVk+qIyMjmp2djbvEqhsfH1d/f78aGxt1zTXX6FOf+pSuueYaNTY2qr+/X+Pj43GXWFW+Pwd8zy/RA9/zS/Sg5vJbxAqFgkmyQqFQ0e/JZDKWSqUsCAILgsAkrW7nfpZKpSyTyTiqvPYMDAysyb3eNjAwEHepVeH7c8D3/Gb0wPf8ZvQgyvzFnr835GCRyWSsoaHBEonEBU+oiUTCGhoa6vIJVexQUa/Dhe/PAd/zm9ED3/Ob0YOo81dlsBgcHLQdO3ZYU1OTdXZ22t13323//d//XZXC1pPP5y2VSl20ke9vaCqVsnw+X9b+alEulytpqDi35XK5uEt3wvfngO/5zeiB7/nN6EEc+Ys9f5d0jcXJkye1b98+Pf/88zpx4oSWl5d155136uzZs6X8moocPXpUi4uLCsOwqMeHYajFxUUdO3asypVF5+DBg2WtO3TokONK4uH7c8D3/BI98D2/RA9qOn/Zo4uZvfPOOybJTp486XziOZ8wDC2dTn/g70gX24IgsHQ6bWEYlrzPWrOysmLJZLKsdyySyaStrKzEHaEivj8HfM9vRg98z29GD+LKX5V3LP6vQqEgSWpvb1/3MUtLS5qbm1uzlWt6elrZbFZmVtI6M1M2m9XMzEzZ+64VY2NjWl5eLmvt8vKystms44qi5ftzwPf8Ej3wPb9ED2o9f9mDhZnp0KFDuu2223Tdddet+7ihoSG1traubt3d3eXuUgsLC2WvlaT5+fmK1teCqampitZPTk46qiQevj8HfM8v0QPf80v0oNbzlz1Y3H///Xr55Zf1j//4jxd83OHDh1UoFFa3iYmJcneppqamstdKUnNzc0Xra0FnZ2dF67u6uhxVEg/fnwO+55foge/5JXpQ6/nLGiz279+vp59+Ws8995y2bt16wcc2NjaqpaVlzVaujo4OpdPpku8kFgSB0un0Bf9ks1H09vYqmSzvhqnJZFLpdNpxRdHy/Tnge36JHvieX6IHtZ6/pMHCzHT//ffr+PHj+sEPfqBt27ZVq67zCoJA+/fvL2vtgQMH6uLWrolEQrt27Spr7e7duzf8bb59fw74nl+iB77nl+hBrecPrISrP774xS/qO9/5jr773e/q6quvXv15a2urNm3aVNTvmJubU2trqwqFQlnvXszOzmrr1q169913i/qYTSKR0KZNm/TGG2+ora2t5P3VovHxcV155ZUlr8vlcurp6XFfUMR8fw74nl+iB77nl+hBHPmLPn+X8lETrfMRlkcffdT5x1UupNS7jY2Ojpa9r1pV6p03BwcH4y7ZKd+fA77nN6MHvuc3owdR56/rW3qbFX9/9Hp7Ir1fscNFvQ0V5/j+HPA9vxk98D2/GT2IMn/dDxZmv7ql6cjIiKXT6TXNTKfTNjIyYrOzsw4qrm25XM76+/s/cNOsZDJp/f39Nj4+HneJVeX7c8D3/Gb0wPf8ZvQgqvzFnr9LusbChUqvsTgfM9PMzIzm5+fV3Nys9vb2DX9xTqnCMFQ2m9Xk5KS6urqUTqc3/IWapfD9OeB7foke+J5fogfVzl/s+bsuBgsAAFBdxZ6//flPWgAAUHUMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOAMgwUAAHCGwQIAADjDYAEAAJxhsAAAAM4wWAAAAGcYLAAAgDMMFgAAwBkGCwAA4AyDBQAAcIbBAgAAOMNgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcScZdgAtmpunpaS0sLKipqUkdHR0KgiDusiLlew/I73d+iR74nl+iB7WSf0O/YzE7O6uRkRH19vaqs7NT27ZtU2dnp3p7ezUyMqLZ2dm4S6w633tAfr/zS/TA9/wSPai5/BaxQqFgkqxQKFT0ezKZjKVSKQuCwIIgMEmr27mfpVIpy2QyjiqvPb73gPx+5zejB77nN6MHUeYv9vy9IQeLTCZjDQ0Nlkgk1jTx/26JRMIaGhrq8gnlew/I73d+M3rge34zehB1/qoNFidPnrTPfvaztnnzZpNkTz75ZFUKW08+n7dUKnXRRr6/oalUyvL5fFn7q0W+94D8fuc3owe+5zejB3HkL/b8XfI1FmfPntX27dv1jW98o9SlThw9elSLi4sKw7Cox4dhqMXFRR07dqzKlUXH9x6Q3+/8Ej3wPb9ED2o6f9mji1nk71iEYWjpdPoDf0e62BYEgaXTaQvDsOR91hrfe0B+v/Ob0QPf85vRg7jyV+0di1ItLS1pbm5uzVau6elpZbNZmVlJ68xM2WxWMzMzZe+7VvjeA/L7nV+iB77nl+hBreev+mAxNDSk1tbW1a27u7vs37WwsFBRLfPz8xWtrwW+94D8fueX6IHv+SV6UOv5qz5YHD58WIVCYXWbmJgo+3c1NTVVVEtzc3NF62uB7z0gv9/5JXrge36JHtR6/qoPFo2NjWppaVmzlaujo0PpdLrkO4kFQaB0Oq329vay910rfO8B+f3OL9ED3/NL9KDW82+oO28GQaD9+/eXtfbAgQN1cWtX33tAfr/zS/TA9/wSPaj1/IGVePXHwsKCfvGLX0iSbrzxRv3t3/6tfud3fkft7e264oorLrp+bm5Ora2tKhQKZb17MTs7q61bt+rdd98t6mM2iURCmzZt0htvvKG2traS91eLfO8B+f3OL9ED3/NL9CCO/EWfv0v9uMlzzz133o+x7Nmzx+nHVS6k1LuNjY6Olr2vWuV7D8jvd34zeuB7fjN6EHX+ur6lt1nx90evtyfS+/neA/L7nd+MHvie34weRJm/7gcLs1/d0nRkZMTS6fSaZqbTaRsZGbHZ2VkHFdc233tAfr/zm9ED3/Ob0YOo8hd7/i75GotKVXqNxfmYmWZmZjQ/P6/m5ma1t7dv+ItzSuV7D8jvd36JHvieX6IH1c5f7Pm7LgYLAABQXcWevzfUx00BAEBtY7AAAADOMFgAAABnGCwAAIAzDBYAAMAZBgsAAOBMMuodnvt069zcXNS7BgAAZTp33r7YXSoiHyzm5+clSd3d3VHvGgAAVGh+fl6tra3r/vPIb5AVhqHefPNNNTc3O78jWHd3tyYmJry98ZbvPSC/3/kleuB7fokeVDO/mWl+fl5btmxRIrH+lRSRv2ORSCS0devWqv3+lpYWL59M7+d7D8jvd36JHvieX6IH1cp/oXcqzuHiTQAA4AyDBQAAcKZuBovGxkY99NBDamxsjLuU2PjeA/L7nV+iB77nl+hBLeSP/OJNAABQv+rmHQsAABA/BgsAAOAMgwUAAHCGwQIAADhTN4PFN7/5TW3btk2XXXaZbrrpJv34xz+Ou6TI/OhHP9KuXbu0ZcsWBUGgp556Ku6SIjU0NKTf+I3fUHNzsy6//HL93u/9nl577bW4y4rMww8/rBtuuGH1hjg333yznnnmmbjLis3Q0JCCINCDDz4YdymR+epXv6ogCNZsv/7rvx53WZH6n//5H/3xH/+xOjo69KEPfUgf//jH9dJLL8VdVmSuvPLKDzwHgiDQvn37Iq+lLgaLf/qnf9KDDz6ov/7rv9Z//dd/6VOf+pR27typ119/Pe7SInH27Flt375d3/jGN+IuJRYnT57Uvn379Pzzz+vEiRNaXl7WnXfeqbNnz8ZdWiS2bt2qr33ta3rxxRf14osv6jOf+YzuvvtuvfLKK3GXFrkXXnhBR44c0Q033BB3KZH72Mc+prfeemt1O336dNwlRSafz+vWW2/VJZdcomeeeUavvvqqvv71r6utrS3u0iLzwgsvrPn3f+LECUnSPffcE30xVgd+8zd/0/bu3bvmZ9dcc419+ctfjqmi+EiyJ598Mu4yYvXOO++YJDt58mTcpcTmwx/+sH3rW9+Ku4xIzc/PW29vr504ccI+/elP2wMPPBB3SZF56KGHbPv27XGXEZsvfelLdtttt8VdRk154IEHLJ1OWxiGke97w79j8d577+mll17SnXfeuebnd955p/7jP/4jpqoQp0KhIElqb2+PuZLorays6PHHH9fZs2d18803x11OpPbt26e77rpLd9xxR9ylxGJsbExbtmzRtm3b9LnPfU5nzpyJu6TIPP3009qxY4fuueceXX755brxxhv1yCOPxF1WbN577z099thjuu+++5x+2WexNvxg8b//+79aWVlRV1fXmp93dXXp7bffjqkqxMXMdOjQId1222267rrr4i4nMqdPn1ZTU5MaGxu1d+9ePfnkk7r22mvjLisyjz/+uH76059qaGgo7lJi8clPflLHjh3T6OioHnnkEb399tu65ZZbND09HXdpkThz5owefvhh9fb2anR0VHv37tWBAwd07NixuEuLxVNPPaXZ2Vl94QtfiGX/kX+7abX836nMzGKZ1BCv+++/Xy+//LL+7d/+Le5SInX11Vfr1KlTmp2d1RNPPKE9e/bo5MmTXgwXExMTeuCBB/Tss8/qsssui7ucWOzcuXP1f19//fW6+eablU6ndfToUR06dCjGyqIRhqF27NihwcFBSdKNN96oV155RQ8//LD+5E/+JObqovftb39bO3fu1JYtW2LZ/4Z/x+LXfu3X1NDQ8IF3J955550PvIuB+rZ//349/fTTeu6557R169a4y4nUpZdeqo985CPasWOHhoaGtH37do2MjMRdViReeuklvfPOO7rpppuUTCaVTCZ18uRJ/d3f/Z2SyaRWVlbiLjFyqVRK119/vcbGxuIuJRKbN2/+wBD90Y9+1JsL+N9vfHxc3//+9/Wnf/qnsdWw4QeLSy+9VDfddNPqFbDnnDhxQrfccktMVSFKZqb7779fx48f1w9+8ANt27Yt7pJiZ2ZaWlqKu4xI3H777Tp9+rROnTq1uu3YsUP33nuvTp06pYaGhrhLjNzS0pJ+9rOfafPmzXGXEolbb731Ax8x//nPf66enp6YKorPo48+qssvv1x33XVXbDXUxZ9CDh06pM9//vPasWOHbr75Zh05ckSvv/669u7dG3dpkVhYWNAvfvGL1f//y1/+UqdOnVJ7e7uuuOKKGCuLxr59+/Sd73xH3/3ud9Xc3Lz67lVra6s2bdoUc3XV95WvfEU7d+5Ud3e35ufn9fjjj+uHP/yhMplM3KVForm5+QPX06RSKXV0dHhznc1f/dVfadeuXbriiiv0zjvvaGBgQHNzc9qzZ0/cpUXi4MGDuuWWWzQ4OKg/+IM/0E9+8hMdOXJER44cibu0SIVhqEcffVR79uxRMhnj6T3yz6FUyd///d9bT0+PXXrppfaJT3zCq48aPvfccybpA9uePXviLi0S58suyR599NG4S4vEfffdt/rc7+zstNtvv92effbZuMuKlW8fN/3DP/xD27x5s11yySW2ZcsW6+/vt1deeSXusiL1ve99z6677jprbGy0a665xo4cORJ3SZEbHR01Sfbaa6/FWgdfmw4AAJzZ8NdYAACA2sFgAQAAnGGwAAAAzjBYAAAAZxgsAACAMwwWAADAGQYLAADgDIMFAABwhsECAAA4w2ABAACcYbAAAADOMFgAAABn/h8JsrOjZ9I0fwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_Sa_ = D_Sa.project_out(isl.dim_type.SET, 2, 1)\n",
    "D_Sc_ = D_Sc.project_out(isl.dim_type.SET, 2, 1)\n",
    "D_ij = isl.union_set(D_Sa_).union(D_Sb).union(D_Sc_)\n",
    "print(D_ij)\n",
    "plot_set_points(D_ij)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ [j, k] : j <= 7 and 0 <= k < j }\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAJblJREFUeJzt3XtsW3f9//HXcbJlaeI4xLijIdlIjOnYldIO1G3wRWy0qqKJiYibxtZpgCjrkkKFBB0SFwHLJIQUD0ZFB0KNgJU/ko4haEwRI9w0WDeiVRsaIUE06TbKktqx6yiI+vP7Y0p+y9qkn2N/fEueD+lIzPXH560XTc6rx8fHnjHGCAAAwIFAuQcAAACrB8UCAAA4Q7EAAADOUCwAAIAzFAsAAOAMxQIAADhDsQAAAM5QLAAAgDO1pd5hLpfTCy+8oGAwKM/zSr17AACQB2OM0um0WltbFQgsf16i5MXihRdeUHt7e6l3CwAAHJicnFRbW9uyf17yYhEMBiW9MlhTU1Opdw8AAPIwOzur9vb2xeP4ckpeLBbe/mhqaqJYAABQZS50GQMXbwIAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZ0r+qZBiMMZoenpamUxGjY2NCofD3HxrGWRlj6z8IS97ZGWPrOxVSlZVfcYimUwqHo8rFospEomoo6NDkUhEsVhM8XhcyWSy3CNWDLKyR1b+kJc9srJHVvYqLitTYqlUykgyqVSqoNcZHh42DQ0NxvM843mekbS4LTzW0NBghoeHHU1evcjKHln5Q172yMoeWdkrZVa2x++qLBbDw8OmpqbGBAKBJSG+dgsEAqampmZN/+UjK3tk5Q952SMre2Rlr9RZ2R6/PWOMsT278aY3vUn/+te/znn8nnvu0UMPPWT1GrOzswqFQkqlUnndeTOZTKqtrU1zc3PK5XIXfH4gEFB9fb2mpqbU3Nzse3/VjKzskZU/5GWPrOyRlb1yZGV7/PZ1jcWTTz6pF198cXE7evSoJOmDH/xgXkPm4+DBg8pms1ZBSq98m2o2m9XAwECRJ6s8ZGWPrPwhL3tkZY+s7FV0VoWcFtmzZ4+JRqMml8s5P5VyPrlczkSj0XPeR7rQ5nme7zmrHVnZIyt/yMseWdkjK3vlyqoob4W82n//+1+1trZq7969uu+++5Z93vz8vObn5xf/e+Hb0fJ5K+Tll19WJBLJZ9zF9eFwOO/11YSs7JGVP+Rlj6zskZW9cmVVlLdCXu3RRx9VMpnUXXfdteLz+vr6FAqFFrf29vZ8d6lMJpP3WklKp9MFra8mZGWPrPwhL3tkZY+s7FV6VnkXix/84AfasWOHWltbV3zevn37lEqlFrfJycl8d6nGxsa810q64HfIryZkZY+s/CEve2Rlj6zsVXpWeRWLf/3rX/r1r3+tT3ziExd8bl1dnZqampZs+QqHw4pGo77vJOZ5nqLRqFpaWvLed7UhK3tk5Q952SMre2Rlr9KzyqtY/PCHP9T69evV1dXlep4VeZ6nnp6evNb29vauqdvAkpU9svKHvOyRlT2yslfpWfm+eDOXy6mjo0Mf/ehH9cADD/jeIfexKB2yskdW/pCXPbKyR1b2Vs19LCTp17/+tU6cOKG77747r8EK1dzcrMHBQXmep0Bg5fEDgYA8z9PQ0NCa+0snkZUfZOUPedkjK3tkZa+is8rrw6wFKPV3hSQSCUeTVy+yskdW/pCXPbKyR1b2SpnVqv6ukAWnT5828XjcRKPRJWFGo1ETj8dNMpl0MPHqQFb2yMof8rJHVvbIyl6psir6DbLyVeg1FudjjNHMzIzS6bSCwaBaWlrW1IU8fpCVPbLyh7zskZU9srJX7Kxsj9+rolgAAIDiKvqdNwEAAF6LYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJypLfcALhhjND09rUwmo8bGRoXDYXmeV+6xKhJZ2SMrf8jLHlnZI6vqU9VnLJLJpOLxuGKxmCKRiDo6OhSJRBSLxRSPx5VMJss9YsUgK3tk5Q952SMre2RVxUyJpVIpI8mkUqmCXmd4eNg0NDQYz/OM53lG0uK28FhDQ4MZHh52NHn1Iit7ZOUPedkjK3tkVZlsj99VWSyGh4dNTU2NCQQCS/7CvXYLBAKmpqZmTf/lIyt7ZOUPedkjK3tkVblsj9+eMcb4OcNx8uRJff7zn9eRI0c0Nzent7zlLfrBD36gzZs3W62fnZ1VKBRSKpVSU1OTn11LeuX0WFtbm+bm5pTL5S74/EAgoPr6ek1NTam5udn3/qoZWdkjK3/Iyx5Z2SOrymZ7/PZ1jcXp06d144036qKLLtKRI0f03HPP6Vvf+lZJ/w89ePCgstms1V86ScrlcspmsxoYGCjyZJWHrOyRlT/kZY+s7JHV6uDrjMUXvvAF/fGPf9Tvf//7vHdYyBkLY4xisZgmJibk50SL53nq7OzU2NjYmrmamKzskZU/5GWPrOyRVeWzPX77KhZXXnmltm/frqmpKY2MjOiNb3yj7rnnHn3yk59cds38/Lzm5+eXDNbe3p5XsXj55ZcViUR8rXnt+nA4nPf6akJW9sjKH/KyR1b2yKryFeWtkImJCe3fv1+xWEyJREK7du1Sb2/viqeh+vr6FAqFFrf29nY/u1wik8nkvVaS0ul0QeurCVnZIyt/yMseWdkjq9XDV7HI5XJ6+9vfrvvvv1+bNm3Spz71KX3yk5/U/v37l12zb98+pVKpxW1ycjLvYRsbG/NeK0nBYLCg9dWErOyRlT/kZY+s7JHV6uGrWGzYsEFXXnnlksfe+ta36sSJE8uuqaurU1NT05ItX+FwWNFo1Pf7aJ7nKRqNqqWlJe99VxuyskdW/pCXPbKyR1arh69iceONN+r5559f8tjf//53XX755U6HWo7neerp6clrbW9v75q6sIes7JGVP+Rlj6zskdUq4ufmGH/5y19MbW2t+cY3vmHGxsbMj3/8Y7Nu3Trzox/9yPkNNpZz+vRp09DQcMGbp+hVN1FpaGgwp0+fzmt/1Yys7JGVP+Rlj6zskVVlsz1++zpjcf311+vw4cN65JFHdPXVV+trX/ua+vv7dfvtt/t5mYI0NzdrcHBQnucpEFh5/EAgIM/zNDQ0tCZvnkJW9sjKH/KyR1b2yGqVKFHRWVTq7wpJJBKOJq9eZGWPrPwhL3tkZY+sKtOq/q6QBadPnzbxeNxEo9Elf/Gi0aiJx+MmmUw6mHh1ICt7ZOUPedkjK3tkVXmK9l0hhSr0u0LOxxijmZkZpdNpBYNBtbS0cCHPMsjKHln5Q172yMoeWVWOotx504ViFAsAAFBcRbnzJgAAwEooFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwJnacg/ggjFG09PTymQyamxsVDgclud55R6rIpGVPbLyh7zskRVWs6o+Y5FMJhWPxxWLxRSJRNTR0aFIJKJYLKZ4PK5kMlnuESsGWdkjK3/Iyx5ZYU0wJZZKpYwkk0qlCnqd4eFh09DQYDzPM57nGUmL28JjDQ0NZnh42NHk1Yus7JGVP+Rlj6xQ7WyP31VZLIaHh01NTY0JBAJLfjhfuwUCAVNTU7Omf1DJyh5Z+UNe9sgKq4Ht8dszxhjbsxtf+cpX9NWvfnXJY5deeqleeukl25fQ7OysQqGQUqmUmpqarNctSCaTamtr09zcnHK53AWfHwgEVF9fr6mpKTU3N/veXzUjK3tk5Q952SMrrBa2x2/f11hcddVVevHFFxe348ePFzSoXwcPHlQ2m7X6AZWkXC6nbDargYGBIk9WecjKHln5Q172yAprje8zFo8++qhGR0fz3mEhZyyMMYrFYpqYmJCPseV5njo7OzU2NrZmrrwmK3tk5Q952SMrrCZFO2MxNjam1tZWdXR06CMf+YgmJiZWfP78/LxmZ2eXbPmanp7W+Pi4rx9Q6ZUf7vHxcc3MzOS972pDVvbIyh/yskdWWIt8FYt3vvOdGhgYUCKR0MMPP6yXXnpJN9xwg6anp5dd09fXp1AotLi1t7fnPWwmk8l7rSSl0+mC1lcTsrJHVv6Qlz2ywlrkq1js2LFD3d3duuaaa3TLLbfoF7/4haRX3kNczr59+5RKpRa3ycnJvIdtbGzMe60kBYPBgtZXE7KyR1b+kJc9ssJaVNANshoaGnTNNddobGxs2efU1dWpqalpyZavcDisaDTq+z1Hz/MUjUbV0tKS976rDVnZIyt/yMseWWEtKqhYzM/P629/+5s2bNjgap4VeZ6nnp6evNb29vauqYugyMoeWflDXvbICmuRr0+FfO5zn9Ott96qyy67TKdOndLXv/51jYyM6Pjx47r88sutXoP7WJQOWdkjK3/Iyx5ZYbUoyqdCpqam9NGPflQbN27UBz7wAV188cV64oknrEuFC83NzRocHJTneQoEVh4/EAjI8zwNDQ2tyR9QsrJHVv6Qlz2ywprj9H6fFkr9XSGJRMLR5NWLrOyRlT/kZY+sUO1W9XeFLDh9+rSJx+MmGo0u+SGNRqMmHo+bZDLpYOLVgazskZU/5GWPrFDNivJdIS4Ueo3F+RhjNDMzo3Q6rWAwqJaWFi56WgZZ2SMrf8jLHlmhGtkev1dFsQAAAMVVtFt6AwAALIdiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKkt9wAuGGM0PT2tTCajxsZGhcNheZ5X7rEqElnZIyt/yAuAVOVnLJLJpOLxuGKxmCKRiDo6OhSJRBSLxRSPx5VMJss9YsUgK3tk5Q95AVjClFgqlTKSTCqVKuh1hoeHTUNDg/E8z3ieZyQtbguPNTQ0mOHhYUeTVy+yskdW/pAXsHbYHr+r8oxFIpFQV1eX5ubmZIyRMWbJny88Njc3p66uLiUSiTJNWn5kZY+s/CEvAOfjmdf+NvChr69P9913n/bs2aP+/n6rNbOzswqFQkqlUmpqavK9z2Qyqba2Ns3NzSmXy13w+YFAQPX19ZqamlJzc7Pv/VUzsrJHVv6QF7D22B6/8z5j8eSTT+rAgQO69tpr832JvBw8eFDZbNbql5kk5XI5ZbNZDQwMFHmyykNW9sjKH/ICsJy8ikUmk9Htt9+uhx9+WK973etcz7QsY4y+/e1v57X2wQcfPOdU7WpGVvbIyh/yArCSvIrF7t271dXVpVtuueWCz52fn9fs7OySLV/T09MaHx/3/YvJGKPx8XHNzMzkve9qQ1b2yMof8gKwEt/F4tChQ3r66afV19dn9fy+vj6FQqHFrb293feQCzKZTN5rJSmdThe0vpqQlT2y8oe8AKzEV7GYnJzUnj179KMf/UiXXHKJ1Zp9+/YplUotbpOTk3kNKkmNjY15r5WkYDBY0PpqQlb2yMof8gKwEl/F4qmnntKpU6e0efNm1dbWqra2ViMjI3rwwQdVW1urs2fPnrOmrq5OTU1NS7Z8hcNhRaNR33fz8zxP0WhULS0tee+72pCVPbLyh7wArMRXsbj55pt1/PhxjY6OLm5btmzR7bffrtHRUdXU1BRrTkmv/GLq6enJa21vb++aur0wWdkjK3/IC8BKCrqPhSS95z3v0dve9jbuY1GByMoeWflDXsDaU/T7WJRLc3OzBgcH5XmeAoGVxw8EAvI8T0NDQ2vylxlZ2SMrf8gLwLIc3D7cl1J/V0gikXA0efUiK3tk5Q95AWvHqv6uEEnavn27pqam1N/fr87OziV/1tnZqf7+fp08eVLbtm0r04SVg6zskZU/5AXgtQq+xsKvQq+xOB9jjGZmZpROpxUMBtXS0sIFYssgK3tk5Q95Aaub7fF7VRQLAABQXKv24k0AAFC5KBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMCZ2nIP4IIxRtPT08pkMmpsbFQ4HJbneeUeqyKRlT2yAgD/qvqMRTKZVDweVywWUyQSUUdHhyKRiGKxmOLxuJLJZLlHrBhkZY+sACB/njHGlHKHs7OzCoVCSqVSampqyvt1EomEuru7lc1mJb3yr8sFC/+qXLdunQYHB7V9+/bChq5yZGWPrADg/GyP31V5xiKRSKirq0tzc3Myxui13Wjhsbm5OXV1dSmRSJRp0vIjK3tkBQCF81Us9u/fr2uvvVZNTU1qamrS1q1bdeTIkWLNdl7JZFLd3d0yxiiXy6343FwuJ2OMuru71+Tpa7KyR1YA4IavYtHW1qYHHnhAx44d07Fjx/Te975X73//+/Xss88Wa75zHDx4UNls9oK//Bfkcjlls1kNDAwUebLKQ1b2yAoA3Cj4GouWlhZ985vf1Mc//nGr5xdyjYUxRrFYTBMTE+ecpl6J53nq7OzU2NjYmrmqn6zskRUAXFjRr7E4e/asDh06pDNnzmjr1q3LPm9+fl6zs7NLtnxNT09rfHzc1y9/6ZUDx/j4uGZmZvLed7UhK3tkBQDu+C4Wx48fV2Njo+rq6rRr1y4dPnxYV1555bLP7+vrUygUWtza29vzHjaTyeS9VpLS6XRB66sJWdkjKwBwx3ex2Lhxo0ZHR/XEE0/o05/+tHbu3Knnnntu2efv27dPqVRqcZucnMx72MbGxrzXSlIwGCxofTUhK3tkBQDuFHyNxS233KJoNKrvfe97Vs/nGovSICt7ZAUAF1ay+1gYYzQ/P1/oy1jxPE89PT15re3t7V1Tv/zJyh5ZAYA7vs5Y3HfffdqxY4fa29uVTqd16NAhPfDAAxoeHtb73vc+q9co9M6byWRSbW1tmpubs/poYCAQUH19vaamptTc3Ox7f9WMrOyRFQCsrChnLP7973/rjjvu0MaNG3XzzTfrz3/+s69S4UJzc7MGBwfleZ4CgZXHDwQC8jxPQ0NDa/KXP1nZIysAcGPVf1fI0NCQtm3bVtjQVY6s7JEVAJzfqv6uEEnavn27pqam1N/fr87OziV/1tnZqf7+fp08eZJf/iIrP8gKAApTtWcsXs0Yo5mZGaXTaQWDQbW0tHBB3TLIyh5ZAcD/Z3v8ri3hTEXjeZ7C4bDC4XC5R6l4ZGWPrADAv6p9KwQAAFQeigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCGYgEAAJyhWAAAAGcoFgAAwBmKBQAAcIZiAQAAnKFYAAAAZygWAADAGYoFAABwhmIBAACcoVgAAABnKBYAAMAZigUAAHCmttwDuGCM0fT0tDKZjBobGxUOh+V5XrnHqkhkBQAopqo+Y5FMJhWPxxWLxRSJRNTR0aFIJKJYLKZ4PK5kMlnuESsGWQEASsEzxphS7nB2dlahUEipVEpNTU15v04ikVB3d7ey2aykV/4lvmDhX+Dr1q3T4OCgtm/fXtjQVY6sAACFsj1+V+UZi0Qioa6uLs3NzckYo9d2o4XH5ubm1NXVpUQiUaZJy4+sAACl5KtY9PX16frrr1cwGNT69et122236fnnny/WbOeVTCbV3d0tY4xyudyKz83lcjLGqLu7e02e6icrAECp+SoWIyMj2r17t5544gkdPXpU//vf/7Rt2zadOXOmWPOd4+DBg8pmsxc8UC7I5XLKZrMaGBgo8mSVh6wAAKVW0DUW//nPf7R+/XqNjIzo3e9+t9WaQq6xMMYoFotpYmLinFP6K/E8T52dnRobG1szn4AgKwCASyW5xiKVSkmSWlpaln3O/Py8Zmdnl2z5mp6e1vj4uK8DpfTKQXZ8fFwzMzN577vakBUAoBzyLhbGGO3du1c33XSTrr766mWf19fXp1AotLi1t7fnu0tlMpm810pSOp0uaH01ISsAQDnkXSzuvfdePfPMM3rkkUdWfN6+ffuUSqUWt8nJyXx3qcbGxrzXSlIwGCxofTUhKwBAOeRVLHp6evTYY4/p8ccfV1tb24rPraurU1NT05ItX+FwWNFo1Pd7/57nKRqNrviWzWpDVgCAcvBVLIwxuvfeezU0NKTf/OY36ujoKNZc5+V5nnp6evJa29vbu6YuRiQrAEA5+PpUyD333KOf/OQn+tnPfqaNGzcuPh4KhVRfX2/1GoXeeTOZTKqtrU1zc3NWH6MMBAKqr6/X1NSUmpubfe+vmpEVAMCVonwqZP/+/UqlUnrPe96jDRs2LG4//elPCx7YVnNzswYHB+V5ngKBlccPBALyPE9DQ0Nr8kBJVgCAUvP9Vsj5trvuuqtI453f9u3b9Ytf/EL19fXyPO+c0/YLj9XX1+uXv/yltm3bVtL5KglZAQBKqSq/K0R65YA5NTWl/v5+dXZ2Lvmzzs5O9ff36+TJkxwoRVYAgNKp2m83fTVjjGZmZpROpxUMBtXS0sLFh8sgKwBAPmyP37UlnKloPM9TOBxWOBwu9ygVj6wAAMVUtW+FAACAykOxAAAAzlAsAACAMxQLAADgDMUCAAA4Q7EAAADOUCwAAIAzFAsAAOAMxQIAADhDsQAAAM5QLAAAgDMUCwAA4AzFAgAAOEOxAAAAzlAsAACAMxQLAADgDMUCAAA4Q7EAAADOUCwAAIAzFAsAAOAMxQIAADhDsQAAAM5QLAAAgDMUCwAA4AzFAgAAOEOxAAAAzlAsAACAMxQLAADgDMUCAAA4Q7EAAADOUCwAAIAzFAsAAOAMxQIAADhDsQAAAM5QLAAAgDMUCwAA4AzFAgAAOEOxAAAAztSWewAXjDGanp5WJpNRY2OjwuGwPM8r91gAAKw5VX3GIplMKh6PKxaLKRKJqKOjQ5FIRLFYTPF4XMlkstwjAgCwpnjGGFPKHc7OzioUCimVSqmpqSnv10kkEuru7lY2m5X0ylmLBQtnK9atW6fBwUFt3769sKEBAFjjbI/fVXnGIpFIqKurS3NzczLG6LXdaOGxubk5dXV1KZFIlGlSAADWFt/F4ne/+51uvfVWtba2yvM8Pfroo0UYa3nJZFLd3d0yxiiXy6343FwuJ2OMuru7eVsEAIAS8F0szpw5o+uuu07f+c53ijHPBR08eFDZbPaCpWJBLpdTNpvVwMBAkScDAAAFXWPheZ4OHz6s2267zXpNIddYGGMUi8U0MTFxztsfF5qzs7NTY2NjfFoEAIA82B6/i/5x0/n5ec3Pzy8ZLF/T09MaHx/3vc4Yo/Hxcc3MzCgcDue9fwAAsLKiX7zZ19enUCi0uLW3t+f9WplMpqBZ0ul0QesBAMDKil4s9u3bp1QqtbhNTk7m/VqNjY0FzRIMBgtaDwAAVlb0t0Lq6upUV1fn5LXC4bCi0Wje11i0tLQ4mQMAAJxfVd3HwvM89fT05LW2t7eXCzcBACgy38Uik8lodHRUo6OjkqR//vOfGh0d1YkTJ1zPdl47d+7UunXrFAjYjR4IBLRu3TrdeeedRZ4MAAD4LhbHjh3Tpk2btGnTJknS3r17tWnTJn3pS19yPtz5NDc3a3BwUJ7nXbBcBAIBeZ6noaEhNTc3l2Q+AADWslX/XSFDQ0Patm1bYUMDALDGrervCpGk7du3a2pqSv39/ers7FzyZ52dnerv79fJkycpFQAAlFDVnrF4NWOMZmZmlE6nFQwG1dLSwoWaAAA4VDF33iwFz/MUDoe5qyYAAGVWtW+FAACAykOxAAAAzlAsAACAMxQLAADgDMUCAAA4Q7EAAADOlPzjpgu3zZidnS31rgEAQJ4WjtsXuv1VyYtFOp2WJLW3t5d61wAAoEDpdFqhUGjZPy/5nTdzuZxeeOEFBYNBp3fHnJ2dVXt7uyYnJ53d0XO1Iit7ZOUPedkjK3tkZa+YWRljlE6n1drauuKXgJb8jEUgEFBbW1vRXr+pqYm/eJbIyh5Z+UNe9sjKHlnZK1ZWK52pWMDFmwAAwBmKBQAAcGbVFIu6ujp9+ctfVl1dXblHqXhkZY+s/CEve2Rlj6zsVUJWJb94EwAArF6r5owFAAAoP4oFAABwhmIBAACcoVgAAABnVk2x+O53v6uOjg5dcskl2rx5s37/+9+Xe6SK87vf/U633nqrWltb5XmeHn300XKPVLH6+vp0/fXXKxgMav369brtttv0/PPPl3usirR//35de+21izfk2bp1q44cOVLusapCX1+fPM/TZz7zmXKPUpG+8pWvyPO8Jdsb3vCGco9VsU6ePKmPfexjCofDWrdund72trfpqaeeKvkcq6JY/PSnP9VnPvMZffGLX9Rf//pXvetd79KOHTt04sSJco9WUc6cOaPrrrtO3/nOd8o9SsUbGRnR7t279cQTT+jo0aP63//+p23btunMmTPlHq3itLW16YEHHtCxY8d07Ngxvfe979X73/9+Pfvss+UeraI9+eSTOnDggK699tpyj1LRrrrqKr344ouL2/Hjx8s9UkU6ffq0brzxRl100UU6cuSInnvuOX3rW99Sc3Nz6Ycxq8A73vEOs2vXriWPXXHFFeYLX/hCmSaqfJLM4cOHyz1G1Th16pSRZEZGRso9SlV43eteZ77//e+Xe4yKlU6nTSwWM0ePHjX/93//Z/bs2VPukSrSl7/8ZXPdddeVe4yq8PnPf97cdNNN5R7DGGNM1Z+x+O9//6unnnpK27ZtW/L4tm3b9Kc//alMU2G1SaVSkqSWlpYyT1LZzp49q0OHDunMmTPaunVrucepWLt371ZXV5duueWWco9S8cbGxtTa2qqOjg595CMf0cTERLlHqkiPPfaYtmzZog9+8INav369Nm3apIcffrgss1R9sXj55Zd19uxZXXrppUsev/TSS/XSSy+VaSqsJsYY7d27VzfddJOuvvrqco9TkY4fP67GxkbV1dVp165dOnz4sK688spyj1WRDh06pKefflp9fX3lHqXivfOd79TAwIASiYQefvhhvfTSS7rhhhs0PT1d7tEqzsTEhPbv369YLKZEIqFdu3apt7dXAwMDJZ+l5N9uWiyv/Qp2Y4zTr2XH2nXvvffqmWee0R/+8Idyj1KxNm7cqNHRUSWTSQ0ODmrnzp0aGRmhXLzG5OSk9uzZo1/96le65JJLyj1OxduxY8fi/77mmmu0detWRaNRHTx4UHv37i3jZJUnl8tpy5Ytuv/++yVJmzZt0rPPPqv9+/frzjvvLOksVX/G4vWvf71qamrOOTtx6tSpc85iAH719PToscce0+OPP662trZyj1OxLr74Yr35zW/Wli1b1NfXp+uuu07xeLzcY1Wcp556SqdOndLmzZtVW1ur2tpajYyM6MEHH1Rtba3Onj1b7hErWkNDg6655hqNjY2Ve5SKs2HDhnOK/Fvf+tayfIih6ovFxRdfrM2bN+vo0aNLHj969KhuuOGGMk2FameM0b333quhoSH95je/UUdHR7lHqirGGM3Pz5d7jIpz88036/jx4xodHV3ctmzZottvv12jo6Oqqakp94gVbX5+Xn/729+0YcOGco9ScW688cZzPhL/97//XZdffnnJZ1kVb4Xs3btXd9xxh7Zs2aKtW7fqwIEDOnHihHbt2lXu0SpKJpPRP/7xj8X//uc//6nR0VG1tLTosssuK+NklWf37t36yU9+op/97GcKBoOLZ8RCoZDq6+vLPF1lue+++7Rjxw61t7crnU7r0KFD+u1vf6vh4eFyj1ZxgsHgOdfpNDQ0KBwOc/3OeXzuc5/Trbfeqssuu0ynTp3S17/+dc3Ozmrnzp3lHq3ifPazn9UNN9yg+++/Xx/60If0l7/8RQcOHNCBAwdKP0x5P5TizkMPPWQuv/xyc/HFF5u3v/3tfCzwPB5//HEj6Zxt586d5R6t4pwvJ0nmhz/8YblHqzh333334s9eJBIxN998s/nVr35V7rGqBh83Xd6HP/xhs2HDBnPRRReZ1tZW84EPfMA8++yz5R6rYv385z83V199tamrqzNXXHGFOXDgQFnm4GvTAQCAM1V/jQUAAKgcFAsAAOAMxQIAADhDsQAAAM5QLAAAgDMUCwAA4AzFAgAAOEOxAAAAzlAsAACAMxQLAADgDMUCAAA4Q7EAAADO/D/gx+ZPjUkUZwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_Sa_ = D_Sa.project_out(isl.dim_type.SET, 0, 1)\n",
    "D_Sc_ = D_Sc.project_out(isl.dim_type.SET, 0, 1)\n",
    "D_jk = isl.union_set(D_Sa_).union(D_Sc_)\n",
    "print(D_jk)\n",
    "plot_set_points(D_jk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ci",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}