CayleyNet.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import time, shutil\n",
    "import numpy as np\n",
    "import os, collections, sklearn\n",
    "import joblib\n",
    "\n",
    "import graph, coarsening\n",
    "import scipy.sparse as sp\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Graph definition and coarsening"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Definition of some flags useful later in the code\n",
    "\n",
    "flags = tf.app.flags\n",
    "FLAGS = flags.FLAGS\n",
    "\n",
    "# Graphs.\n",
    "flags.DEFINE_integer('number_edges', 8, 'Graph: minimum number of edges per vertex.')\n",
    "flags.DEFINE_string('metric', 'euclidean', 'Graph: similarity measure (between features).')\n",
    "flags.DEFINE_bool('normalized_laplacian', True, 'Graph Laplacian: normalized.')\n",
    "flags.DEFINE_integer('coarsening_levels', 4, 'Number of coarsened graphs.')\n",
    "\n",
    "# Directories.\n",
    "flags.DEFINE_string('dir_data', 'data_mnist', 'Directory to store data.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Here we proceed at computing the original grid where the images live and the various coarsening that are applied\n",
    "#for each level\n",
    "\n",
    "def grid_graph(m):\n",
    "    z = graph.grid(m)\n",
    "    dist, idx = graph.distance_sklearn_metrics(z, k=FLAGS.number_edges, metric=FLAGS.metric) \n",
    "    #dist contains the distance of the 8 nearest neighbors for each node sorted in ascending order\n",
    "    #idx contains the indexes of the 8 nearest for each node sorted in ascending order by distance\n",
    "\n",
    "    A = graph.adjacency(dist, idx)\n",
    "    return A\n",
    "\n",
    "def coarsen(A, levels):\n",
    "    graphs, parents = coarsening.metis(A, levels) #Coarsen a graph multiple times using the METIS algorithm. \n",
    "                                                  #Everything starts with a random point and then decides how to \n",
    "                                                  #combine the points.\n",
    "                                                  #Construction is done a priori, so we have one graph\n",
    "                                                  #for all the samples!\n",
    "                    \n",
    "                                                  #graphs = list of spare adjacency matrices (it contains in position \n",
    "                                                  #          0 the original graph)\n",
    "                                                  #parents = list of numpy arrays (every array in position i contains \n",
    "                                                  #           the mapping from graph i to graph i+1, i.e. the idx of\n",
    "                                                  #           node i in the coarsed graph) \n",
    "    perms = coarsening.compute_perm(parents) #Return a list of indices to reorder the adjacency and data matrices so\n",
    "                                             #that the union of two neighbors from layer to layer forms a binary tree.\n",
    "                                             #Fake nodes are appended at the end of the current graph\n",
    "    laplacians = []\n",
    "    for i,A in enumerate(graphs):\n",
    "        M, M = A.shape\n",
    "\n",
    "        # We remove any possible self-connection.\n",
    "        A = A.tocoo()\n",
    "        A.setdiag(0)\n",
    "\n",
    "        if i < levels: #if we have to pool the graph \n",
    "            A = coarsening.perm_adjacency(A, perms[i]) #matrix A is here extended with the fakes nodes\n",
    "                                                       #in order to do an efficient pooling operation\n",
    "                                                       #in tensorflow as it was a 1D pooling\n",
    "\n",
    "        A = A.tocsr()\n",
    "        A.eliminate_zeros()\n",
    "        Mnew, Mnew = A.shape\n",
    "        print('Layer {0}: M_{0} = |V| = {1} nodes ({2} added), |E| = {3} edges'.format(i, Mnew, Mnew-M, A.nnz//2))\n",
    "\n",
    "        L = graph.laplacian(A, normalized=FLAGS.normalized_laplacian)\n",
    "        laplacians.append(L)\n",
    "    return laplacians, perms[0] if len(perms) > 0 else None\n",
    "\n",
    "t_start = time.time()\n",
    "\n",
    "np.random.seed(0)\n",
    "A = grid_graph(28)\n",
    "L, perm = coarsen(A, FLAGS.coarsening_levels)\n",
    "\n",
    "print('Execution time: {:.2f}s'.format(time.time() - t_start))\n",
    "\n",
    "graph.plot_spectrum(L)\n",
    "del A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Normalize Laplacian\n",
    "L_norm = []\n",
    "for k in range(len(L)):\n",
    "    L_norm.append(L[k] - sp.eye(L[k].shape[0]))\n",
    "graph.plot_spectrum(L_norm, ymin=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#loading of MNIST dataset\n",
    "\n",
    "from tensorflow.examples.tutorials.mnist import input_data\n",
    "mnist = input_data.read_data_sets(FLAGS.dir_data, one_hot=False)\n",
    "\n",
    "train_data = mnist.train.images.astype(np.float32)\n",
    "val_data = mnist.validation.images.astype(np.float32) #the first 5K samples of the training dataset \n",
    "                                                      #are used for validation\n",
    "test_data = mnist.test.images.astype(np.float32)\n",
    "train_labels = mnist.train.labels\n",
    "val_labels = mnist.validation.labels\n",
    "test_labels = mnist.test.labels\n",
    "\n",
    "t_start = time.time()\n",
    "train_data = coarsening.perm_data(train_data, perm)\n",
    "val_data = coarsening.perm_data(val_data, perm)\n",
    "test_data = coarsening.perm_data(test_data, perm)\n",
    "print('Execution time: {:.2f}s'.format(time.time() - t_start))\n",
    "del perm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Model definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "class CayleyNet:\n",
    "    \"\"\"\n",
    "    The neural network model.\n",
    "    \"\"\"\n",
    "    \n",
    "    #Helper functions used for constructing the model\n",
    "    def _weight_variable(self, shape, regularization=True, name=\"\"): \n",
    "        \"\"\"Initializer for the weights\"\"\"\n",
    "        \n",
    "        initial = tf.truncated_normal_initializer(0, 0.1)\n",
    "        var = tf.get_variable('weights'+name, shape, tf.float32, initializer=initial)\n",
    "        if regularization: #append the loss of the current variable to the regularization term \n",
    "            self.regularizers.append(tf.nn.l2_loss(var))\n",
    "        return var\n",
    "    \n",
    "    def _bias_variable(self, shape, regularization=True):\n",
    "        \"\"\"Initializer for the bias\"\"\"\n",
    "        \n",
    "        initial = tf.constant_initializer(0.1)\n",
    "        var = tf.get_variable('bias', shape, tf.float32, initializer=initial)\n",
    "        if regularization:\n",
    "            self.regularizers.append(tf.nn.l2_loss(var))\n",
    "        return var\n",
    "    \n",
    "    def _h_variable(self, shape, regularization=False, name=''):\n",
    "        \"\"\"Initializer for the zoom parameter h\"\"\"\n",
    "        \n",
    "        initial = tf.random_uniform_initializer()\n",
    "        var = tf.get_variable('h'+name, shape, tf.float32, initializer=initial)\n",
    "        if regularization:\n",
    "            self.regularizers.append(tf.nn.l2_loss(var))\n",
    "        return var\n",
    "\n",
    "    def frobenius_norm(self, tensor): \n",
    "        \"\"\"Computes the frobenius norm for a given laplacian\"\"\"\n",
    "        \n",
    "        square_tensor = tf.square(tensor)\n",
    "        tensor_sum = tf.reduce_sum(square_tensor)\n",
    "        frobenius_norm = tf.sqrt(tensor_sum)\n",
    "        return frobenius_norm\n",
    "    \n",
    "    def compute_sparse_D_inv_indices(self, M):\n",
    "        \"\"\"Computes the indices required for constructing a sparse version of D^-1.\"\"\"\n",
    "        \n",
    "        idx_main_diag = np.tile(np.expand_dims(np.arange(0, 2*M),1), [1, 2])\n",
    "        idx_diag_ur = np.concatenate([np.expand_dims(np.arange(0, M),1), np.expand_dims(np.arange(0, M)+M,1)], 1)\n",
    "        idx_diag_ll = np.concatenate([np.expand_dims(np.arange(0, M)+M,1), np.expand_dims(np.arange(0, M),1)], 1)\n",
    "        idx = np.concatenate([idx_main_diag, idx_diag_ur, idx_diag_ll], 0)\n",
    "        return idx  \n",
    "    \n",
    "    def compute_sparse_R_indices(self, L_off_diag, M):\n",
    "        \"\"\"Computes the indices required for constructing a sparse version of R.\"\"\"\n",
    "        \n",
    "        idx_L = np.asarray(np.where(L_off_diag)).T\n",
    "        idx_L_sh = idx_L + np.expand_dims(np.asarray([M,M]),0)\n",
    "        idx = np.concatenate([idx_L, idx_L_sh])\n",
    "        return idx\n",
    "    \n",
    "    def compute_sparse_numerator_projection_indices(self, L, M):\n",
    "        \"\"\"Computes the indices required for constructing the numerator projection sparse matrix.\"\"\"\n",
    "        \n",
    "        idx_L = np.asarray(np.where(L)).T\n",
    "        idx_L_sh = idx_L + np.expand_dims(np.asarray([M,M]),0)\n",
    "        idx_diag_ur = np.concatenate([np.expand_dims(np.arange(0, M),1), np.expand_dims(np.arange(0, M)+M,1)], 1)\n",
    "        idx_diag_ll = np.concatenate([np.expand_dims(np.arange(0, M)+M,1), np.expand_dims(np.arange(0, M),1)], 1)\n",
    "        idx = np.concatenate([idx_L, idx_L_sh, idx_diag_ur, idx_diag_ll])\n",
    "        return idx\n",
    "    \n",
    "    def cayleyConv(self, x, L_np, Fout, K): \n",
    "        \"\"\"Applies chebyshev polynomials over the graph.\"\"\"\n",
    "        \n",
    "        M, Fin = x.get_shape()[1:] # M the number of samples in the images, Fin the number of features\n",
    "        M, Fin = int(M), int(Fin)\n",
    "        N = tf.shape(x)[0] # N is the number of images\n",
    "        \n",
    "        # Applies cayley transform by means of Jacobi method.\n",
    "        diag_L_np = np.diag(L_np)  # vector containing the diagonal of L\n",
    "        L_off_diag_np = L_np - np.diag(diag_L_np) # off-diagonal entries of L \n",
    "        \n",
    "        list_x_pos_exp = [tf.cast(tf.expand_dims(x,0), 'complex64')] # 1 x N x M x F\n",
    "        \n",
    "        for iii in range(self.n_h):  # for every zoom parameter we want to use (typically one).\n",
    "            h = self._h_variable([1,1], regularization=False, name='_h%f' % iii)\n",
    "            self.list_h.append(h)\n",
    "            \n",
    "            # Computes matrices required by Jacobi (https://en.wikipedia.org/wiki/Jacobi_method)\n",
    "            \n",
    "            # To make things more efficient we reprent a complex vector of shape M as real vector of shape 2*M\n",
    "            # where the first M values represent real coefficients while the second M the imaginary ones.\n",
    "            # All the matrices here defined are computed according to such notation (it allows to use sparse matrices\n",
    "            # with TF with complex values).\n",
    "            \n",
    "            # ************************** COMPUTES numerator projection **************************\n",
    "            idx = self.compute_sparse_numerator_projection_indices(L_np, M)\n",
    "            \n",
    "            vals_L = tf.squeeze(h*L_np[np.where(L_np)])\n",
    "            vals = tf.concat([vals_L, vals_L, tf.ones([M,]), -tf.ones([M,])], 0)\n",
    "            \n",
    "            cayley_op_neg_sp = tf.SparseTensor(idx, vals, [M*2, M*2])\n",
    "            cayley_op_neg_sp = tf.sparse_reorder(cayley_op_neg_sp)\n",
    "        \n",
    "            # ************************** COMPUTES D **************************\n",
    "            D_real = tf.squeeze(h*diag_L_np)\n",
    "            D = tf.complex(D_real, tf.ones_like(D_real))\n",
    "            D_inv = tf.pow(D, -tf.ones_like(D)) # vector of M elements <- diagonal of D^-1\n",
    "            \n",
    "            idx = self.compute_sparse_D_inv_indices(M)\n",
    "            vals = tf.concat([tf.real(D_inv), tf.real(D_inv), -tf.imag(D_inv), tf.imag(D_inv)], 0)\n",
    "            \n",
    "            D_inv_ext_sp = tf.SparseTensor(idx, vals, [M*2, M*2])\n",
    "            D_inv_ext_sp = tf.sparse_reorder(D_inv_ext_sp)\n",
    "            \n",
    "            # ************************** COMPUTES R **************************\n",
    "            idx = self.compute_sparse_R_indices(L_off_diag_np, M)\n",
    "            \n",
    "            vals_L = tf.squeeze(h*L_off_diag_np[np.where(L_off_diag_np)])\n",
    "            vals = tf.concat([vals_L, vals_L], 0)\n",
    "            \n",
    "            R_sp = tf.SparseTensor(idx, vals, [M*2, M*2])\n",
    "            R_sp = tf.sparse_reorder(R_sp)\n",
    "            \n",
    "            # Applies Jacobi method\n",
    "            c_transform = tf.transpose(x, [1,0,2]) # shape = M, N, F\n",
    "            c_transform = tf.reshape(c_transform, [M, -1]) # shape = M, N*F\n",
    "            last_sol = tf.concat([c_transform, tf.zeros_like(c_transform)],0)\n",
    "            for k in range(K):  # for every order of our polynomial\n",
    "                \n",
    "                # Jacobi initialization\n",
    "                b = tf.sparse_tensor_dense_matmul(cayley_op_neg_sp, last_sol) # shape = M, N*F\n",
    "                a = tf.sparse_tensor_dense_matmul(D_inv_ext_sp, b) # shape = M, N*F\n",
    "                \n",
    "                # Jacobi iterations\n",
    "                cond = lambda i, _: tf.less(i, self.num_jacobi_iter)\n",
    "                body = lambda i, c_sol: [tf.add(i, 1), a  - tf.sparse_tensor_dense_matmul(D_inv_ext_sp, \n",
    "                                                                                          tf.sparse_tensor_dense_matmul(R_sp, c_sol))]\n",
    "                \n",
    "                c_sol = tf.while_loop(cond, body, [0, a], parallel_iterations=1, swap_memory=True)\n",
    "                c_sol = c_sol[-1]\n",
    "                    \n",
    "                # Constructs and saves the final complex matrices\n",
    "                c_sol_complex = tf.complex(c_sol[:M,:], c_sol[M:, :]) #M x N*F\n",
    "                c_sol_reshaped = tf.reshape(c_sol_complex, [M, -1, Fin])\n",
    "                c_sol_reshaped = tf.transpose(c_sol_reshaped, [1, 0, 2]) #N x M x F\n",
    "                list_x_pos_exp.append(tf.expand_dims(c_sol_reshaped,0)) #1 x N x M x Flist_x_pos_exp\n",
    "                \n",
    "                last_sol = c_sol\n",
    "        x_pos_exp = tf.concat(list_x_pos_exp, 0) # shape = n_h*K x N x M x Fin\n",
    "        x_pos_exp = tf.transpose(x_pos_exp, [1,2,0,3])  #N x M x n_h*K x Fin\n",
    "        x_pos_exp = tf.reshape(x_pos_exp, [N*M, -1]) #N*M x 2*K*Fin\n",
    "        \n",
    "        real_conv_weights = self._weight_variable([Fin*(self.n_h*K+1), Fout], regularization=False, name='_real')#tf.ones([Fin*(self.n_h*K+1), Fout])#self._weight_variable([Fin*(self.n_h*K+1), Fout], regularization=False, name='_real')\n",
    "        imag_conv_weights = self._weight_variable([Fin*(self.n_h*K+1), Fout], regularization=False, name='_imag')#tf.ones([Fin*(self.n_h*K+1), Fout])#self._weight_variable([Fin*(self.n_h*K+1), Fout], regularization=False, name='_imag')\n",
    "        \n",
    "        W_pos_exp = tf.complex(real_conv_weights, -imag_conv_weights)\n",
    "        \n",
    "        x_pos_exp_filt = tf.matmul(x_pos_exp, W_pos_exp)\n",
    "        \n",
    "        x_filt = 2*tf.real(x_pos_exp_filt)\n",
    "        return tf.reshape(x_filt, [N, M, Fout])\n",
    "\n",
    "\n",
    "    def b1relu(self, x): #sums a bias and applies relu\n",
    "        \"\"\"Bias and ReLU. One bias per filter.\"\"\"\n",
    "        N, M, F = x.get_shape()\n",
    "        b = self._bias_variable([1, 1, int(F)], regularization=False)\n",
    "        return tf.nn.relu(x + b) #add the bias to the convolutive layer\n",
    "\n",
    "\n",
    "    def mpool1(self, x, p): #efficient pooling realized thanks to the reordering of the laplacians we have done a priori\n",
    "        \"\"\"Max pooling of size p. Should be a power of 2.\"\"\"\n",
    "        if p > 1:\n",
    "            x = tf.expand_dims(x, 3)  # N x M x F x 1\n",
    "            x = tf.nn.max_pool(x, ksize=[1,p,1,1], strides=[1,p,1,1], padding='SAME')\n",
    "            return tf.squeeze(x, [3])  # N x M/p x F\n",
    "        else:\n",
    "            return x\n",
    "    \n",
    "\n",
    "    def b1relu(self, x): #sums a bias and applies relu\n",
    "        \"\"\"Bias and ReLU. One bias per filter.\"\"\"\n",
    "        N, M, F = x.get_shape()\n",
    "        b = self._bias_variable([1, 1, int(F)], regularization=False)\n",
    "        return tf.nn.relu(x + b) #add the bias to the convolutive layer\n",
    "\n",
    "\n",
    "    def mpool1(self, x, p): #efficient pooling realized thanks to the reordering of the laplacians we have done a priori\n",
    "        \"\"\"Max pooling of size p. Should be a power of 2.\"\"\"\n",
    "        if p > 1:\n",
    "            x = tf.expand_dims(x, 3)  # N x M x F x 1\n",
    "            x = tf.nn.max_pool(x, ksize=[1,p,1,1], strides=[1,p,1,1], padding='SAME')\n",
    "            return tf.squeeze(x, [3])  # N x M/p x F\n",
    "        else:\n",
    "            return x\n",
    "\n",
    "    def fc(self, x, Mout, relu=True):\n",
    "        \"\"\"Fully connected layer with Mout features.\"\"\"\n",
    "        N, Min = x.get_shape()\n",
    "        W = self._weight_variable([int(Min), Mout], regularization=True)\n",
    "        b = self._bias_variable([Mout], regularization=True)\n",
    "        x = tf.matmul(x, W) + b\n",
    "        return tf.nn.relu(x) if relu else x\n",
    "    \n",
    "    #function used for extracting the result of our model\n",
    "    def _inference(self, x, dropout): #definition of the model\n",
    "        \n",
    "        # Graph convolutional layers.\n",
    "        x = tf.expand_dims(x, 2)  # N x M x F=1\n",
    "        j = 0\n",
    "        self.list_h = list()\n",
    "        for i in range(len(self.p)):\n",
    "            with tf.variable_scope('cgconv{}'.format(i+1)):\n",
    "                with tf.name_scope('filter'):\n",
    "                    x = self.cayleyConv(x, self.L_np[i*2], self.F[i], self.K[i])\n",
    "                    if (i==0):\n",
    "                        self.debug = x\n",
    "                with tf.name_scope('bias_relu'):\n",
    "                    x = self.b1relu(tf.cast(tf.real(x), 'float32'))\n",
    "                with tf.name_scope('pooling'):\n",
    "                    x = self.mpool1(x, self.p[i])\n",
    "                    \n",
    "            j += int(np.log2(self.p[i])) if self.p[i] > 1 else 0\n",
    "        \n",
    "        # Fully connected hidden layers.\n",
    "        _, M, F = x.get_shape()\n",
    "        x = tf.reshape(x, [-1, int(M*F)])  # N x M\n",
    "        for i,M in enumerate(self.M[:-1]): #apply a fully connected layer for each layer defined in M\n",
    "                                           #(we discard the last value in M since it contains the number of classes we have\n",
    "                                           #to predict)\n",
    "            with tf.variable_scope('fc{}'.format(i+1)):\n",
    "                x = self.fc(x, M)\n",
    "                x = tf.nn.dropout(x, dropout)\n",
    "        \n",
    "        # Logits linear layer, i.e. softmax without normalization.\n",
    "        with tf.variable_scope('logits'):\n",
    "            x = self.fc(x, self.M[-1], relu=False)\n",
    "        return x\n",
    "    \n",
    "    def __init__(self, p, K, F, M, M_0, batch_size, num_jacobi_iter, L,\n",
    "                 decay_steps, decay_rate, learning_rate=1e-4, momentum=0.9, regularization=5e-4, clip_norm=1e1,\n",
    "                 idx_gpu = '/gpu:0'):\n",
    "        self.regularizers = list() #list of regularization l2 loss for multiple variables\n",
    "        self.n_h = 1\n",
    "        self.num_jacobi_iter = num_jacobi_iter\n",
    "        self.p = p #dimensions of the pooling layers\n",
    "        self.K = K #List of polynomial orders, i.e. filter sizes or number of hops\n",
    "        self.F = F #Number of features of convolutional layers\n",
    "        \n",
    "        self.M = M #Number of neurons in fully connected layers\n",
    "        \n",
    "        self.M_0 = M_0 #number of elements in the first graph \n",
    "        \n",
    "        self.batch_size = batch_size\n",
    "        \n",
    "        #definition of some learning parameters\n",
    "        self.decay_steps = decay_steps\n",
    "        self.decay_rate = decay_rate\n",
    "        self.learning_rate = learning_rate\n",
    "        self.regularization = regularization\n",
    "        \n",
    "        with tf.Graph().as_default() as g:\n",
    "                self.graph = g\n",
    "                tf.set_random_seed(0)\n",
    "                with tf.device(idx_gpu):\n",
    "                        #definition of placeholders\n",
    "                        self.L_np = [c_L.toarray().astype('float32') for c_L in L]\n",
    "                        self.ph_data = tf.placeholder(tf.float32, (self.batch_size, M_0), 'data')\n",
    "                        self.ph_labels = tf.placeholder(tf.int32, (self.batch_size), 'labels')\n",
    "                        self.ph_dropout = tf.placeholder(tf.float32, (), 'dropout')\n",
    "                    \n",
    "                        #Model construction\n",
    "                        self.logits = self._inference(self.ph_data, self.ph_dropout)\n",
    "                        \n",
    "                        #Definition of the loss function\n",
    "                        with tf.name_scope('loss'):\n",
    "                            self.cross_entropy = tf.nn.sparse_softmax_cross_entropy_with_logits(logits=self.logits, labels=self.ph_labels)\n",
    "                            self.cross_entropy = tf.reduce_mean(self.cross_entropy)\n",
    "                        with tf.name_scope('regularization'):\n",
    "                            self.regularization *= tf.add_n(self.regularizers)\n",
    "                        self.loss = self.cross_entropy + self.regularization\n",
    "                        \n",
    "                        #Solver Definition\n",
    "                        with tf.name_scope('training'):\n",
    "                            # Learning rate.\n",
    "                            global_step = tf.Variable(0, name='global_step', trainable=False) #used for counting how many iterations we have done\n",
    "                            if decay_rate != 1: #applies an exponential decay of the lr wrt the number of iterations done\n",
    "                                learning_rate = tf.train.exponential_decay(\n",
    "                                        learning_rate, global_step, decay_steps, decay_rate, staircase=True)\n",
    "                            # Optimizer.\n",
    "                            if momentum == 0:\n",
    "                                optimizer = tf.train.GradientDescentOptimizer(learning_rate)\n",
    "                            else: #applies momentum for increasing the robustness of the gradient \n",
    "                                optimizer = tf.train.MomentumOptimizer(learning_rate, momentum)\n",
    "                            #grads = optimizer.compute_gradients(self.loss)\n",
    "                            tvars = tf.trainable_variables()\n",
    "                            #grads, _ = tf.clip_by_global_norm(tf.gradients(self.loss, tvars), clip_norm)\n",
    "                            grads, variables = zip(*optimizer.compute_gradients(self.loss))\n",
    "                            grads, _ = tf.clip_by_global_norm(grads, clip_norm)\n",
    "                            self.op_gradients = optimizer.apply_gradients(zip(grads, variables), \n",
    "                                                                          global_step=global_step)\n",
    "                            \n",
    "                        #Computation of the norm gradients (useful for debugging)\n",
    "                        self.var_grad = tf.gradients(self.loss, tf.trainable_variables())\n",
    "                        self.norm_grad = self.frobenius_norm(tf.concat([tf.reshape(g, [-1]) for g in self.var_grad], 0))\n",
    "\n",
    "                        #Extraction of the predictions and computation of accuracy\n",
    "                        self.predictions = tf.cast(tf.argmax(self.logits, dimension=1), tf.int32)\n",
    "                        self.accuracy = 100 * tf.contrib.metrics.accuracy(self.predictions, self.ph_labels)\n",
    "        \n",
    "                        # Create a session for running Ops on the Graph.\n",
    "                        config = tf.ConfigProto(allow_soft_placement = True)\n",
    "                        config.gpu_options.allow_growth = True\n",
    "                        self.session = tf.Session(config=config)\n",
    "\n",
    "                        # Run the Op to initialize the variables.\n",
    "                        init = tf.global_variables_initializer()\n",
    "                        self.session.run(init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Training & testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Convolutional parameters\n",
    "p = [4, 4]   # Dimensions of the pooling layers\n",
    "K = [12, 12] # List of polynomial orders, i.e. filter sizes or number of hops\n",
    "F = [32, 64] # Number of features of convolutional layers\n",
    "\n",
    "#FC parameters\n",
    "C = max(train_labels) + 1 # Number of classes we have\n",
    "M = [512, C] # Number of neurons in fully connected layers\n",
    "\n",
    "#Solver parameters\n",
    "batch_size = 100\n",
    "decay_steps = train_data.shape[0] / batch_size # number of steps to do before decreasing the learning rate\n",
    "decay_rate = 0.95\n",
    "learning_rate = 0.01\n",
    "momentum = 0.9\n",
    "regularization = 5e-4\n",
    "\n",
    "# Definition of keep probabilities for dropout layers\n",
    "dropout_training = 0.5\n",
    "dropout_val_test = 1.0\n",
    "\n",
    "num_jacobi_iter = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Construction of the learning obj\n",
    "M_0 = L[0].shape[0] # number of elements in the first graph\n",
    "learning_obj = CayleyNet(p, K, F, M, M_0, batch_size, num_jacobi_iter, L,\n",
    "                         decay_steps, decay_rate,\n",
    "                         learning_rate=learning_rate, regularization=regularization,\n",
    "                         momentum=momentum)#, clip_norm=100)\n",
    "\n",
    "# definition of overall number of training iterations and validation frequency\n",
    "num_iter_val = 600\n",
    "num_total_iter_training = 21000\n",
    "\n",
    "num_iter = 0\n",
    "\n",
    "list_training_loss = list()\n",
    "list_training_norm_grad = list()\n",
    "list_val_accuracy = list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#training and validation\n",
    "indices = collections.deque() # queue containing a permutation of the training indexes\n",
    "for k in range(num_iter, num_total_iter_training):\n",
    "\n",
    "    #Construction of the training batch\n",
    "    if len(indices) < batch_size: # Be sure to have used all the samples before using one a second time.\n",
    "        indices.extend(np.random.permutation(train_data.shape[0])) #reinitialize the queue of indices\n",
    "    idx = [indices.popleft() for i in range(batch_size)] #extract the current batch of samples\n",
    "\n",
    "    #data extraction\n",
    "    batch_data, batch_labels = train_data[idx,:], train_labels[idx] \n",
    "\n",
    "    feed_dict = {learning_obj.ph_data: batch_data, \n",
    "                 learning_obj.ph_labels: batch_labels, \n",
    "                 learning_obj.ph_dropout: dropout_training}\n",
    "\n",
    "    #Training\n",
    "    tic = time.time()\n",
    "    _, current_training_loss, norm_grad = learning_obj.session.run([learning_obj.op_gradients, \n",
    "                                                                    learning_obj.loss, \n",
    "                                                                    learning_obj.norm_grad], feed_dict = feed_dict) \n",
    "    training_time = time.time() - tic\n",
    "\n",
    "    list_training_loss.append(current_training_loss)\n",
    "    list_training_norm_grad.append(norm_grad)\n",
    "    if (np.mod(num_iter, num_iter_val)==0): #validation\n",
    "        msg = \"[TRN] iter = %03i, cost = %3.2e, |grad| = %.2e (%3.2es)\" \\\n",
    "                    % (num_iter, list_training_loss[-1], list_training_norm_grad[-1], training_time)\n",
    "        print msg\n",
    "\n",
    "        #Validation Code\n",
    "        tic = time.time()\n",
    "        val_accuracy = 0\n",
    "        for begin in range(0, val_data.shape[0], batch_size):\n",
    "            end = begin + batch_size\n",
    "            end = min([end, val_data.shape[0]])\n",
    "\n",
    "            #data extraction\n",
    "            batch_data = np.zeros((end-begin, val_data.shape[1]))\n",
    "            batch_data = val_data[begin:end,:]\n",
    "            batch_labels = np.zeros(batch_size)\n",
    "            batch_labels[:end-begin] = val_labels[begin:end]\n",
    "\n",
    "            feed_dict = {learning_obj.ph_data: batch_data, \n",
    "                         learning_obj.ph_labels: batch_labels,\n",
    "                         learning_obj.ph_dropout: dropout_val_test}\n",
    "\n",
    "            batch_accuracy = learning_obj.session.run(learning_obj.accuracy, feed_dict)\n",
    "            val_accuracy += batch_accuracy*batch_data.shape[0]\n",
    "        val_accuracy = val_accuracy/val_data.shape[0]\n",
    "\n",
    "        val_time = time.time() - tic\n",
    "        msg = \"[VAL] iter = %03i, acc = %4.2f (%3.2es)\" % (num_iter, val_accuracy, val_time)\n",
    "        print msg\n",
    "    num_iter += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Test code\n",
    "tic = time.time()\n",
    "test_accuracy = 0\n",
    "for begin in range(0, test_data.shape[0], batch_size):\n",
    "    end = begin + batch_size\n",
    "    end = min([end, test_data.shape[0]])\n",
    "\n",
    "    batch_data = np.zeros((end-begin, test_data.shape[1]))\n",
    "    batch_data = test_data[begin:end,:]\n",
    "\n",
    "    feed_dict = {learning_obj.ph_data: batch_data, learning_obj.ph_dropout: 1}\n",
    "\n",
    "    batch_labels = np.zeros(batch_size)\n",
    "    batch_labels[:end-begin] = test_labels[begin:end]\n",
    "    feed_dict[learning_obj.ph_labels] = batch_labels\n",
    "\n",
    "    batch_accuracy = learning_obj.session.run(learning_obj.accuracy, feed_dict)\n",
    "    test_accuracy += batch_accuracy*batch_data.shape[0]\n",
    "test_accuracy = test_accuracy/test_data.shape[0]\n",
    "test_time = time.time() - tic\n",
    "msg = \"[TST] iter = %03i, acc = %4.2f (%3.2es)\" % (num_iter, test_accuracy, test_time)\n",
    "print msg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}