{
 "cells": [
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import cycle\n",
    "from sklearn.cluster import AffinityPropagation\n",
    "import sklearn\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import imageio\n",
    "from io import BytesIO"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_gif(figures, filename, fps=10, **kwargs):\n",
    "    images = []\n",
    "    for fig in figures:\n",
    "        output = BytesIO()\n",
    "        fig.savefig(output)\n",
    "        plt.close(fig)  \n",
    "        output.seek(0)\n",
    "        images.append(imageio.imread(output))\n",
    "    imageio.mimsave(filename, images, fps=fps, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 50\n",
    "size = (n, 2)\n",
    "np.random.seed(1)\n",
    "x = np.random.normal(0, 1, size)\n",
    "x = np.append(x, np.random.normal(5, 1, size), axis=0)\n",
    "\n",
    "c = ['r' for _ in range(n)] + ['b' for _ in range(n)]\n",
    "plt.scatter(x[:, 0], x[:, 1], c=c)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "def similarity(xi, xj):\n",
    "    return -((xi - xj)**2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def create_matrices():\n",
    "    S = np.zeros((x.shape[0], x.shape[0]))\n",
    "    R = np.array(S)\n",
    "    A = np.array(S)\n",
    "    # when looking in row i, the value means you should compare to column i - value\n",
    "    for i in range(x.shape[0]):\n",
    "        for j in range(x.shape[0]):\n",
    "            S[i, j] = similarity(x[i], x[j])\n",
    "            \n",
    "    return A, R, S\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ r(i, k) \\leftarrow s(i, k) - \\max\\limits_{k' s.t. k' \\neq k}\\{ a(i, k') + s(i, k') \\}$"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_r(lmda=0.9, slow=False):\n",
    "    global R\n",
    "    if slow:\n",
    "        for i in range(x.shape[0]):\n",
    "            for k in range(x.shape[0]):\n",
    "                v = S[i, :] + A[i, :]\n",
    "                v[k] = -np.inf\n",
    "                v[i]= -np.inf\n",
    "                R[i, k] = R[i, k] * lmda + (1 - lmda) * (S[i, k] - np.max(v))\n",
    "    \n",
    "    else:\n",
    "        # For every column k, except for the column with the maximum value the max is the same.\n",
    "        # So we can subtract the maximum for every row, and only need to do something different for k == argmax\n",
    "        \n",
    "        v = S + A\n",
    "        rows = np.arange(x.shape[0])\n",
    "        # We only compare the current point to all other points, so the diagonal can be filled with -infinity\n",
    "        np.fill_diagonal(v, -np.inf)\n",
    "\n",
    "        # max values\n",
    "        idx_max = np.argmax(v, axis=1)\n",
    "        first_max = v[rows, idx_max]\n",
    "\n",
    "        # Second max values. For every column where k is the max value.\n",
    "        v[rows, idx_max] = -np.inf\n",
    "        second_max = v[rows, np.argmax(v, axis=1)]\n",
    "\n",
    "        # Broadcast the maximum value per row over all the columns per row.\n",
    "        max_matrix = np.zeros_like(R) + first_max[:, None]\n",
    "        max_matrix[rows, idx_max] = second_max\n",
    "\n",
    "        new_val = S - max_matrix\n",
    "\n",
    "        R = R * lmda + (1 - lmda) * new_val"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "A, R, S = create_matrices()\n",
    "\n",
    "%timeit update_r(slow=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ a(i, k) \\leftarrow \\min\\{0, r(k,k) + \\sum\\limits_{i' s.t. i' \\notin \\{i, k\\}}{\\max\\{0, r(i', k)\\}}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ a(k, k) \\leftarrow \\sum\\limits_{i' \\neq k}\\max(0, r(i', k)) $"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def update_a(lmda=0.9, slow=False):\n",
    "    global A\n",
    "    \n",
    "    if slow:\n",
    "        for i in range(x.shape[0]):\n",
    "            for k in range(x.shape[0]):\n",
    "                v = np.array(R[:, k])\n",
    "                if i != k:\n",
    "                    v[i] = -np.inf\n",
    "                    v[k] = - np.inf\n",
    "                    v[v < 0] = 0\n",
    "\n",
    "                    A[i, k] = A[i, k] * lmda + (1 - lmda) * min(0, R[k, k] + v.sum())\n",
    "                else:\n",
    "                    v[k] = -np.inf\n",
    "                    v[v < 0] = 0\n",
    "                    A[k, k] = A[k, k] * lmda + (1 - lmda) * v.sum()\n",
    "    else:\n",
    "        k_k_idx = np.arange(x.shape[0])\n",
    "        # set a(i, k)\n",
    "        v = np.array(R)\n",
    "        v[v < 0] = 0\n",
    "        np.fill_diagonal(v, 0)\n",
    "        v = v.sum(axis=0) # columnwise sum\n",
    "        v = v + R[k_k_idx, k_k_idx]\n",
    "\n",
    "        # broadcasting of columns 'r(k, k) + sum(max(0, r(i', k))) to rows.\n",
    "        v = np.ones(A.shape) * v\n",
    "\n",
    "        # For every column k, subtract the positive value of k. \n",
    "        # This value is included in the sum and shouldn't be\n",
    "        v -= np.clip(R, 0, np.inf)\n",
    "        v[v > 0] = 0\n",
    "        \n",
    "        # set(a(k, k))\n",
    "        v_ = np.array(R)\n",
    "        np.fill_diagonal(v_, 0)\n",
    "\n",
    "        v_[v_ < 0] = 0\n",
    "\n",
    "        v[k_k_idx, k_k_idx] = v_.sum(axis=0) # column wise sum\n",
    "        A = A * lmda + (1 - lmda) * v\n",
    "\n",
    "def plot_iteration(A, R):\n",
    "    fig = plt.figure(figsize=(12, 6))\n",
    "    sol = A + R\n",
    "    \n",
    "    labels = np.argmax(sol, axis=1)\n",
    "\n",
    "    exemplars = np.unique(labels)\n",
    "    colors = dict(zip(exemplars, cycle('bgrcmyk')))\n",
    "    \n",
    "    for i in range(len(labels)):\n",
    "        X = x[i][0]\n",
    "        Y = x[i][1]\n",
    "        \n",
    "        if i in exemplars:\n",
    "            exemplar = i\n",
    "            edge = 'k'\n",
    "            ms = 10\n",
    "        else:\n",
    "            exemplar = labels[i]\n",
    "            ms = 3\n",
    "            edge = None\n",
    "            plt.plot([X, x[exemplar][0]], [Y, x[exemplar][1]], c=colors[exemplar])\n",
    "        plt.plot(X, Y, 'o', markersize=ms,  markeredgecolor=edge, c=colors[exemplar])\n",
    "        \n",
    "\n",
    "    plt.title('Number of exemplars: %s' % len(exemplars))\n",
    "    return fig, labels, exemplars\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.median(S)"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "A, R, S = create_matrices()\n",
    "preference = np.median(S)\n",
    "preference =  -1000\n",
    "\n",
    "np.fill_diagonal(S, preference)\n",
    "damping = 0.5\n",
    "\n",
    "figures = []\n",
    "for i in range(50):\n",
    "    update_r(damping)\n",
    "    update_a(damping, 0)\n",
    "    \n",
    "    if i % 5 == 0:\n",
    "        fig, labels, exemplars = plot_iteration(A, R)\n",
    "        figures.append(fig)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "make_gif(figures, 'test.gif', 2)"
   ]
  }
 ],
 "metadata": {
  "git": {
   "suppress_outputs": true
  },
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}