notebooks/005-tricount3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-324a96e6c478b2b2",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Triangle counting via sparse linear algebra\n",
    "\n",
    "This example considers different ways to count triangles in a graph. It illustrates some of the constructs available in Intrepydd, although for this particular problem, a pure Numpy/Scipy version is likely faster in the sequential case. (On a multicore system, an Intrepydd implementation may be faster or more scalable.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b4ca6f642fa17c3b",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Background: Counting triangles in a social network\n",
    "\n",
    "A social network may be modeled as an undirected graph, like the one shown below.\n",
    "\n",
    "![An example of an undirected graph](http://cse6040.gatech.edu/datasets/graph-triangles-example.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-73c50c3d950b4e65",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "The _nodes_ (or _vertices_) of this graph, shown as numbered circles, might represent people, and the _edges_ (or links connecting them) might represent who is friends with whom. In this case, person 0 is friends with all the \"odd birds\" of this network, persons 1, 3, and 5, but has no direct connection to persons 2 and 4.\n",
    "\n",
    "One type of analysis one might perform on such a graph is _counting triangles_, that is, the number of relationships of the form $a$ knows $b$, $b$ knows $c$, and $c$ knows $a$. In the graph shown above, there are two such triangles: (0, 1, 3) and (0, 3, 5)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A conventional graph-based approach.** A commonly used Python-based library for graph computations is [NetworkX](https://networkx.github.io/). While it turns out to be much slower than the method described in this notebook, if you are interested in seeing how to use NetworkX, see our [supplemental notebook](./tricount-networkx.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A matrix-based method for counting triangles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a4404c6fd84c6a38",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Adjacency matrix.** Let $A$ be the _adjacency matrix_ representation of the graph, defined as follows. The entries of $A$ are either 0 or 1; and $a_{i,j}$ equals 1 if and only if there is an edge connecting nodes $i$ and $j$. For instance, for the graph shown above,\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "        0 & 1 & 0 & 1 & 0 & 1 \\\\\n",
    "        1 & 0 & 0 & 1 & 0 & 0 \\\\\n",
    "        0 & 0 & 0 & 1 & 0 & 0 \\\\\n",
    "        1 & 1 & 1 & 0 & 1 & 1 \\\\\n",
    "        0 & 0 & 0 & 1 & 0 & 0 \\\\\n",
    "        1 & 0 & 0 & 1 & 0 & 0\n",
    "    \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Observe that the relationships are symmetric. For instance, 0 and 1 are mutually connected; therefore, both $a_{0,1}$ and $a_{1, 0}$ equal 1, and in general, $A = A^T$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 1, 0, 1, 0, 1],\n",
       "       [1, 0, 0, 1, 0, 0],\n",
       "       [0, 0, 0, 1, 0, 0],\n",
       "       [1, 1, 1, 0, 1, 1],\n",
       "       [0, 0, 0, 1, 0, 0],\n",
       "       [1, 0, 0, 1, 0, 0]])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy\n",
    "\n",
    "A = numpy.array([[0, 1, 0, 1, 0, 1],\n",
    "                 [1, 0, 0, 1, 0, 0],\n",
    "                 [0, 0, 0, 1, 0, 0],\n",
    "                 [1, 1, 1, 0, 1, 1],\n",
    "                 [0, 0, 0, 1, 0, 0],\n",
    "                 [1, 0, 0, 1, 0, 0]])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a3ab862d327066fb",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Matrix-based graph analysis via linear algebra.** Given such a matrix-based representation of a graph, here is one way to count triangles using linear algebra operations.\n",
    "\n",
    "First, let $A \\cdot B$ denote matrix multiplication. That is, $C = A \\cdot B$ means $c_{i,j} = \\sum_k a_{i,k} b_{k, j}$.\n",
    "\n",
    "Next, let $A \\odot B$ denote _elementwise_ multiplication. That is, $E = A \\odot B$ means $e_{i, j} = a_{i, j} b_{i, j}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-4b7f330b5d10a7fc",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Then, here is a two-step method to compute the number of triangles, which we'll denote as $t(A)$:\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "       C & = & (A \\cdot A) \\odot A \\\\\n",
    "    t(A) & = & \\frac{1}{6} \\sum_{i, j} c_{i,j}.\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step computes a \"count\" matrix $C$. Each element, $c_{i,j}$, counts the number of triangles in which both $i$ and $j$ appear. You can see this fact algebraically; the scalar formula for $c_{i,j}$ is given by\n",
    "\n",
    "$$\n",
    "    c_{i,j} = \\underbrace{\\sum_k a_{i,k} \\cdot a_{k, j}}_{A \\cdot A} \\cdot \\underbrace{a_{i,j}}_{\\odot A}\n",
    "$$\n",
    "\n",
    "Observe that the summand involves the product $a_{i,k} a_{k, j} a_{i,j}$, which is nonzero only if all of the $(i,k)$, $(k,j)$, and $(i,j)$ edges exist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the example shown above, $c_{0, 1} = c_{1,0} = 1$ since there is only one triangle that uses the edge $(0, 1)$, whereas $c_{0, 3} = c_{3, 0} = 2$ because the edge $(0, 3)$ appears in two triangles.\n",
    "\n",
    "The second step sums all the elements of $C$. However, the sum alone will overcount the number of unique triangles by a factor of six (6), hence the additional factor of $\\frac{1}{6}$. (Why?)\n",
    "\n",
    "> Instead of summing all the entries of $A$, one can exploit symmetry and consider just the upper- or lower-triangle, which is a good additional optimization you might try once you've finished the main exercises contained herein."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-779775a9d973ec2f",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "The function, **`count_triangles(A, ret_C)`**, below, implements the above formula. That is, given a symmetric Numpy array `A` representing the adjacency matrix of a graph, this function will return the number of triangles. It assumes the input matrix `A` is a (dense!) Numpy array.\n",
    "\n",
    "The `ret_C` flag is a boolean, and when `True`, the function returns **both** the triangle count **and** the matrix $C$; otherwise, it returns just the triangle count."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "count_triangles",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, array([[0, 1, 0, 2, 0, 1],\n",
      "       [1, 0, 0, 1, 0, 0],\n",
      "       [0, 0, 0, 0, 0, 0],\n",
      "       [2, 1, 0, 0, 0, 1],\n",
      "       [0, 0, 0, 0, 0, 0],\n",
      "       [1, 0, 0, 1, 0, 0]]))\n"
     ]
    }
   ],
   "source": [
    "def count_triangles(A, ret_C=False):\n",
    "    assert (type(A) is numpy.ndarray) and (A.ndim == 2) and (A.shape[0] == A.shape[1])\n",
    "    C = (A @ A) * A\n",
    "    n_tri = C.sum() // 6\n",
    "    if ret_C:\n",
    "        return n_tri, C\n",
    "    return n_tri\n",
    "\n",
    "print(count_triangles(A, ret_C=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell does some extensive testing based on random sparse graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "count_triangles_test2",
     "locked": true,
     "points": 2,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 16 16.0\n",
      "1 20 18.0\n",
      "2 30 22.0\n",
      "3 23 22.25\n",
      "4 14 20.6\n",
      "5 19 20.333333333333332\n",
      "6 22 20.571428571428573\n",
      "7 33 22.125\n",
      "8 27 22.666666666666668\n",
      "9 25 22.9\n",
      "10 32 23.727272727272727\n",
      "11 23 23.666666666666668\n",
      "12 14 22.923076923076923\n",
      "13 16 22.428571428571427\n",
      "14 24 22.533333333333335\n",
      "15 20 22.375\n",
      "16 29 22.764705882352942\n",
      "17 23 22.77777777777778\n",
      "18 24 22.842105263157894\n",
      "19 19 22.65\n",
      "20 22 22.61904761904762\n",
      "21 24 22.681818181818183\n",
      "22 17 22.434782608695652\n",
      "23 22 22.416666666666668\n",
      "24 20 22.32\n",
      "\n",
      "==> The statistics on triangle counts for random sparse graphs are consistent with expectations!\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/richie/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4HOW1+PHvUbeK1V1kyZZt2cbd4AJuYEoIEMDkUh0S4OIESCCNlEt+ueESQnJDGslNSEINJYAxkGKCwRCKG7jIDXdblovkJsmSJVmy+vn9MSNYFjVbWs2udD7Po0e7M+/MnCk7Z2bemXdEVTHGGGM6I8zrAIwxxoQ+SybGGGM6zZKJMcaYTrNkYowxptMsmRhjjOk0SybGGGM6zZKJHxEZLCInRCS8A2XniEhhd8QVbEQkW0RURCK8jqUl/utRRPqLyDIRqRSRX3sdX1frzdviqRCRG0XkTa/jCAYiMlNEdru/k6tE5HURudntd4uIrPApqyKS09b4gjqZiEi0iDwhIvvdncAGEbnUp3/zDu2Ez9+P2hjfPhG5qK1pquoBVY1X1caunJe2iMh9IlLvNx/DfPpPEpF1IlLt/p/k009E5EEROeb+/UJEpLtiD1YtrMfbgBKgr6p+x8PQTDdp6YBHVZ9T1Ys9jutREdkpIk0icks7ZaNF5EkRqRCRIyJyt0+/LBFZJSKl/gdIIvKGiExpJ5T7gT+4v5N/qOqlqvr06c5XUCcTIAIoAM4DEoEfAQtFJNuvXJK7QOJV9SenOzGPj7Jf9JmHeFXNd2OKAv4J/BVIBp4G/ul2B2cneRUwEZgAXA7c3u3RB78hwDY9jad0A7FddOTM1/RYm4CvAes7UPY+YATO9ns+8H0RucTt9wOc/cFQ4Krm5CEi1wP5qprbzriHAFtPOfrWqGpI/QEfAle7n7MBBSI6MNyzQBNwEjgBfN9n+PnAAWCZ/ziB/wS2A5VAPnC7zzjnAIU+3/8LOOiW3Qlc2MF5ug/4ayv9LnbHKT7dDgCXuJ/fB27z6TcfWNXKuOYAhcB3gCLgMPCfPv0TgWeAYmA/8N9AmNsvHPgVztF9PnCn33JKBJ5wx3kQeAAId/vlAEuBcnf4F09nubSwbt4DfgKsdJf5m0Caf1ngKaAeqHPX/UVANPBb4JD791sg2m85/RdwxN12mrt932fZXQVcBuwCSoH/18a8PAX8CVgMVLkxfA7YAFTgHDTd18K83uyu7xLghz79+7jjLAO2Ad/jk9viaHf5HMfZYVzpF8sfgdfd5bESGOAugzJgB3BmG/MyFnjLneejzfPdwWXa2rb3FPAw8Jq7LlcDw336n+EzzZ3AdX7L4tc422w5sMLtdsBdhifcv+nALcAKn2FnAGvd4dYCM3z6vUfr21cMzgHeMXcZrwX6n+K+bAVwSztlDgIX+3z/CbDA/fw6MMr9vAC4DujrblNJ7Yx3D5/cH0a78/tlt7//clIgp81xnsrMe/0H9AdqgDP8fnAH3Q31L80ru5Xh9wEXtfCDfQaIczfA5m7NO6zPAcMBwTlDqgbO8v2BuJ9H4ewQMnzGPdz9PAs43kZc97kbcynOD/+rPv2+DbzuV/5fwHfcz+XA2T79pgCVrUxnDtCAc3obibMjrAaS3f7P4JwFJbjx7wLmu/3uwNnJZAEpwLt+y+kfwCPucuwHrMFNvMALwA9xzoRjgFkdXN/30X4y2QOMdNfde8DPWyn7FPCAz7juB1a5sabjJOWf+C2nB3F+ZH18ut3rLruv4CTd593lNRZn2xzWyrw85a6rmT7LYQ4w3v0+AWfHfJVf/I+5058I1AKj3f4/B5a76yIL2MLH22IkkAf8PyAKuABnZzjKJ5YSYLIbxzvAXuAmnIOGB4B3W5mPBJxE8B132ATc7a+Dy7S1be8pnO1/Gs4BwHN8vNOMw/lt/afb7yw3/rFu/4dx1v0gN/4Z7nprXoYRPvHfgruTdJddGfAld7zz3O+pHdi+bgdeBWLdaU7GuYTaZckE50qE4pOkgGuAze7nXwJ3AUnu+h4H/A64uYPT38cn94fv0RuSibsB/ht4xKdbPM7OMwIn0bwMLDmFhde8sQ1roVuLZzs4O81v+vxAmn/AOThHXBcBkac4b2OADJ8fwmFgntvvR80/Kp/yz+EexQKNuMnV/T7CjV9amM4cnCMR3x9XEXCOO+1aYIxPv9uB99zP7wB3+PS7mI+P/Pu7w/bx6T8Pd4eEk6QeBTJPcbncR/vJ5L99+n8NeKOVsk/xyWSyB7jM5/tngX0+y6kOiGlh2TWfbSW44/dN5Otwk0EL8/IU8Ew78/tb4CG/+DN9+q8BbnA/5+Oenbrfb/PZFmfjnFGF+fR/wWebeQp4zKff14HtPt/H08rBj7teN7TSr71l2uK25xPT4z79LgN2uJ+vB5b7TesR4H9wEvFJYGIL8XxiG3C73cLHyeRLwBq/YT7A3cG3s33dipMsJ5zKNu03rfaSSZYbv+92+BmfZZoCvIhz2ezbwJk4B3kpOAc5y4C72hj/ProwmQR7nQkAIhKGc6mhDicTA6CqJ1Q1V1UbVPWo2+9iEel7ipMoaGPal/pUch3H2cjT/Mupah7wLZwdYJGILBCRjI5MXFW3qeohVW1U1fdxji6ucXufwDl19dUX50izpf59gRPqbgEtOKaqDT7fq3GSchrOUex+n377cY72wEl2BX79mg3BSfaHReS4u5wewTlCBefSkABrRGSriNzaSmyn44jP5+Z56YgMPj2vvuurWFVr/IY5ph9X6J90/x/16X+ynel/YjsTkbNF5F0RKRaRcpyzP/9tq7X5a2t9ZAAFqtrk13+Qz3f/uDs6H1k4SaMl7S3T1ra9Zq3N6xDg7OZty92+bsS5NJeGc4bUWkxt8Y+3OWbf5dRaTM8CS4AFInLIvfEl8jRiaMsJ97//77sSQFVLVfV6VZ2Is8/4Pc6BwT04Z6oXAXeIyJgujqtFQZ9M3DuTnsA5+r1aVevbKN68A23tbqbWdrAtdheRaOAVnLqC/qqahHPNu8Xxq+rzqjoLZ+NXnMskp0N9prEVmOB3h9YEPq4424pzCaTZRE6vUq0Ep15hiE+3wTiXEME5W8ry69esAOfMJE1Vk9y/vqo6FkBVj6jqV1Q1A+ds54/t3WboqsK5jNBswCnNUdsO8el5PeTzvbVtpTP8x/k8sAjIUtVE4M+0vu36a2t9HAKy3IMw3/4H6bwCnMu+LWlvmXZmmkt9tq3mG26+irPd1rQSU3vr0D9e6OByUtV6Vf2xqo7BuZpwOc5lwi6jqmU467kjv+/bcOpKt+CcWeaqah2wGefyV8AFfTLBqbQcDVyhqid9e7hHdqNEJExEUoH/w7ksU97KuI4Cw1rp15IonGuvxUCDe1tyi7cVunFc4CagGpyjuw7dXiwic0Uk2b3NdxrwDZy6C3BOPRuBb7i3CTafmb3j/n8GuFtEBrlnQt/BuWRwStwj7oXAT0UkQUSGAHfjVDLi9vuGiGSKSDLO0U/zsIdxKid/LSJ93fUxXETOc+fvWhHJdIuX4fzIG91+74nIfa2EtRE4V5xnRhJx7l7pKi8A/y0i6SKShlMX8td2hulqCUCpqta46/0LpzDsQuAH7naTiXNE2mw1TiL+vohEisgc4AqcStrO+hcwQES+5W6PCSJyttsvUMv0X8BIEfmSOz+RIjJVREa7Z19PAr8RkQwRCReR6e7vsBinkrm13/xid7xfEJEI9y6oMe702iQi54vIePeuvAqcA7Hmbfo+EXmvjWGjRCQG58AhUkRi/BK/r2dwlmmyiJyBU1f3lN/4+uHcEHOf22kvcL6INFcD5Lc3P10hqJOJu0O7HZgEHJGPn8G40S0yDHgD57RvC87R8bw2Rvm/OCvmuIh8t73pq2olzo59Ic5O8As4R5IticapFC3BOTXuh1MBiojMFpETrQwHcANOBVolzsbzoLr3e7tHF1fhHPUcx7lWe5XbHZzLSa/iHIFswbkb5pH25q0VX8fZCeXjXM99HueHCk5F8BKc67Prgb/5DXsTTvLdhrOsXgYGuv2mAqvdZbAIp85pr9svC+dumU9R1bdwrgl/iFMf0e6P/BQ8AOS6496MM08PdOH4O+JrwP0iUomz4114CsP+GOeSzF6cRP5scw9327gSuBRne/wjcJOq7uhswO5v4jM4yekIsBvnllUI0DJ1p3kxzu/kkDvd5psjAL7rTm8tTiX+gzj1RdXAT4GV7m/+HL/xHsM5o/gOzl1Z3wcuV9WSDoQ1AGcbr8C523MpHyfOVrdp15s4B5szcOoSTwLnwkcPVfqeefwPziW8/e40fqmqb/iN71fA/aravI/5X5ybLgqARdr+LcJdQlq/tG5MYLlH1C+p6nSvYzGmq4jIRpzHAo55HUt3smRijDGm04L6MpcxxpjQYMnEGGNMp1kyMcYY02lB2Xz46UhLS9Ps7GyvwzDGmJCybt26ElVN7+x4ekwyyc7OJje3W+6AM8aYHkNE/FsBOC12mcsYY0ynWTIxxhjTaZZMjDHGdJolE2OMMZ1mycQYY0ynWTIxxhjTaZZMjDHGdJolE2NOQ1PTR68zNcbQgx5aNCbQyqvrefqDfazZW8rGguOEhwlnDk5ixvBUbpqeTUxkuNchGuMZSybGtENV+deHh/nxq9s4VlXL6AF9uerMDBoalfUHyvjZ4h28su4gD10/iTEZfdsfoTE9kCUTY9rQ2KTc88qHvLSukPGDEnnqP6cyblDiJ8q8t7OI7738IVc9vJJfXTeRKydmeBStMd4JaJ2JiFwiIjtFJE9E7mmhf7SIvOj2Xy0i2W73SBF5WkQ2i8h2EenKd38b0yGNTcr3Xt7ES+sK+foFOfz9azM+lUgA5ozqx5JvncukrCS++9Im1u0v8yBaY7wVsGQiIuHAwzjvoR4DzBORMX7F5gNlqpoDPITz7maAa4FoVR0PTAZub040xnQHVeX7L3/I39Yf5DufGcl3Lh5FRHjrP5eUuCge+dJkBibGcPuzuRSWVXdjtMZ4L5BnJtOAPFXNV9U6YAEw16/MXOBp9/PLwIUiIoACcSISAfQB6oCKAMZqzCc8sWIvr6wv5FsXjeDrF47o0DDJcVE8cfMUauubuOOv62hobApwlMYEj0Amk0FAgc/3Qrdbi2VUtQEoB1JxEksVcBg4APxKVUsDGKsxH1m7r5T/fX0Hl4wdwDc7mEia5fRL4H+vHs+WgxU8v+ZAgCI0JvgEMplIC938b8xvrcw0oBHIAIYC3xGRYZ+agMhtIpIrIrnFxcWdjdcYjp2o5a7n15OV3IdfXDsB50T51Hxu/EBmDE/l12/uorSqLgBRGhN8AplMCoEsn++ZwKHWyriXtBKBUuALwBuqWq+qRcBKYIr/BFT1UVWdoqpT0tM7/aIwY7h30VbKqup5+Maz6BsTeVrjEBHuu3IsJ2ob+OWSnV0coTHBKZDJZC0wQkSGikgUcAOwyK/MIuBm9/M1wDvqPFZ8ALhAHHHAOcCOAMZqDG9vP8prHx7m6xfkMDbj03dtnYqR/RO4eXo2C9YeYNshq+4zPV/AkolbB3IXsATYDixU1a0icr+IXOkWewJIFZE84G6g+fbhh4F4YAtOUvqLqn4YqFiNOVHbwH//Ywuj+idw+3nDu2Sc37xoBHFRETyybE+XjM+YYBbQhxZVdTGw2K/bvT6fa3BuA/Yf7kRL3Y0JlF8t2cmRihr+8IWziIrommOsxD6RzJuWxZMr9/Hdi0eRlRLbJeM1JhhZQ4+m18srquTZVfv5wrTBTB6S3KXjvnXWUATnVmNjejJLJqbX+9niHcRGhnP3Z0Z2+bgHJvbhykkZvLi2gDK7s8v0YJZMTK+2YncJ7+wo4s4LckiNjw7ING47dxgn6xt5dtX+gIzfmGBgycT0Wo1NygOvbSMzuQ+3zMgO2HTOGNCX2SPSWLDmAI1N9g4U0zNZMjG91qJNB9lxpJLvX3JGwN9Fcv3ULA6V17AyrySg0zHGK9YEvek1nl/9cfMmjU3Kb/+9iwF9Y6g4Wf+JfoHwH2cNIrFPJC+tK+TckfaArel57MzE9EobDpRxrKqOi0b3J+w0mkw5VTGR4Vw1KYMlW49QXl0f8OkZ090smZhep6GpiXd2FjEoqQ+jByZ023SvnZJFXUMTizYd7LZpGtNdLJmYXmfd/jKOV9dz0ej+p9WQ4+kam9GX0QP78tK6wm6bpjHdxZKJ6VUam5Slu4rJSu7DyP7x3TptEeGayZl8WFhOXtGJbp22MYFmycT0KhsLjnO8up7zR/Xr1rOSZp8bPxCA1zcf7vZpGxNIlkxMr9GkytJdRQxMjGHUgO6rK/E1IDGGKUOSec2SielhLJmYXmPLwXJKTtRx3sh0T85Kml06fiA7jlSSX2yXukzPYcnE9AqqTl1JWnw04wZ17l0lnXXpuAEAvL7liKdxGNOVLJmYXmFFXgmHy2s4d0RatzxX0paMpD6cOTiJ1z60S12m57BkYnqFx5bvJSE6gklZSV6HAjgV8dsOV7CvpMrrUIzpEpZMTI+380gly3YVc87wVCLCg2OTv8QudZkeJqC/LBG5RER2ikieiNzTQv9oEXnR7b9aRLLd7jeKyEafvyYRmRTIWE3P9fjyfGIiwzg7O8XrUD6SmRzL2Iy+vLPjqNehGNMlApZMRCQc513ulwJjgHkiMsav2HygTFVzgIeABwFU9TlVnaSqk4AvAftUdWOgYjU9V1FlDf/ceIhrJmcSGx1c7ZpeeEY/1u0vs5dmmR4hkGcm04A8Vc1X1TpgATDXr8xc4Gn388vAhfLpezbnAS8EME7Tgz37wX7qm5q4deZQr0P5lAtG96dJYemuYq9DMabTAplMBgEFPt8L3W4tllHVBqAcSPUrcz2tJBMRuU1EckUkt7jYfpDmk07WNfLXVfu58Iz+DEvv3qZTOmLCoETS4qN4e0eR16EY02mBTCYt3X/p/5q5NsuIyNlAtapuaWkCqvqoqk5R1Snp6faOCPNJr6wvpKy6nq/MDr6zEoCwMOH8Uf1YurOI+sYmr8MxplMCmUwKgSyf75nAodbKiEgEkAiU+vS/AbvEZU5DU5Py5Iq9TMhMZNrQ4Kl493fh6H5U1DSwbn+Z16EY0ymBTCZrgREiMlREonASwyK/MouAm93P1wDvqKoCiEgYcC1OXYsxp+TtHUXkl1Tx5dnDPG06pT2zRqQTFR7GO3apy4S4gCUTtw7kLmAJsB1YqKpbReR+EbnSLfYEkCoiecDdgO/tw+cChaqaH6gYTc/12PJ8MhJjPmq6JFjFR0dw9rAU3t5utwib0BbQeyVVdTGw2K/bvT6fa3DOPloa9j3gnEDGZ3qmDwuPs2ZvKT+8bDSRQfKQYlvOG5nOA69t5+DxkwxK6uN1OMacluD/pRlzih5fvpf46Aiun5bVfuEgcN5I5+aRZXaLsAlhlkxMj3Lw+Ele23yYG6Zm0Tcm0utwOiSnXzwD+sZYMjEhzZKJ6VGeWrkXgFtmZnsbyCkQEc4dmcbKvBIa7BZhE6IsmZgeo7KmngVrCrh03AAyk2O9DueUnDsynYqaBjYVlnsdijGnxZKJ6TFeXFtAZW0DX5k9zOtQTtmsnDRErN7EhC5LJqZHaGhs4i8r9zEtO4WJQfLOklORFBvFhMwklu+2ZGJCkyUT0yO8sfUIB4+fZH6QNp3SEeeNSGNjwXHKq+u9DsWYU2bJxIQ8VeWx5XvJTo3lotH9vQ7ntM0emU6Twvt7SrwOxZhTZsnEhLzc/WVsKjjO/FlDCQ8L3qZT2jMpK4m4qHBWWjIxIciSiQl5jy3LJyk2kqsnZ3odSqdEhocxbWgK7+cd8zoUY06ZJRMT0vaVVPHW9qPcePZgYqOC602Kp2NmThr5JVUcOn7S61CMOSWh/+szp+X51Qe8DqFLLNp0kDAREmIie8Q8zcxJA2BlXgnXTgmN5mCMATszMSGsus55D8jEzKSQaTqlPaP6J5AaF8X7e+xSlwktlkxMyFqzt5T6RmWWezTfE4SFCdOHp7IyrwT31T7GhARLJiYkNTQ28UH+MaeRxMQYr8PpUrNy0iiqrCWv6ITXoRjTYZZMTEjaWHCcypoGZvegs5JmvvUmxoSKgCYTEblERHaKSJ6I3NNC/2gRedHtv1pEsn36TRCRD0Rkq4hsFpGedfhpTluTKst2l5CRGENOv3ivw+lyWSmxZKX0YaXVm5gQErBkIiLhwMPApcAYYJ6IjPErNh8oU9Uc4CHgQXfYCOCvwB2qOhaYA1gbEwaAHYcrKDlRy+yR6UH9fvfOmJWTxqo9x6xJehMyAnlmMg3IU9V8Va0DFgBz/crMBZ52P78MXCjO3uFi4ENV3QSgqsdUtTGAsZoQoaos3VVMcmwk4zISvQ4nYGYMT6OytoHNB61JehMaAplMBgEFPt8L3W4tllHVBqAcSAVGAioiS0RkvYh8v6UJiMhtIpIrIrnFxdbaam+w71g1BWUnmTUiPaSbTmnPjOGpAHaLsAkZgUwmLf3S/e91bK1MBDALuNH9/3kRufBTBVUfVdUpqjolPT29s/GaELBsVzGxUeFMHpzsdSgBlRofzRkDEqwS3oSMQCaTQsD3Ed5M4FBrZdx6kkSg1O2+VFVLVLUaWAycFcBYTQg4UlHDzqOVzBieSlREz78RcVZOGrn7y6iptyu8JvgF8he5FhghIkNFJAq4AVjkV2YRcLP7+RrgHXWe1FoCTBCRWDfJnAdsC2CsJgQs31VMZLhwztBUr0PpFjNz0qhraCJ3X5nXoRjTroAlE7cO5C6cxLAdWKiqW0XkfhG50i32BJAqInnA3cA97rBlwG9wEtJGYL2qvhaoWE3wO15dx6bC40zNTiE2unc0KTdtaAoRYWJN0puQENBfpaouxrlE5dvtXp/PNcC1rQz7V5zbg435qO5gZg98SLE1cdERnDk4ifet3sSEgJ5/4dmEvMqaetbsK2ViZhLJsVFeh9OtZgxP48OD5fYqXxP0LJmYoLd8dwkNjcr5o/p5HUq3m5mThip8kG+3CJvgZsnEBLXKmnpW7z3GpKwk0hKivQ6n203KSqJPZLi9F94EPUsmJqgt21VMY5Ny/hm976wEICoijLOHpbDC6k1MkLNkYoKWc1ZS6pyVxPe+s5JmM4enkV9cxZHyGq9DMaZVlkxM0Fq2q5gm7Z11Jb5m5DjP1djT8CaYWTIxQanio7OSZFJ78VkJwOgBfUmJi7LnTUxQs2RigtLHZyXW5pq9yteEAksmJuhU1NSzZm8pZ9pZyUdmDk/jaEUte4qrvA7FmBZZMjFB572dRTSpMsfOSj4yM6e5SXq71GWCkyUTE1RKTtSyZm8pU4ak2FmJj8EpsWQm92HFbksmJjhZMjFB5c1tR4kIC+PC0b37Di5/IsLM4Wmsyj9GY5PVm5jgY8nEBI2C0mq2HCxn1og0EmIivQ4n6MzISaWipoEt9ipfE4QsmZigoKq8vuUw8dERzO5FLQOfihnDneVitwibYGTJxASFzQfL2XesmgtH9yM6MtzrcIJSeoK9ytcEL0smxnN1DU28vuUIAxNjmJqd4nU4QW3G8DRy99mrfE3wCWgyEZFLRGSniOSJyD0t9I8WkRfd/qtFJNvtni0iJ0Vko/v350DGaby1bHcx5SfruXxCBmEiXocT1GbmpFLb0MT6/fYqXxNcApZMRCQceBi4FBgDzBORMX7F5gNlqpoDPAQ86NNvj6pOcv/uCFScxltl1XUs21XM+EGJDE2L8zqcoHf2sFTCw8RaETZBJ5Cv7Z0G5KlqPoCILADmAtt8yswF7nM/vwz8QcQOTXsLVeXVTYcQgUvHDfA6nIB6fvWBLhvXoKQ+LNp0iMzk2HbLfuHswV02XWPaEsjLXIOAAp/vhW63FsuoagNQDqS6/YaKyAYRWSois1uagIjcJiK5IpJbXFzctdGbgNt2uIIdRyq5aHR/knrZ63g7Y3h6PAfLTnKyzupNTPAIZDJp6QzD/2mr1socBgar6pnA3cDzItL3UwVVH1XVKao6JT3dmt4IJTX1jby66RADE2M+uuXVdExOv3gU2Fti7XSZ4BHIZFIIZPl8zwQOtVZGRCKARKBUVWtV9RiAqq4D9gAjAxir6WZvbTtKZU0DV00aRHiYXdk8FVkpfYgMF/KKT3gdijEfCWQyWQuMEJGhIhIF3AAs8iuzCLjZ/XwN8I6qqoikuxX4iMgwYASQH8BYTTfaW1LFqvxjnD0slayU9q/7m0+KCAsjOzWOPZZMTBAJWDJx60DuApYA24GFqrpVRO4XkSvdYk8AqSKSh3M5q/n24XOBD0VkE07F/B2qWhqoWE33qWto4pX1hSTHRfHZsf29DidkDU+Pp7iyloqT9V6HYgzQwbu5RGSmqq5sr5s/VV0MLPbrdq/P5xrg2haGewV4pSOxmdCyZNsRSqvq+PLsoURH2JPupyunXzxshT3FJzhzcLLX4RjT4TOT33ewmzGtyis6wQd7jjF9WCrD0uK9DiekDUiMITYq3C51maDR5pmJiEwHZgDpInK3T6++gB1Wmg6rrKlnYW4B/RKi+ezYnv1MSXcIE2FYWhx7iqtQVezxLOO19s5MooB4nKST4PNXgVNhbky7mlR5eV0hNfWN3DBtMFER1iRcVxjeL57yk/WUnKjzOhRj2j4zUdWlwFIReUpV93dTTKaHWbG7hN1FJ7hq0iAG9I3xOpweY0S/BAB2Ha0kPcHeSmm81dHmVKJF5FEg23cYVb0gEEGZnuNAaTVvbjvC+EGJTM22iuKulBIXRXp8NDuPVjLT3gFjPNbRZPIS8GfgccDacDAdcrKukQVrD5DYJ5LPnznIrusHwKgBCXyQf4zahka7O854qqMXrxtU9U+qukZV1zX/BTQyE9JUlb9tKKTiZD03TB1MjL3wKiBGDUigsUnZU2R3dRlvdTSZvCoiXxORgSKS0vwX0MhMSFuRV8LWQxVcPGaAPeUeQNmpcURHhLHjSKXXoZherqOXuZqbPPmeTzcFhnVtOKYnyCs6wRtbjjAuoy+zR9i1/EAKDxNG9Itn19FKu0XYeKpDyURVhwY6ENMzlFXVsWDtAdITorl6cqbt3LrBqAF92XKogsPlNWQk9fE3luoqAAAgAElEQVQ6HNNLdbQ5lZta6q6qz3RtOCaU1Tc28dzq/TSp8sVzhliFcDcZ2d9pTWDn0UpLJsYzHb3MNdXncwxwIbAesGRiAKfC/R8bDnK4vIYvTR9CWrw999BdEmIiyUzuw/bDFZw/qp/X4ZheqqOXub7u+11EEoFnAxKRCUkf5B9jQ8FxLhrdjzMGfOo9ZibAxg7sy5JtRzleXWdvrTSeON12Lapx3jFiDHtLqli8+TCjByQwx46MPTEmIxFwXoVsjBc6WmfyKh+/cjccGA0sDFRQJnSUn6zn+TUHSImL4topWYRZhbsn0hOi6ZcQzbZDFfYaZOOJjtaZ/MrncwOwX1ULAxCPCSENjU08v3o/9Y1NfGXWUHsw0WNjMvqydGcxVbUNxEV39KdtTNfo0GUut8HHHTgtBicDHWqmVEQuEZGdIpInIve00D9aRF50+68WkWy//oNF5ISIfLcj0zPda/GWIxSUneSaszLpZw04em5sRiIK7Dhil7pM9+tQMhGR64A1OG9FvA5YLSJtNkHvvsP9YeBSYAwwT0TG+BWbD5Spag7wEPCgX/+HgNc7EqPpXh8WHmdV/jFm5aQxblCi1+EYICMxhqTYSLYesmRiul9Hz4V/CExV1SIAEUkH/o3zfvbWTAPyVDXfHWYBMBfY5lNmLnCf+/ll4A8iIqqqInIVkA9UdTBG001KKmv5+4aDDE6JtRddBRERYezAvqzeW0pNfaNddjTdqqN3c4U1JxLXsQ4MOwgo8Ple6HZrsYyqNgDlQKqIxAH/Bfy4rQmIyG0ikisiucXFxe3Phem0+sYmnl9zgPAw4YapWYSHWYV7MBk/KJGGJmWbnZ2YbtbRZPKGiCwRkVtE5BbgNWBxO8O0tJfRDpb5MfCQqrbZFKqqPqqqU1R1Snp6ejvhmK7w6qZDHK2o4bopWfY8QxDKSoklOTaSTYXHvQ7F9DLtvQM+B+ivqt8Tkf8AZuEkgA+A59oZdyGQ5fM9EzjUSplCEYkAEoFS4GzgGhH5BZAENIlIjar+oWOzZQJh/f4ycveXcf6odEb2T/A6HNMCEWFiVhJLdxZTWVPvdTimF2nvzOS3QCWAqv5NVe9W1W/jnJX8tp1h1wIjRGSoiEQBNwCL/Mos4uMWia8B3lHHbFXNVtVsdzo/s0TirZITtfxz00GGpcVx4ej+Xodj2jApMwkFNh8s9zoU04u0VwGfraof+ndU1Vz/23hbKNMgIncBS3AedHxSVbeKyP1ArqouAp4AnhWRPJwzkhtOYx5MgDU2KQtzC4gIC7MHE0NAv74xDEyMYVOBXeoy3ae9ZNLWwwPtNk+qqovxq1tR1Xt9Ptfg3G7c1jjua286JrDe21lEYdlJbpiaRWKfSK/DMR0wMTOJN7YeYV9JFdlpcZ7E8PzqA55MF+ALZw/2bNq9VXuXudaKyFf8O4rIfMBe29sLFJZV8+7OIiZlJTEhM8nrcEwHTcxKQoC/rbeGKkz3aO/M5FvA30XkRj5OHlOAKODzgQzMeK+hqYlX1hcSHx3BFRMyvA7HnILEPpGM6B/Pi7kFfOPCEUSEn26brsZ0TJtbmKoeVdUZOLfq7nP/fqyq01X1SODDM15auquYoxW1zJ00iD5R9gBcqJmWncLRilre3WnPYJnA6+j7TN4F3g1wLCaIHKmo4b0dxUzITGT0QHs/SSgaNaAv/RKieX71fj4zxu7AM4Fl577mU5rctyZGR4ZxuV3eClnhYcL1U7N4b1cxB4+f9Doc08NZMjGfsn5/GQdKq7ls3EDirSnzkHb9VOe54RfXeHdnlekdLJmYT6iubeCNrUcYkhrLmYPt7q1Ql5kcy/mj+vH8mgJq6hu9Dsf0YJZMzCcs2XaUmvpG5k4chNjDiT3Cl2cNpeSE09KzMYFiycR85ODxk+TuK2X6sFQGJNrLrnqK6cNTGT8okUeX5dPY5N/WqjFdw5KJAUBVee3Dw8RGhVvbWz2MiHDHecPZW1LFm1vtjn4TGJZMDADbDlew71gVF43pby9V6oEuGTeA7NRY/rx0D6p2dmK6niUTQ0NjE69vOUK/hGimDEnxOhwTAOFhwlfOHcamwnKW7rKHGE3Xs2RiWLW3lNKqOi4bP9DenNiDXTM5k6yUPvz89R1Wd2K6nCWTXq6mvpH3dhaRkx5vL7zq4aIjwvn+Z89gx5FKXrEGIE0Xs2TSy63IK6G6rpGLx1qle29w+YSBTMpK4tdv7uRknT13YrqOJZNe7ERtAyt2lzAuoy+ZybFeh2O6gYjww8+N5mhFLY8uy/c6HNODBDSZiMglIrJTRPJE5J4W+keLyItu/9XNb28UkWkistH92yQi1tx9ALy3s4j6xiYuskYAe5Wp2SlcPmEgD7+bx84jlV6HY3qIgCUTEQkHHgYuBcYA80RkjF+x+UCZquYADwEPut23AFNUdRJwCfCIiFgjUV2o/GQ9q/eWctaQZPol2AOKvc2PrxxLQkwE331pE/WNTV6HY3qAQJ6ZTAPyVDVfVeuABcBcvzJzgafdzy8DF4qIqGq1qja43WMAu/Wkiy3fXYyqcv6ofl6HYjyQGh/NA1eNY/PBcv783h6vwzE9QCCTySCgwOd7odutxTJu8igHUgFE5GwR2QpsBu7wSS4fEZHbRCRXRHKLi+3e+Y4qqqxhzd5SJmUlkRIX5XU4xiOXjh/IFRMz+N3bu/lgzzGvwzEhLpDJpKUHFvzPMFoto6qrVXUsMBX4gYh86lqMqj6qqlNUdUp6enqnA+4tHl++l8YmZc5IOyvp7R64ahxDUmP56nPr2H+syutwTAgLZDIpBLJ8vmcCh1or49aJJAKlvgVUdTtQBYwLWKS9SGlVHX9dtZ8JmYmkJUR7HY7xWGKfSJ64eSoAtz61lvKT9R5HZEJVIJPJWmCEiAwVkSjgBmCRX5lFwM3u52uAd1RV3WEiAERkCDAK5/3zppOeWJHPyfpG5lhdiXFlp8Xxpxsnc6C0mhsfX8WxE7Veh2RCUMCSiVvHcRewBNgOLFTVrSJyv4hc6RZ7AkgVkTzgbqD59uFZwCYR2Qj8HfiaqpYEKtbe4nh1HU+/v5/Lxg2kf1+7g8t8bPrwVB790hR2Hz3BdY98wCF7za85RQF9zkRVF6vqSFUdrqo/dbvdq6qL3M81qnqtquao6jRVzXe7P6uqY1V1kqqepar/CGScvcVfVu7jRG0Dd12Q43UoJgidf0Y/np1/NkUVtcx9eCXv7SzyOiQTQuwJ+F6ioqaev6zcy8Vj+jN6YF+vwzFBatrQFF7+6gySYyO55S9r+Z9/bqGixupRTPssmfQSz36wn4qaBr5+wQivQzFBbtSABBbdNYtbZw7l6Q/2M/vBd3n43TwqLamYNlgy6QVq6ht5csVe5oxKZ3xmotfhmBAQExnOvVeM4V9fn8WUIcn8cslOpv7033xzwQbe2XGUqtpPPfZlejlroqQXeGV9Iceq6rj93OFeh2JCzLhBiTxxy1Q2FRxnYW4Br246xD83HiIiTBifmcjogX0ZlhbH8PR4hqXHkZkca+/E6aUsmfRwjU3K48v3Mn5QIucMs7comtMzMSuJiVlJ3HvFGNbsLWVV/jHW7C3ltQ8Pf+LZlKjwMDKSYshKiaWmvonk2EiS46JIjo0iJS6KuKhwRCzZ9ESWTHq4f28/yt6SKn4/70z7EZtOi44IZ/aIdGaPcFqcUFVKq+rIL6kiv/gE+SVVFJaepLCsmt1FJ6j2e2dKXHQEg5JiyEqOZXh6PJkpfYgIs6vtPYElkx7u0WX5ZCb34dJxA7wOxXjg+dUHunV6Q1LiGJIS99H32oZGjlfXU1ZVx7GqOg6X13Do+EneOVrE2zuKiIoIY/SABCZkJjGif7wllhBmyaQHW7e/lHX7y7jvijFEhNuP1HS/6Ihw+vcN/9RDstV1DeQXV7HraCVbD1WwqbCchOgIpg1LYVp2CgkxkR5FbE6XJZMe7JGl+ST2ieS6qVntFzamG8VGRTBuUCLjBiUyd5Kyu6iSVfnHeHt7Ect2FTN9WBrnjUynT1S416GaDrJk0kPlF5/gre1HuXNODrFRtppN8AoPE84Y0JczBvSluLKWd3cWsXx3MWv2HePiMQOYNjSFMKvvC3p27aOHemz5XiLDw7h5RrbXoRjTYekJ0Vw3JYu7LshhUFIfFm06xGPL8ymutMYng50lkx6o5EQtr6wv5OqzBpFuzcybEDQwsQ+3zhzK1WdlUlRRyx/e3c2GA2Veh2XaYMmkB3rm/X3UNzbx5dnDvA7FmNMmIkweksw3LxpBZnIsL60r5G/rC2mwd9YHJUsmPUx1XQPPrNrPRaP7Mzw93utwjOm0vjGR3DpzKOeNTCd3fxlPrtxHdZ015xJsLJn0MC/lFnK8up7bzrWzEtNzhIcJnx07gOumZFFQVs0jS/Mpq6rzOizjw5JJD9LQ2MTjK/I5c3ASU4Ykex2OMV1uUlYSt84cyonaBh5dnm9vhQwiAU0mInKJiOwUkTwRuaeF/tEi8qLbf7WIZLvdPyMi60Rks/v/gkDG2VMs2XqUgtKT3H7uMGs6xfRYQ9Pi+PLsodQ3NvHY8nxKLKEEhYAlExEJBx4GLgXGAPNEZIxfsflAmarmAA8BD7rdS4ArVHU8zjvinw1UnD2FqvLosj1kp8bymTHWdIrp2QYm9uHLs4a5DZnaJa9gEMgzk2lAnqrmq2odsACY61dmLvC0+/ll4EIREVXdoKqH3O5bgRgRsXtc27B6bymbCsuZP3uYNQFueoUBiTHMnzWMusYm/vL+Xk7YO1Y8FchkMggo8Ple6HZrsYyqNgDlQKpfmauBDar6qXNZEblNRHJFJLe4uLjLAg9Fjy7LJyUuimsnZ3odijHdZkBiDDdPz+Z4dT1Pv7+P2vrG9gcyARHIZNLS4bGeShkRGYtz6ev2liagqo+q6hRVnZKenn7agYa63UcreWdHETdNH0JMpLVlZHqXIalxfGHaYA6Xn+TF3AKa1H83Y7pDIJNJIeDbwmAmcKi1MiISASQCpe73TODvwE2quieAcYa8x5bnExMZxk3Ts70OxRhPnDGwL5dPyGDHkUqWbDnidTi9UiCTyVpghIgMFZEo4AZgkV+ZRTgV7ADXAO+oqopIEvAa8ANVXRnAGENeUUUN/9hwiGsnZ5ESF+V1OMZ45pxhqZwzLIXleSUsXFvQ/gCmSwUsmbh1IHcBS4DtwEJV3Soi94vIlW6xJ4BUEckD7gaabx++C8gBfiQiG92/foGKNZQ9sWIvDU1NzJ811OtQjPHc58ZnkJMez3//cwtbDpZ7HU6vItpDri9OmTJFc3NzvQ6jW5VX1zPj529zwej+/H7emac0bHe/gc+Y7lJV28CTK/cSES786+uzSexjL9pqi4isU9UpnR2PPQEfwp7+YB9VdY18bc5wr0MxJmjERUfwhy+cxeHjNXzvpU30lAPmYGfJJERV1zXwl5V7ueCMfowe2NfrcIwJKpOHJHPPpWfw5rajPL58r9fh9AqWTELUC2sKKKuu587z7azEmJbMnzWUS8YO4Odv7CB3X6nX4fR4lkxCUF1DE48ty2fa0BQmD0nxOhxjgpKI8ItrJ5CZ3Ie7nt9gjUIGmCWTEPT3DYUcqajhzvNzvA7FmKDWNyaSP954FqXVddy90OpPAsmSSYhpbFL+vDSfsRl9OXdEmtfhGBP0xmYk8qPPjWbprmKeen+f1+H0WJZMQswbW46wt6SKO8/PsWbmjemgL54zhItG9+N/F+9g++EKr8PpkSyZhBBV5Y/v5TEsLY7PjrVm5o3pKBHhwasnkBgbyTde2ECNNQjZ5SyZhJC3txex9VAFd8wZbs3MG3OKUuOj+c11E9lddIKfvrbd63B6HEsmIaKpSfnNW7sYkhrL58/0b8nfGNMRs0ek85XZQ3l21X7e2nbU63B6FEsmIeLNbUfYdriCb144gshwW23GnK7vfnYUYwb25fsvb+JoRY3X4fQYtlcKAU1NykNv7WZYehxXTszwOhxjQlp0RDj/N+9MTtY38p2Fm2hqstuFu4IlkxDw2ubD7DxayTcvHEGEnZUY02k5/eK59/KxrMgr4YkV1txKV7A9U5Cra2jiV2/uZFT/BC6fYGclxnSVedOy+OzY/vxiyQ5rrr4LWDIJci+sOcD+Y9Xcc9kZdgeXMV1IRPj5f0wgNS6ab7ywgeq6Bq9DCmmWTIJYZU09v3t7N9OHpTJnZO99x70xgZIcF8VvrpvI3mNV/ORf27wOJ6QFNJmIyCUislNE8kTknhb6R4vIi27/1SKS7XZPFZF3ReSEiPwhkDEGs0eX5VNaVccPLjvDnnY3JkBm5KRxx3nDeWFNAa9vPux1OCErYMlERMKBh4FLgTHAPBEZ41dsPlCmqjnAQ8CDbvca4EfAdwMVX7ArLKvm0WX5XDExgwmZSV6HY0yP9u2LRjIhM5F7/raZw+UnvQ4nJAXyzGQakKeq+apaBywA5vqVmQs87X5+GbhQRERVq1R1BU5S6ZV+tng7InDPpWd4HYoxPV5URBi/u+FM6hub+PaLG2m024VPWSCTySCgwOd7odutxTKq2gCUA6kBjCkkrMwrYfHmI9w5J4dBSX28DseYXmFoWhz3XTmWVfml/HnpHq/DCTmBTCYtXeT3T/cdKdP6BERuE5FcEcktLi4+peCCVX1jE/ct2srglFi+cu4wr8Mxple5dnImn5swkIfe2sXGguNehxNSAplMCoEsn++ZwKHWyohIBJAIdPj9mqr6qKpOUdUp6ek9426nx5fvZXfRCX50+RhiIsO9DseYXkVE+NlV4+nfN4Y7n1tPaVWd1yGFjEAmk7XACBEZKiJRwA3AIr8yi4Cb3c/XAO9oL34VWn7xCR769y4uHTeAz4zp73U4xvRKibHO2xmLK2v55oINVn/SQQFLJm4dyF3AEmA7sFBVt4rI/SJypVvsCSBVRPKAu4GPbh8WkX3Ab4BbRKSwhTvBepSmJuW/XvmQmIgwfjx3rNfhGNOrTcxK4r4rx7J8dwm/+/cur8MJCRGBHLmqLgYW+3W71+dzDXBtK8NmBzK2YPPX1ftZu6+MX14zgX4JMV6HY0yvN29aFhsOlPF/7+QxemBfLh0/0OuQgpo9AR8Edh2t5KevbefckelcMznT63CMMTj1Jz+5ahxnDU7i2ws3Wvtd7bBk4rGa+ka+8cIGEmIi+NW1E+xJd2OCSExkOI98aQopsVF8+elce/9JGyyZeOxni7ez40glv7x2ol3eMiYIpSdE89jNU6ioqeeWv6yloqbe65CCkiUTD/1tfSHPfLCf+bOGcv6ofl6HY4xpxdiMRP78xcnsPlrJbc/kUtvQ6HVIQceSiUc2HCjjnr9t5pxhKdZkijEh4NyR6fzq2omsyi/lGy9soL6xyeuQgoolEw8cLj/J7c+uo3/faP5442R7p7sxIeKqMwfxP1eMYcnWo3xzgSUUXwG9Ndh8WmlVHV96Yg3VdY08O/9sUuKivA7JGHMK/nPmUBqblAde246wkYeun0RUhB0QWjLpRhU19dz05GoKSqt5+tZpjBqQ4HVIxpjT8OXZTrt5D7y2nYqaev70xcnER/fu3aml025yvLqOW55cw84jlfz5S5M5Z1ivbxzZmJD25dnD+MU1E3h/zzHmPbqK4spar0PylCWTbnDo+Emu/fMHbDlYwe/nnWV3bhnTQ1w3JYvHbprM7qJKrvj9CjYcKPM6JM9YMgmwzYXlXP2n9zlSXsPTt07jknEDvA7JGNOFLjijP698dQaREcL1j6zi2VX76Y3t1VoyCRBV5dlV+7n6T+8jwIu3T2f6cLu0ZUxPNDYjkVfvmsX04an86B9buPWptb3uaXlLJgFQVFHDV/+6nh/9YwszclJ57RuzGZPR1+uwjDEBlBQbxV9umcp9V4zhg/xjXPzQMp5bvb/XNGHfu28/6GL1jU0sWHOAX7yxk9rGJu659Axumz2MsDBrb8uY3iAsTLhl5lDOHZnOPX/bzA//voW/rjrADy8bzcyc1B7d9p4lky7Q2KS8uukQv/33LvYdq2bG8FR++vnxDE2L8zo0Y4wHhqXH8+Jt57B48xF+tng7X3xiNWcNTuLO83M4f1S/HnmAacmkE4oqa3gpt5DnVx/g4PGTjB7Yl8dvmsKFo/v16CMQY0z7RITPTRjIhaP78fK6Qv703h7mP51LZnIf5k0bzNxJGWQmx3odZpexZHIKVJX8kipW7C7h9S2HWbO3lCaFmTmp/OjyMVw8pn+PPOIwxpy+mMhwvnjOEK6fmsUbW47wwpoD/HLJTn65ZCfjByVy8Zj+zMhJY0JmYkg3rRTQZCIilwC/A8KBx1X15379o4FngMnAMeB6Vd3n9vsBMB9oBL6hqksCGau/2oZGCkpPsq+kih1HKth6qIJ1+8soch9MyukXz13n53DlpEHk9IvvztCMMSEoMjyMKyZmcMXEDPYfq+L1LUd4fcsRfv3WLn791i5io8IZl5HIuEGJjOwfT1ZKLFnJsQxMigmJJBOwZCIi4cDDwGeAQmCtiCxS1W0+xeYDZaqaIyI3AA8C17vve78BGAtkAP8WkZGq2uXtPu8rqeIP7+ZxvLqO0qo6jlfXU1pdR/nJenxvFR+SGsvZw1KZPiyV6cNTrT7EGHPahqTGccd5w7njvOGUVtWxKv8Yq/OP8eHBcp5bvZ/aho8bkAwPEwYmxpAcG0V8dATxMREkREcQGx1ORFgYIjAuI5GrPX5LayDPTKYBeaqaDyAiC4C5gG8ymQvc535+GfiDOJUNc4EFqloL7BWRPHd8H3R1kDUNjazMKyE5NorkuEgykvqQHBtFanwUQ1JjGZwSx8j+8STERHb1pI0xhpS4KC4bP5DL3HfMNzQ2cbi8hoLSagrKqikoPUlhWTXlJ+s5UdtAQWk1J2obqKptoEmhSZXKmoYenUwGAQU+3wuBs1sro6oNIlIOpLrdV/kNO8h/AiJyG3Cb+/WEiOzsmtC7TBpQ4nUQnRDq8UPoz0Ooxw8ezMONXTu6oF8HW4Bft12krXkY0hUxBDKZtFQT7f/0TmtlOjIsqvoo8Oiph9Y9RCRXVad4HcfpCvX4IfTnIdTjh9Cfh1CPH7pnHgJZq1MIZPl8zwQOtVZGRCKARKC0g8MaY4wJEoFMJmuBESIyVESicCrUF/mVWQTc7H6+BnhHnRbSFgE3iEi0iAwFRgBrAhirMcaYTgjYZS63DuQuYAnOrcFPqupWEbkfyFXVRcATwLNuBXspTsLBLbcQp7K+AbgzEHdydYOgvQTXQaEeP4T+PIR6/BD68xDq8UM3zIP0xqaSjTHGdK3gfxLGGGNM0LNkYowxptMsmXQREXlSRIpEZItPtxQReUtEdrv/k72MsS2txH+fiBwUkY3u32VextgWEckSkXdFZLuIbBWRb7rdQ2kdtDYPIbEeRCRGRNaIyCY3/h+73YeKyGp3Hbzo3pATlNqYh6dEZK/POpjkdaxtEZFwEdkgIv9yvwd8HVgy6TpPAZf4dbsHeFtVRwBvu9+D1VN8On6Ah1R1kvu3uJtjOhUNwHdUdTRwDnCn2yxPKK2D1uYBQmM91AIXqOpEYBJwiYicg9NM0kPuOijDaUYpWLU2DwDf81kHG70LsUO+CWz3+R7wdWDJpIuo6jKcO9J8zQWedj8/DVzVrUGdglbiDxmqelhV17ufK3F+SIMIrXXQ2jyEBHWccL9Gun8KXIDTXBIE/zpobR5ChohkAp8DHne/C92wDiyZBFZ/VT0Mzo4C6OdxPKfjLhH50L0MFrSXiHyJSDZwJrCaEF0HfvMAIbIe3MsrG4Ei4C1gD3BcVRvcIi02jRRM/OdBVZvXwU/ddfCQ2+J5sPot8H2gubXIVLphHVgyMW35EzAc53T/MO02/+M9EYkHXgG+paoVXsdzOlqYh5BZD6raqKqTcFqtmAaMbqlY90Z1avznQUTGAT8AzgCmAinAf3kYYqtE5HKgSFXX+XZuoWiXrwNLJoF1VEQGArj/izyO55So6lH3h9UEPIazcwhaIhKJsxN+TlX/5nYOqXXQ0jyE2noAUNXjwHs4dT9JbnNJEEJNI/nMwyXuJUh1WzL/C8G7DmYCV4rIPmABzuWt39IN68CSSWD5NhdzM/BPD2M5Zc07YdfncRonDUrudeEngO2q+hufXiGzDlqbh1BZDyKSLiJJ7uc+wEU49T7v4jSXBMG/Dlqahx0+BySCU98QlOtAVX+gqpmqmo3Tosg7qnoj3bAO7An4LiIiLwBzcJp6Pgr8D/APYCEwGDgAXKuqQVnJ3Ur8c3AurSiwD7i9uf4h2IjILGA5sJmPrxX/P5w6h1BZB63NwzxCYD2IyAScyt1wnAPVhap6v4gMwzlKTgE2AF90j/CDThvz8A6QjnPJaCNwh09FfVASkTnAd1X18u5YB5ZMjDHGdJpd5jLGGNNplkyMMcZ0miUTY4wxnWbJxBhjTKdZMjHGGNNplkxM0BMRFZFf+3z/rojc10XjfkpErmm/ZKenc63bGvC7ft2zReQLbQyXISIvt9a/E/HcIiJ/6Orxmt7LkokJBbXAf4hImteB+BKR8FMoPh/4mqqe79c9G2gxmYhIhKoeUtWAJztjOsuSiQkFDTjvsP62fw//MwsROeH+nyMiS0VkoYjsEpGfi8iN7rsqNovIcJ/RXCQiy91yl7vDh4vIL0Vkrdu43+0+431XRJ7HebjQP5557vi3iMiDbrd7gVnAn0Xkl36D/ByY7b4j49vuGcNLIvIq8KZ75rLFHU+2G+d692+GT0zvicjLIrJDRJ5zn9RGRC5zu60Qkf8T9/0WfjGni8gr7ryuFZGZbvfz5OP3d2wQkYQOrS3TK0W0X8SYoPAw8KGI/OIUhpmI09BgKZAPPK6q08R56dTXgW+55bKB83AaU3xXRHKAm4ByVZ3qthC7UkTedMtPA8ap6l7fiYlIBs57IybjvDPiTRG5yn2C+skUuZEAAALCSURBVAKcp5Fz/WK8x+3enMRuAaYDE1S1VJzWg5sVAZ9R1RoRGQG8AExx+50JjMVpc2klMFNEcoFHgHNVda/bykFLfofzrosVIjIYWOIut+8Cd6rqSnEan6xpZXhjLJmY0KCqFSLyDPAN4GQHB1vb3OyIiOwBmpPBZsD3ctNCtxHF3SKSj9M67MXABJ+znkRgBFAHrPFPJK6pwHuqWvz/27t71iiiKIzj/xNIpWjjJ4iQwsZGCDbapLIXsRArRUVFS7+BFr4UKgRrESzsJL5UgkIMsrgRQ4IgFoK1GMEiybE4Z+FmndlZM0XYzfPrZmf2zr2wzN1z7nBu3vMxcIwoq/M/XteUfJkE7lvs8rcBTBfnFt39e973IzFBrgFfi74+Ac5XtDsLHMpgBmBfRiHvgDs5jme99kWqaDKRUXIP6BBVW3vWyXRtpnbK7UjL2kObxfEmW3/7/TWFnKjBdMXdX5Ynst7R75r+VZX63o669q8TddMOE2MuI4VyrBvE+IbtzwRw1N37J+mbZvYcOAEsmNmsu68M2absMlozkZGR/9afsnXL0W9EWgliV8XJbTR90swmch1lClglUj0XLUrCY2bTZranoZ33wHEzO5CL86eBNw3f+QUMuxaxH/iRUdQZohjhICvAVJEqO1Vz3Svgcu8gIx/M7KC7f3L3W8AHImITqaTJREbNbaKycc8j4gG+CMxQ/69+kFXioT9PVIP9Q2x5ugx0cgF8joZIPlNqN4hy312g4+5Npb6XgHUz65rZPy8Y9HkInDWzBSLFNXCsGWlcAl6Y2VsiqvlZcelV4Ei+aLAMXMjPr+WLBF0itTjf0D/ZxVQ1WGSMmdled1/LFOAD4Iu7393pfsn4UWQiMt7O5YL8ZyJNNrfD/ZExpchERERaU2QiIiKtaTIREZHWNJmIiEhrmkxERKQ1TSYiItLaX3PEO/0kzfrMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test cell: `count_triangles_test2`\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "\n",
    "def check_count_triangles_large(n, den=1e-2):\n",
    "    U_large = numpy.triu(numpy.random.rand(n, n) <= den).astype(int)\n",
    "    numpy.fill_diagonal(U_large, 0)\n",
    "    A_large = U_large + U_large.T\n",
    "    return count_triangles(A_large)\n",
    "\n",
    "n, den, k_max, mu, sd = 500, 1e-2, 25, 21, 5\n",
    "nts = numpy.zeros(k_max, dtype=int)\n",
    "for k in range(k_max):\n",
    "    nts[k] = check_count_triangles_large(n, den)\n",
    "    print(k, nts[k], numpy.sum(nts[:k+1])/(k+1))\n",
    "sns.distplot(nts)\n",
    "plt.xlabel(\"Number of triangles\")\n",
    "plt.ylabel(\"Count\")\n",
    "plt.title(\"{} trials: {} nodes, uniform random connections, {}% fill\".format(k_max, n, den*100))\n",
    "\n",
    "assert (mu-sd) <= numpy.mean(nts) <= (mu+sd), \\\n",
    "       f\"mean={numpy.mean(nts)}: The statistics on triangle counts for random graphs do not match expectations.\"\n",
    "\n",
    "print(\"\\n==> The statistics on triangle counts for random sparse graphs are consistent with expectations!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b4ec96c4fe4e412e",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## The IMDB actor network dataset\n",
    "\n",
    "Later, we will apply the matrix-based triangle counting method to a real \"social\" network, namely, the graph of actors who have appeared in the same movie. That is, let each actor be a node in the graph, with an edge between two actor nodes $i$ and $j$ if they appeared in the same movie together. If there are $n$ actors, then this graph can be represented by an $n \\times n$ adjacency matrix, where $a_{i,j}=1$ if and only if $i$ and $j$ appeared in at least one movie together, and 0 otherwise. (We will only allow 0 or 1 values in the matrix, effectively ignoring the number of movies in which the actors \"co-occurred.)\n",
    "\n",
    "The dataset in this problem uses data collected on a crawl of the [Top 250 movies on the Internet Movie Database](https://github.com/napsternxg/IMDB-Graph/tree/master/Data/tutorial/tutorial) (circa 2012). If you are interested in the raw data to construct the sparse matrix corresponding to the graph, refer to the [NetworkX-based supplemental notebook](./tricount-networks.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-38fd73d5b41c7559",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<12975x12975 sparse matrix of type '<class 'numpy.int64'>'\n",
       "\twith 1246812 stored elements in Compressed Sparse Row format>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pickle\n",
    "import scipy\n",
    "\n",
    "with open('imdb_tricounts.pickle', 'rb') as fp:\n",
    "    ntri_nx = pickle.load(fp)\n",
    "with open('imdb_tricounts_matrix.pickle', 'rb') as fp:\n",
    "    tricounts_nx = pickle.load(fp)\n",
    "with open('imdb_actors.pickle', 'rb') as fp:\n",
    "    V_actors = pickle.load(fp)\n",
    "A_casts = scipy.sparse.load_npz('imdb_matrix.npz')\n",
    "A_casts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observation: sparsity.** The matrix of actor co-occurrences is quite sparse. There are $12,975$ actors, which means as many as $\\approx 12,\\!975^2 / 2 \\approx 84$ million possible co-occurences; yet, the matrix (`A_casts`) only appears to have $1.2$ million nonzeros, which is about $600,\\!000$ observed co-occurrences, or less than 1% of the possible co-occurrences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-173b9f679f64c689",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## A sparse matrix-based approach\n",
    "\n",
    "The first matrix-based example we gave above stores the graph as a _dense_ matrix. However, as noted above, the graph itself is actually quite sparse. Thus, any method that scales like $\\mathcal{O}(n^2)$ or worse, where $n$ is the number of actors, will be much slower than a method that exploits sparsity.\n",
    "\n",
    "Luckily, Numpy/Scipy has built-in support for sparse matrices and basic operations on them. Therefore, we can use the simple implementation below to implement the counting formula, and trust that Numpy/Scipy works reasonably well if the input matrix is sparse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found ~ 24.6 million triangles\n"
     ]
    }
   ],
   "source": [
    "# Triangle counter for a *sparse* adjacency matrix\n",
    "def count_triangles_spmat__baseline(A, ret_C=False): # We'll try to speed this up using Intrepydd later\n",
    "    C = (A @ A).multiply(A)\n",
    "    n_tri = C.sum() // 6\n",
    "    if ret_C:\n",
    "        return n_tri, C\n",
    "    return n_tri\n",
    "\n",
    "ntri_spmat, C_spmat = count_triangles_spmat__baseline(A_casts, ret_C=True)\n",
    "print(f\"Found ~ {ntri_spmat*1e-6:.1f} million triangles\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Passed!\n"
     ]
    }
   ],
   "source": [
    "# Test the results by building `tricounts` as described above and comparing it to NetworkX's result.\n",
    "def spmat_to_tricounts(C):\n",
    "    tricounts_ind = C.dot(numpy.ones(len(V_actors), dtype=int)) // 2\n",
    "    tricounts = {actor: count for actor, count in zip(V_actors, tricounts_ind)}\n",
    "    return tricounts\n",
    "\n",
    "def test_spmat(tricounts_spmat, tricounts_nx):\n",
    "    for actor in tricounts_nx:\n",
    "        assert tricounts_nx[actor] == tricounts_spmat[actor], \\\n",
    "               \"For actor '{}', NetworkX saw {} triangles whereas you saw {}.\".format(actors[actor]\n",
    "                                                                                      , tricounts_nx[actor]\n",
    "                                                                                      , tricounts[actor])\n",
    "    return True\n",
    "\n",
    "tricounts = spmat_to_tricounts(C_spmat)\n",
    "test_spmat(tricounts, tricounts_nx)\n",
    "print(\"Passed!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "560 ms ± 21.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "\n",
      "==> Execution time: ~ 560.5 ms\n"
     ]
    }
   ],
   "source": [
    "# Lastly, let's see if the matrix approach yields any speedup...\n",
    "\n",
    "t_spmat = %timeit -o count_triangles_spmat__baseline(A_casts)\n",
    "print(f\"\\n==> Execution time: ~ {t_spmat.average*1e3:.1f} ms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observation.** In experiments we have done on conventional desktop machines, this sparse matrix-based implementation is 10-30x faster than a NetworkX-based implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An Intrepydd version: Domain-specific wrappers\n",
    "\n",
    "The linear algebraic algorithm has a \"straightforward\" translation to Intrepydd's domain-specific wrappers. In this example, the `calc_C()` Intrepydd kernel defined below implements the Numpy/Scipy statement, `C = (A @ A).multiply(A)`, where `A` and `C` are sparse matrices. In the current version of Intrepydd, sparse matrix operations are wrappers for the [CombBLAS](https://people.eecs.berkeley.edu/~aydin/CombBLAS/html/) library."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a moment to read this implementation, and observe the following.\n",
    "- In the current version of Intrepydd, it is not possible to pass full sparse matrix objects between Python and Intrepydd. One workaround is to \"extract\" the underlying arrays representing the sparse matrix and pass those instead. Here, we assume the input `A` and the output `C` are [compressed sparse row (CSR) matrices](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html), which are represented by an array of floating-point values, integer column indices, and an integer array of row pointers.\n",
    "- Sparse matrices in Intrepydd, at present, must store the values as floating-point. By contrast, the array `A_casts` constructed above has only integer values. That will require the caller to do an extra conversion, which might affect execution time.\n",
    "- The Intrepydd call `csr_to_spm()` translates from these CSR arrays to an Intrepydd sparse matrix object. Refer to the [Intrepydd Guide function list](https://hpcgarage.github.io/intrepyddguide/library/functions.html) for additional details.\n",
    "- Given an Intrepydd sparse matrix object `A`, the Intrepydd call `spmm(A, A)` corresponds to `A @ A` in Numpy/Scipy.\n",
    "- The Intrepydd call `spm_mul(A, B)` corresponds to `A.multiply(B)`.\n",
    "- The Intrepydd call `spm_to_csr()` translates from an Intrepydd sparse matrix object back to an array representation for the caller. These are passed-by-reference arguments, which the caller of `calc_C` must supply. (See the modified Python function, `count_triangles_spmat__pydd()`, a couple code cells below.) Furthermore, this operation might also affect overall execution time, just as `csr_to_spm()` can."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting tricount_pydd.pydd\n"
     ]
    }
   ],
   "source": [
    "%%writefile tricount_pydd.pydd\n",
    "\n",
    "def calc_C(a_vals: Array(float64), a_inds: Array(int32), a_ptrs: Array(int32), a_nc: int32,\n",
    "           c_vals: Array(float64), c_inds: Array(int32), c_ptrs: Array(int32)):\n",
    "    A = csr_to_spm(a_vals, a_inds, a_ptrs, a_nc)\n",
    "    C = spm_mul(spmm(A, A), A)\n",
    "    spm_to_csr(C, c_vals, c_inds, c_ptrs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pyddc tricount_pydd.pydd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above defined the Intrepydd kernel; here is how it may be called from a Python driver function (`count_triangles_spmat__pydd`). Pay particular attention to how one gets the array representation for a Numpy/Scipy CSR object, and how one can construct the output.\n",
    "\n",
    "> Regarding the output, this implementation exploits the fact that for `C = A.multiply(A @ A)`, `C` must have the _same_ sparsity pattern as `A`. (Can you see why?) Therefore, when the driver allocates space for `C` (`c_vals`, `c_inds`, and `c_ptrs`), it is sufficient for these to have the same sizes as those for `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found ~ 24.6 million triangles\n"
     ]
    }
   ],
   "source": [
    "import tricount_pydd\n",
    "\n",
    "def count_triangles_spmat__pydd(A, ret_C=False):\n",
    "    a_vals, a_inds, a_ptrs, a_nc = A.data.astype(float), A.indices, A.indptr, numpy.int32(A.shape[1])\n",
    "    c_vals = numpy.empty(len(a_vals), numpy.float64)\n",
    "    c_inds = numpy.empty(len(a_inds), numpy.int32)\n",
    "    c_ptrs = numpy.empty(len(a_ptrs), numpy.int32)\n",
    "    tricount_pydd.calc_C(a_vals, a_inds, a_ptrs, a_nc, c_vals, c_inds, c_ptrs)\n",
    "    n_tri = int(c_vals.sum()) // 6\n",
    "    if ret_C:\n",
    "        from scipy.sparse import csr_matrix\n",
    "        C = csr_matrix((c_vals, c_inds, c_ptrs), shape=A.shape)\n",
    "        return n_tri, C\n",
    "    return n_tri\n",
    "\n",
    "ntri_pydd, C_pydd = count_triangles_spmat__pydd(A_casts, ret_C=True)\n",
    "print(f\"Found ~ {ntri_pydd*1e-6:.1f} million triangles\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Passed!\n"
     ]
    }
   ],
   "source": [
    "# Correctness check:\n",
    "tricounts_pydd = spmat_to_tricounts(C_pydd)\n",
    "test_spmat(tricounts_pydd, tricounts_nx)\n",
    "print(\"Passed!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "958 ms ± 166 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "\n",
      "==> Execution time: ~ 957.9 ms\n"
     ]
    }
   ],
   "source": [
    "t_pydd = %timeit -o count_triangles_spmat__pydd(A_casts)\n",
    "print(f\"\\n==> Execution time: ~ {t_pydd.average*1e3:.1f} ms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observation.** The time is a little bit slower than the pure Numpy/Scipy version. Again, as noted previously, there are extra overheads in the current version of Intrepydd for input/output data transfers due to `csr_to_spm()` and `spm_to_csr()` calls."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another Intrepydd version: \"lowering\" the implementation\n",
    "\n",
    "The current version of Intrepydd is translated into C++ and then compiled. So, you can also try to write \"low-level\" implementations that look like what you might have written in C++, and hope to get C++-level performance from it. The following implementation is an example of \"lowering\" the original implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, lowering makes it possible to optimize _across_ function calls. That is, recall the linear algebraic  operation that triangle counting expresses, written here in a slightly different and more suggestive way (since $A \\odot B = B \\odot A$):\n",
    "\n",
    "$$\n",
    "    c_{i,j} = \\underbrace{a_{i,j} \\, \\cdot}_{A \\odot} \\underbrace{\\sum_k a_{i,k} \\cdot a_{k, j}}_{A \\cdot A}.\n",
    "$$\n",
    "\n",
    "In this form, it should be clear that $C$ only has a nonzero output $c_{i,j}$ when $A$ has a nonzero element $a_{i,j}$. Therefore, rather than compute the entire product $B = A \\cdot A$ first, we can compute only those $b_{i,j}$ where $a_{i,j} \\neq 0$. The computation of a single $b_{i,j}$ is a _sparse dot product_ between a sparse row $i$ of $A$ and a sparse column $j$ of $A$. And since $A$ is symmetric ($A = A^T$), sparse rows and columns are \"interchangeable.\"\n",
    "\n",
    "Of course, it is not clear if such an implementation will be faster or not because there are many ways to implement $A \\cdot A$ to save computation and data movement, but from the point of view of learning Intrepydd, it is instructive to try.\n",
    "\n",
    "> _Aside._ This formula also explains why one expects $C$ to have the same sparsity pattern as $A$, as noted previously.\n",
    "\n",
    "In addition, lowering confers one more potential benefit in this case. Recall that Intrepydd sparse matrix objects must store their values in floating-point; by contrast, if we lower and do not use the sparse interface at all, we can change the implementation to use integer types, too.\n",
    "\n",
    "> _Logical operations._ At present, logical operations (`or`, `and`, `not`) are not supported in Intrepydd. Therefore, the sparse dot product code cannot use the more compact, `has_nzs_left = (ki < ni) and (kj < nj)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting tricount_pydd_lowered.pydd\n"
     ]
    }
   ],
   "source": [
    "%%writefile tricount_pydd_lowered.pydd\n",
    "\n",
    "def sparse_dot_rows_ij(a_vals: Array(int32), a_inds: Array(int32), a_ptrs: Array(int32),\n",
    "                       i: int32, j: int32) -> int32:\n",
    "    # Assumes column indices are sorted within each row\n",
    "    ki = a_ptrs[i]   # type: int32\n",
    "    ni = a_ptrs[i+1] # type: int32\n",
    "    kj = a_ptrs[j]   # type: int32\n",
    "    nj = a_ptrs[j+1] # type: int32\n",
    "    s = 0\n",
    "    has_nzs_left = 1\n",
    "    if ki >= ni: has_nzs_left = 0\n",
    "    if kj >= nj: has_nzs_left = 0\n",
    "    while has_nzs_left:\n",
    "        if a_inds[ki] < a_inds[kj]:\n",
    "            ki += 1\n",
    "        elif a_inds[ki] > a_inds[kj]:\n",
    "            kj += 1\n",
    "        else:\n",
    "            s += a_vals[ki] * a_vals[kj]\n",
    "            ki += 1\n",
    "            kj += 1\n",
    "        if ki >= ni: has_nzs_left = 0\n",
    "        if kj >= nj: has_nzs_left = 0\n",
    "    return s\n",
    "\n",
    "def calc_C(n: int32,\n",
    "           a_vals: Array(int32), a_inds: Array(int32), a_ptrs: Array(int32),\n",
    "           c_vals: Array(int32), c_inds: Array(int32), c_ptrs: Array(int32)):\n",
    "    for i in range(n):\n",
    "        for k in range(a_ptrs[i], a_ptrs[i+1]):\n",
    "            j = a_inds[k]\n",
    "            a_ij = a_vals[k]\n",
    "            c_vals[k] = a_ij * sparse_dot_rows_ij(a_vals, a_inds, a_ptrs, i, j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pyddc tricount_pydd_lowered.pydd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found ~ 24.6 million triangles\n"
     ]
    }
   ],
   "source": [
    "import tricount_pydd_lowered\n",
    "\n",
    "def count_triangles_spmat__pydd2(A, ret_C=False):\n",
    "    assert A.shape[0] == A.shape[1]\n",
    "    n = A.shape[0]\n",
    "    a_vals, a_inds, a_ptrs = A.data.astype(numpy.int32), A.indices, A.indptr\n",
    "    c_vals = a_vals.copy()\n",
    "    c_inds = a_inds.copy()\n",
    "    c_ptrs = a_ptrs.copy()\n",
    "    tricount_pydd_lowered.calc_C(n,\n",
    "                                 a_vals, a_inds, a_ptrs,\n",
    "                                 c_vals, c_inds, c_ptrs)\n",
    "    n_tri = c_vals.sum() // 6\n",
    "    if ret_C:\n",
    "        from scipy.sparse import csr_matrix\n",
    "        C = csr_matrix((c_vals, c_inds, c_ptrs), shape=A.shape)\n",
    "        return n_tri, C\n",
    "    return n_tri\n",
    "\n",
    "ntri_pydd2, C_pydd2 = count_triangles_spmat__pydd2(A_casts, ret_C=True)\n",
    "print(f\"Found ~ {ntri_pydd2*1e-6:.1f} million triangles\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Passed!\n"
     ]
    }
   ],
   "source": [
    "# Correctness check:\n",
    "tricounts_pydd2 = spmat_to_tricounts(C_pydd2)\n",
    "test_spmat(tricounts_pydd2, tricounts_nx)\n",
    "print(\"Passed!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.3 s ± 578 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "\n",
      "==> Execution time: ~ 7303.5 ms\n"
     ]
    }
   ],
   "source": [
    "t_pydd2 = %timeit -o count_triangles_spmat__pydd2(A_casts)\n",
    "print(f\"\\n==> Execution time: ~ {t_pydd2.average*1e3:.1f} ms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observation.** This Intrepydd version is not faster than the earlier Intrepydd version or the native Numpy/Scipy versions. However, having a lowered version means you can try even more optimizations, including parallelization! (All the other versions have been purely sequential.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d3f5c8a0add3fec7",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Fin!** We hope this notebook has been a useful tutorial of how to use Intrepydd's sparse matrix operations, in the context of implementing graph algorithms that exploits the duality between a sparse matrix and a sparse graph. In particular, if you are working on a graph analysis problem, you might consider whether there is a matrix version and try to implement that using the APIs and techniques described above."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "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": 1
}