SBM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import absolute_import\n",
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "\n",
    "import edward as ed\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "from edward.models import Bernoulli, Multinomial, Beta, Dirichlet, PointMass, Normal\n",
    "from sklearn.metrics.cluster import adjusted_rand_score\n",
    "import matplotlib.pyplot as plt\n",
    "ed.set_seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "netw = pd.read_csv('/Users/oliver/Dropbox/EDU/DTU/Blockchain-Transaction-Classification/Data/network_subsample_2.csv', sep=\",\", header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Editing data:\n",
    "netw_uniq1 = netw[netw['userID_send'] != netw['userID_recv']]                                            #Removing all users who sent to themselves\n",
    "netw_uniq = netw_uniq1.drop_duplicates(subset=['userID_send','userID_recv'], keep='last', inplace=False)  #Removing identical pairs, keeping the last ones\n",
    "\n",
    "#Finding x% last made links:\n",
    "netw_uniq=netw_uniq.sort_values(by='unixtime', axis=0, ascending=True)\n",
    "unix = int(np.round(len(netw_uniq)*0.99))\n",
    "sort = netw_uniq.iloc[unix][1]\n",
    "\n",
    "#Creating subset:\n",
    "data = netw_uniq                                                                                         #Creating a subset\n",
    "\n",
    "#Editing matrix:\n",
    "data.drop(['tx_id'], axis=1)                                                                            #Dropping column w. tx_id\n",
    "cols = ['userID_send','userID_recv','unixtime']                                                          #Rearraning columns\n",
    "data = data[cols]                                                                                         #Implementing rearranging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nodes = data.iloc[:, 0].tolist() + data.iloc[:, 1].tolist()\n",
    "nodes = sorted(list(set(nodes)))\n",
    "nodes = [(i,nodes[i]) for i in range(len(nodes))]\n",
    "for i in range(len(nodes)):\n",
    "    data = data.replace(nodes[i][1], nodes[i][0])\n",
    "\n",
    "from scipy.sparse import coo_matrix\n",
    "M = coo_matrix((data.iloc[:,2], (data.iloc[:,0],data.iloc[:,1])), shape=(len(nodes), len(nodes)))\n",
    "M_zeroing = M.todense()\n",
    "M_fullData = M.todense()\n",
    "M_originalZero= M.todense()\n",
    "M_fullData[M_fullData>0]=1                           #Making all non-zeros into ones\n",
    "OnesBeforeZeroing=(M_fullData>0).sum()               #Number of 1'nes before sort= 125944\n",
    "\n",
    "#Before touching:\n",
    "(M_fullData>sort).sum()  #=1258\n",
    "(M_fullData==0).sum()    #=99255017\n",
    "(M_fullData>0).sum()     #=125944\n",
    "\n",
    "#Adjecency matrix with only 99% of links:\n",
    "M_zeroing[M_zeroing>sort]= 0                          #Making all the last 1% into zeros\n",
    "M_zeroing[M_zeroing>0]=1                              #Making all non-zeros into ones\n",
    "OnesAfterZeroing=(M_zeroing>0).sum()                  #Number of 1'nes after sort =124686\n",
    "\n",
    "#Before running 61-63:\n",
    "(M_zeroing>sort).sum()  #=1258\n",
    "(M_zeroing==0).sum()    #=99255017\n",
    "(M_zeroing>0).sum()     #=125944\n",
    "#After:\n",
    "(M_zeroing>sort).sum()  #=0\n",
    "(M_zeroing==0).sum()    #=99256275 (=99255017+1258)\n",
    "(M_zeroing>0).sum()     #=124686 (125944-1258)\n",
    "\n",
    "percentage = OnesAfterZeroing/OnesBeforeZeroing       #Checking the number fits with 1%\n",
    "\n",
    "#Creating dataset consisting of only ones we have removed:\n",
    "M_onesRemoved= M_fullData-M_zeroing\n",
    "(M_onesRemoved>0).sum() #=1258\n",
    "\n",
    "#Creating a dataset consisting of only the correct zero's - same amount as correct ones (=1258)\n",
    "(M_originalZero==0).sum()\n",
    "(M_originalZero==1).sum()\n",
    "M_originalZero[M_originalZero>0]=1\n",
    "\n",
    "#Defining variables and model:\n",
    "x_train= M_zeroing\n",
    "M_fullDataA= M_fullData\n",
    "M_onesRemovedA=M_onesRemoved"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N = x_train.shape[0]  # number of vertices\n",
    "K = 10  # number of clusters\n",
    "gamma = Dirichlet(concentration=tf.ones([K]))\n",
    "Pi = Beta(concentration0=tf.ones([K, K]), concentration1=tf.ones([K, K]))\n",
    "Z = Multinomial(total_count=1.0, probs=gamma, sample_shape=N)\n",
    "Z1 = tf.matmul(Z, tf.matmul(Pi, tf.transpose(Z)))\n",
    "Z1 = Z1-tf.multiply(Z1,tf.diag(tf.ones(N))) + 1e-12*tf.diag(tf.ones(N))\n",
    "X = Bernoulli(probs = Z1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inference (EM algorithm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "qgamma = PointMass(tf.nn.softmax(tf.get_variable(\"qgamma/params\", [K])))\n",
    "qPi = PointMass(tf.nn.sigmoid(tf.get_variable(\"qPi/params\", [K, K])))\n",
    "qZ = PointMass(tf.nn.softmax(tf.get_variable(\"qZ/params\", [N, K])))\n",
    "\n",
    "inference = ed.MAP({gamma: qgamma, Pi: qPi, Z: qZ}, data={X: x_train})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n_iter = 200\n",
    "inference.initialize(n_iter=n_iter)\n",
    "tf.global_variables_initializer().run()\n",
    "info_loss = np.zeros(n_iter)\n",
    "for _ in range(inference.n_iter):\n",
    "  info_dict = inference.update()\n",
    "  inference.print_progress(info_dict)\n",
    "  info_loss[_] = info_dict['loss']\n",
    "inference.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of (minus) Log likelihood "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(info_loss)\n",
    "#plt.savefig('SBM.png',dpi = 400)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Criticism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sess = ed.get_session()\n",
    "x_post = ed.copy(Z1, {gamma: qgamma,\n",
    "                     Pi: qPi,\n",
    "                     Z: qZ,})\n",
    "x_gen = sess.run(x_post)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Last data processing before ROC\n",
    "#Getting pi for M_onesRemoved:\n",
    "#type(M_onesRemovedA)\n",
    "M_onesRemovedtf = tf.convert_to_tensor(M_onesRemovedA, np.float32)\n",
    "pi_onesRemoved = tf.multiply(M_onesRemovedtf,x_gen) #ændret her !!!!!\n",
    "pi_onesRemoved_matrix = pi_onesRemoved.eval()\n",
    "pi_onesRemoved_array = np.asarray(pi_onesRemoved_matrix).reshape(-1)\n",
    "pi_onesRemoved_array = pi_onesRemoved_array[pi_onesRemoved_array!=0]          #All probabilities for ones_removed\n",
    "nrOfZeros = len(pi_onesRemoved_array)\n",
    "\n",
    "#Getting pi for M_fulldata\n",
    "#type(M_fullDataA)\n",
    "M_fullDataA = np.asarray(pi_onesRemoved_matrix).reshape(-1)\n",
    "where_zero = np.where(M_fullDataA==0)[0]\n",
    "where_zero_index = np.random.choice(where_zero,nrOfZeros)\n",
    "pi_zeros = x_gen\n",
    "pi_array = np.asarray(pi_zeros).reshape(-1)\n",
    "pi_originalZeros = pi_array[where_zero_index]                                 #All probabilities for correct zeros\n",
    "\n",
    "#Creating arrays with zeros and ones:\n",
    "zeros = np.zeros(nrOfZeros)\n",
    "ones = np.ones(nrOfZeros)\n",
    "#Setting together:\n",
    "y_test = np.concatenate((zeros, ones), axis=0)\n",
    "p = np.concatenate((pi_originalZeros, pi_onesRemoved_array), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "from scipy.io import loadmat\n",
    "from sklearn import cross_validation\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn import metrics\n",
    "\n",
    "def rocplot(p, y):\n",
    "    fpr, tpr, thresholds = metrics.roc_curve(y.ravel(),p.ravel())\n",
    "    AUC = metrics.roc_auc_score(y.ravel(), p.ravel()\n",
    "    plt.plot(fpr, tpr, 'r', [0, 1], [0, 1], 'k')\n",
    "    plt.xlim([-0.01,1.01]); ylim([-0.01,1.01])\n",
    "    plt.xlabel('False positive rate (1-Specificity)')\n",
    "    plt.ylabel('True positive rate (Sensitivity)')\n",
    "    plt.title('Receiver operating characteristic (ROC)\\n AUC={:.3f}'.format(AUC))\n",
    "    plt.show()\n",
    "    return AUC#, tpr, fpr\n",
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (10, 3)\n",
    "print(rocplot(p, y_test))"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}