Wedges and 4-Cycles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <span style=\"color:blue\"> Alternative Clustering -- motifs</span> \n",
    "Nick Zolman"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:red\"> Wedges </span>\n",
    "In the following algorithm, we compute a \"wedge\" version of the clustering coffecient. Where we have a wedge motif (which can be incoming or outgoing, e.g. $a \\to c, b \\to c$ or $ c \\to a, c \\to b$). Inspiration come from the intensity of motifs discussed in Onnela et al's paper: https://journals.aps.org/pre/pdf/10.1103/PhysRevE.71.065103\n",
    "\n",
    "For a given node $a$ with weights $\\{a_i\\}$ (again, could be incoming or outgoing weights--though not mixed), we calculate the following:\n",
    "\n",
    "$$ C_{\\text{wedge, a}} = \\frac{2}{k_a(k_a - 1)}\\sum_{i<j} \\sqrt{|\\hat{a}_i \\hat{a}_j|} $$\n",
    "\n",
    "\n",
    "Where $ \\hat{a}_i = \\frac{a_i}{ \\max_{j} (a_j)}$ is a normalization and $k_a$ is $a$'s incoming or outgoing degree (though not both). \n",
    "\n",
    "We calculate this below in the following way:\n",
    "\n",
    "Let $a$ be an array of weights. $A = a a^T$. Subtract off the diagonal elements of A: $B = A - \\text{diag}(A)$. Then we flatten the array, apply the square root to each element, and sum over the elements. We have to divide by 2 to get rid of double counting the elements (we take the sum over all $(i,j)$ and not just when $i < j$. This is purely because I'm lazy, and I'm convinced NumPy is better than \"for\" loops.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n",
      "/Users/nzolman/anaconda3/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: compiletime version 3.6 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.5\n",
      "  return f(*args, **kwds)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "conv2d_3 (Conv2D)            (None, 32, 26, 26)        320       \n",
      "_________________________________________________________________\n",
      "conv2d_4 (Conv2D)            (None, 32, 24, 24)        9248      \n",
      "_________________________________________________________________\n",
      "max_pooling2d_1 (MaxPooling2 (None, 32, 12, 12)        0         \n",
      "_________________________________________________________________\n",
      "dropout_1 (Dropout)          (None, 32, 12, 12)        0         \n",
      "_________________________________________________________________\n",
      "flatten_1 (Flatten)          (None, 4608)              0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 128)               589952    \n",
      "_________________________________________________________________\n",
      "dropout_2 (Dropout)          (None, 128)               0         \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 10)                1290      \n",
      "=================================================================\n",
      "Total params: 600,810\n",
      "Trainable params: 600,810\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "from keras.models import load_model\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# load a test model\n",
    "model_filename = 'keras_mnist_modelv1.h5'\n",
    "model = load_model(model_filename)\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wedge Motifs\n",
    "\n",
    "I allow for the normalization to change in case one wanted more control over whether the normalization was globally calculated (within a layer) or locally calculated (at a node)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Given an array of weights, calculate the pairwise geometric mean and sum over it\n",
    "\n",
    "\n",
    "def wedge_cluster(node, norm = 1):\n",
    "    m = norm # maximum value\n",
    "    \n",
    "    node = node/m  # \"normalize\" the vector\n",
    "\n",
    "    nodeMat = np.outer(node, node) # Compute the outer product\n",
    "\n",
    "    nodeMat = nodeMat - np.diag(np.diag(nodeMat)) # subtract off the diagonal\n",
    "\n",
    "    sqrtArr = np.sqrt(np.abs(nodeMat.flatten())) # flatten and take the square root\n",
    "\n",
    "    return 2* np.sum(sqrtArr)/(2* node.size * (node.size - 1)) # normalize properly\n",
    "\n",
    "def wedge_local(weights):\n",
    "    return np.array([wedge_cluster(node, norm = np.max(np.abs(node))) for node in weights])\n",
    "\n",
    "def wedge_global(weights):\n",
    "    m = np.max(np.abs(weights))\n",
    "    return np.array([wedge_cluster(node, norm = m) for node in weights])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# pull in the layer and the incoming weights\n",
    "dense_layer_in= model.layers[-3]\n",
    "dense128_incoming = np.array(dense_layer_in.get_weights())[0].T\n",
    "\n",
    "dense_layer_out= model.layers[-1]\n",
    "dense128_outgoing = np.array(dense_layer_out.get_weights())[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wedge_128_in_local = wedge_local(dense128_incoming)\n",
    "wedge_128_in_global = wedge_global(dense128_incoming)\n",
    "\n",
    "wedge_128_outgoing_local = wedge_local(dense128_outgoing)\n",
    "wedge_128_outgoing_global = wedge_global(dense128_outgoing) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.title('Incoming Weights')\n",
    "plt.plot(wedge_128_in_local, label = 'local')\n",
    "plt.plot(wedge_128_in_global, label = 'global')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.title('Outgoing Weights')\n",
    "plt.plot(wedge_128_outgoing_local, label = 'local')\n",
    "plt.plot(wedge_128_outgoing_global, label = 'global')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.title('Incoming')\n",
    "plt.hist(wedge_128_in_local, 20, alpha = 0.4, label = 'local')\n",
    "plt.hist(wedge_128_in_global, 20, alpha = 0.4, label = 'global')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.title('Outgoing')\n",
    "plt.hist(wedge_128_outgoing_local, 20, alpha = 0.4, label = 'local')\n",
    "plt.hist(wedge_128_outgoing_global, 20, alpha = 0.4, label = 'global')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discussion: \n",
    "\n",
    "As you can see, there is certainly a qualitative difference for both the incoming and outgoing weights. If the norms were computed at the node level, we get larger coefficients. And as one would expect, we can get large peaks for local corresponding to similarly valued low-weights on a node."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"color:red\"> 4-cycle clustering </span>\n",
    "\n",
    "Suppose we have a left layer, $L$ with weights $L_{l,i}$ corresponding from node $l$ in the layer to node $i$ in the middle layer. And suppose we have a right layer $R$ with weights $R_{i,r}$ corresponding to the weight from node $i$ in the middle layer and node $r$ in the right layer. \n",
    "\n",
    "So for nodes $i, j$ in the middle layer, we can calculate the four-cycle clustering coefficient corresponding to it. \n",
    "\n",
    "$$ C_{\\text{4-cylce, i}} = \\frac{1}{n_L n_R (n-1)} \\sum_{\\substack{j \\neq i\\\\l \\leq n_L\\\\r \\leq n_R}} \\hat{L}_{l,i} \\hat{L}_{l,j} \\hat{R}_{i,r} \\hat{R}_{j,r}$$\n",
    "\n",
    "\n",
    "Where $\\hat{L}_{l,i} = \\sqrt[4]{\\frac{|L_{l,i}|}{\\text{max}_l(|L_{l,i}|)}}$ and likewise for $\\hat{R}_{i,r}$.\n",
    "\n",
    "We can calculate this by calculating $\\hat{L}^T \\hat{L}$ and $\\hat{R} \\hat{R}^T $\n",
    "\n",
    "$$ C_{\\text{4-cylce, a}} = \\sum_{j \\neq i} \\left(\\hat{L}^T \\hat{L})\\right)_{ij} \\left(\\hat{R} \\hat{R}^T \\right)_{ij}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def quad_coeff(in_weights, out_weights, norm_in = 0, norm_out = 0):\n",
    "\n",
    "    n_in = in_weights.shape[0]\n",
    "    n = in_weights.shape[1]\n",
    "    n_out = out_weights.shape[1]\n",
    "    motif_poss = (n_in)*(n_out)*(n-1)\n",
    "\n",
    "    # Default norm is the respective layer norms\n",
    "    if (norm_out * norm_in == 0):\n",
    "        norm_in = np.max(np.abs(in_weights))\n",
    "        norm_out = np.max(np.abs(out_weights))\n",
    "    \n",
    "    # normalize the weights\n",
    "    in_weights = in_weights/norm_in\n",
    "    out_weights = out_weights/norm_out\n",
    "    \n",
    "    # transform entries to be to the 1/4th power\n",
    "    in_weights = np.power(np.abs(in_weights), .25)\n",
    "    out_weights = np.power(np.abs(out_weights), .25)\n",
    "  \n",
    "    # Caclulate A^t.A and B.B^t\n",
    "    in_matrix = np.matmul(in_weights.T, in_weights)\n",
    "    out_matrix = np.matmul(out_weights, out_weights.T)\n",
    "    \n",
    "    # subtract off the diagonals\n",
    "    in_matrix = in_matrix - np.diag(np.diag(in_matrix))\n",
    "    out_matrix = out_matrix - np.diag(np.diag(out_matrix))\n",
    "    \n",
    "    # entry-wise multiplication\n",
    "    quad_matrix = np.multiply(in_matrix, out_matrix)\n",
    "    \n",
    "    # Return normalized coeff\n",
    "    return(1/(motif_poss) * np.sum(quad_matrix, axis = 0))\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(quad_coeff(dense128_incoming.T, dense128_outgoing))\n",
    "plt.show()\n",
    "\n",
    "plt.hist(quad_coeff(dense128_incoming.T, dense128_outgoing), 20, alpha = 0.4, label = 'global')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}