{"nbformat_minor": 0, "metadata": {"language_info": {"file_extension": ".py", "codemirror_mode": {"name": "ipython", "version": 3}, "mimetype": "text/x-python", "nbconvert_exporter": "python", "name": "python", "version": "3.6.0", "pygments_lexer": "ipython3"}, "kernelspec": {"language": "python", "display_name": "Python 3", "name": "python3"}}, "cells": [{"outputs": [], "metadata": {"_execution_state": "idle", "collapsed": false, "_uuid": "e09a686b40e7f55237aeafc1a7c70c4af6ec8f35"}, "execution_count": null, "cell_type": "markdown", "source": "After having implemented a simple 2-layer MLP ([here][1]), I present here a simple implementation of a Convolutional Neural Network with just one convolutional layer and one mean pooling layer. \nThe design follows closely some exercise I found online [here][2], and some useful help  on the implementation of the back propagation [e.g. here][3]. To reach the claimed accuracy (98%), I used 40 filters.\n\n\n  [1]: https://www.kaggle.com/aurelienbl/2-layer-nn-for-the-digits-recognition-97\n  [2]: http://ufldl.stanford.edu/tutorial/\n  [3]: https://grzegorzgwardys.wordpress.com/2016/04/22/8/"}, {"outputs": [], "metadata": {"_uuid": "8bad7d9710cc641d96c99a665b504889afcadca9"}, "execution_count": null, "cell_type": "code", "source": "# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as nm # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nfrom subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n\n# Any results you write to the current directory are saved as output.\nimport matplotlib.pyplot as plt\nfrom scipy import optimize as opti\n"}, {"outputs": [], "metadata": {"_execution_state": "idle", "collapsed": false, "_uuid": "3261a8c181f97236ff7e238f010c0827c5889263"}, "execution_count": null, "cell_type": "code", "source": "## Defining some functions.\n\n## This is the sigmoid activation function\ndef sigmoid(z):\n    return 1./(1+ nm.exp(-z));\n\n## This is the gradient of the sigmoid function. Used in the backprop\ndef sigmoidGradient(z):\n     return sigmoid(z)*(1-sigmoid(z));\n\n    \n### these 2 function are used to go back and forth between all the paramters stacked and flatten\n### and arrays of proper shape and size\n### weights and bias are the parameters of the final dense layer (basixally a MLP woth no hidden layer)\n### and Wf and bf are the parameters of the CONV layer\n### all_units is a tuple that contains the number  of units of the hidden layers of the MLP. Here is empty\ndef Wb2theta(weights, bias, Wf,bf, all_units):\n    #for each layer w then b \n    theta=weights[0].flatten()\n    theta=nm.append(theta, bias[0].flatten())\n    for i in range(1,len(all_units)-1):\n        theta=nm.append(theta,weights[i].flatten())\n        theta=nm.append(theta,bias[i].flatten())\n\n    theta=nm.append(theta, Wf.flatten())\n    theta=nm.append(theta, bf.flatten())\n    return theta\n\n\ndef theta2Wb(theta, all_units, ):\n    \n    w=[]\n    b=[]\n    wf=[]\n    bf=[]\n\n    id=0\n    for i in range(len(all_units)-1):\n        size1=all_units[i+1]\n        size2=all_units[i]\n        w.append(nm.array(theta[id:id+size1*size2]).reshape((size1, size2)))\n        id=id+size1*size2\n        b.append(nm.array(theta[id:id+size1]).reshape((size1, 1)))\n        id=id+size1\n\n    size1=filterDim**2 * numFilters\n    wf=nm.array(theta[id:id+size1].reshape(filterDim,filterDim, numFilters))\n    id=id+size1\n    size1=numFilters\n    bf=nm.array(theta[id:id+size1].reshape(size1))\n\n    return w,b, wf,bf"}, {"outputs": [], "metadata": {"_execution_state": "idle", "collapsed": false, "_uuid": "9b493b808455872d4669dfe9c25c665df1e86c98"}, "execution_count": null, "cell_type": "code", "source": "## More functions. These two are the convolution and mean pooling part\n\ndef convolve_images(filterDim, numFilters, images, W, b):\n    N_images=images.shape[0]\n    imageDim =images.shape[1]\n\n    convDim = imageDim -filterDim +1\n    convolvedFeatures = nm.zeros((convDim, convDim, numFilters, N_images))\n    for i in range(N_images):\n        for j in range(numFilters):\n            convolvedFeatures[:,:, j,i]= signal.convolve2d(images[i, :,:],W[:,:,j], 'valid') + b[j]\n        \n    return convolvedFeatures\n\n\ndef cnnPool(poolDim , convolvedFeatures):\n    numImages = convolvedFeatures.shape[3]\n    numFilters = convolvedFeatures.shape[2];\n    convolvedDim = convolvedFeatures.shape[0];\n\n    pooledFeatures = nm.zeros((convolvedDim / poolDim,\n                               (convolvedDim / poolDim), numFilters, numImages));\n    \n    for i in range(convolvedDim / poolDim):\n        for j in range(convolvedDim / poolDim):\n            pooledFeatures[i,j, :,:]= nm.mean(convolvedFeatures[i* poolDim:(i+1)*poolDim   ,j*(poolDim):(j+1)*poolDim ,  :,:], axis=(0,1))\n            \n    return pooledFeatures\n\n    "}, {"outputs": [], "metadata": {"_execution_state": "idle", "collapsed": false, "_uuid": "d79e33bfa067bcbbc135d040977e79ae92e7532c"}, "execution_count": null, "cell_type": "code", "source": "#### this is the main part of the code; the function that compute the cost function\n### and its gradient\n\ndef costFunction(theta, X, y, all_units, pred=False):\n\n    m=X.shape[0]\n    y=y.reshape(y.size, 1)\n    ### reshaping the input paramters in a more convenient form\n    weights, bias, Wf, bf= theta2Wb(theta,all_units)\n    \n    ## 1 convolution + sigmoid + mean pooling\n\n    convolvedFeatures = convolve_images(filterDim, numFilters, X, Wf, bf);\n    pooledFeatures= cnnPool(poolDim , sigmoid(convolvedFeatures))\n    \n    ## agreggating everything in one big vector \n    denseFeatures=pooledFeatures.reshape(pooledDim*pooledDim*numFilters, m)\n \n    ### adding the ouput to a sofmax layer. This part is a code for a N hidden layer MLP. \n    ### Here N=0, but for simplicity I left the (here useless) loops.\n    a0=denseFeatures\n    \n    zl=[]\n    al=[]\n    \n    al.append(a0)\n    zl.append(0)\n        \n    for i in range( Nh):\n        zl.append(nm.dot(weights[i], al[i]) + bias[i])\n        al.append(f(zl[i+1]))\n    \n    zl.append(nm.dot(weights[Nh], al[Nh]) + bias[Nh])\n    \n    ## In case one just want to compute prediction and no cost are gradient\n    if pred:\n        result=nm.exp(zl[-1])/nm.sum(nm.exp(zl[-1]), axis=0)\n        pred=nm.argmax(result, axis=0)\n        return  pred\n    \n   \n    ## defing the big target array\n    Y=nm.zeros((all_units[-1], m))\n    for k in range(all_units[-1]):\n        id=nm.where(y==k)[0]\n        Y[k,id]=1\n        \n    ## finding the mac of the softmax argument to remove it to prevent overflow\n    mm=nm.max(zl[-1], axis=0)\n    \n    ## cost function\n    J=-(Y*((zl[-1]-mm) - nm.log(nm.sum(nm.exp(zl[-1]-mm), axis=0)))).sum()\n    \n    # J=-(Y*nm.log(nm.exp(zl[-1])/nm.sum(nm.exp(zl[-1]), axis=0))).sum()\n    \n    ### back propagation\n\n    delta_L= -( Y - nm.exp(zl[-1]-mm)/nm.sum(nm.exp(zl[-1]-mm), axis=0)   )\n\n    delta_Lm1 = nm.dot(weights[0].T, delta_L)\n    delta_pool=delta_Lm1.reshape(pooledDim,pooledDim,numFilters,m)\n    \n    delta_pool2=nm.zeros((convDim, convDim, numFilters, m))\n\n    for i in range(m):\n        for j in range(numFilters):\n            delta_pool2[:,:,j,i]=nm.kron(delta_pool[:,:,j,i], nm.ones((poolDim,poolDim)))/ (poolDim**2)*sigmoidGradient(convolvedFeatures[:,:,j,i])\n\n    dJdw2=nm.zeros((filterDim, filterDim, numFilters, m))\n\n    for i in range(m):\n        for j in range(numFilters):\n            dJdw2[:,:,j,i]=signal.convolve2d(delta_pool2[:,:,j,i],nm.rot90(X[i, :,:],2), 'valid')\n\n\n    dJdw=nm.sum(dJdw2, axis=3)\n    dJdb=delta_pool2.sum(axis=(0,1,3))\n    \n    dJdw0=[nm.dot(delta_L,denseFeatures.T),]\n    dJdb0=[delta_L.sum(axis=1).reshape(numClasses, 1),]\n\n    grad=Wb2theta(dJdw0,dJdb0,dJdw, dJdb, all_units )\n\n    return J, grad\n"}, {"outputs": [], "metadata": {"_execution_state": "idle", "collapsed": false, "_uuid": "edfabddce42a6be3ca940f93171d8572013ff170"}, "execution_count": null, "cell_type": "code", "source": "A=pd.read_csv('../input/train.csv').values\nB=pd.read_csv('../input/test.csv').values\n\n## defining the target and the feature matrix.\nX=(A[:, 1:]-128)/128.\nB=(B[:,:]-128)/128.\ny=nm.int32(A[:,0])\n\n\nimageDim = 28;       #  % image dimension\nfilterDim = 7;        #  % filter dimension\nnumFilters = 2     #   % number of feature maps\n\npoolDim = 2; \nconvDim = imageDim -filterDim +1 \npooledDim=convDim/poolDim  \nnumClasses=10\n\n### Parameters for the dense MLP layer.\nNh=0\nNl=Nh+2\nall_units=[]\nall_units.append(nm.int32(pooledDim*pooledDim*numFilters))\nNunits=() \n\nfor  l in range(Nh):\n    all_units.append(nm.int32(Nunits[l]))\n\nall_units.append(numClasses)\n\n## defining the initial weights at random.  The choice for the MLP weights comes from Glorot&Bengio, 2010\nWf = nm.random.randn(filterDim,filterDim, numFilters);\nbf = nm.random.rand(numFilters)*0;\n\nweights=[]\nbias=[]\nfor l  in range(Nl-1):\n    s=nm.sqrt(6)/ nm.sqrt(all_units[l+1] +all_units[l])\n    weights.append(nm.random.rand(all_units[l+1] ,all_units[l]   ) *2*s -s)\n    bias.append(nm.random.rand(all_units[l+1],1)*0)\n\ntheta_init=Wb2theta(weights, bias, Wf,bf, all_units)\n\n\nWf = nm.random.randn(filterDim,filterDim, numFilters);\nbf = nm.random.rand(numFilters)*0;\n\nweights=[]\nbias=[]\nfor l  in range(Nl-1):\n    s=nm.sqrt(6)/ nm.sqrt(all_units[l+1] +all_units[l])\n    weights.append(nm.random.rand(all_units[l+1] ,all_units[l]   ) *2*s -s)\n    bias.append(nm.random.rand(all_units[l+1],1)*0)\n\n\n\ntheta_init=Wb2theta(weights, bias, Wf,bf, all_units)\n\nimages=X.copy()\nm=images.shape[0]\ntheta=theta_init.copy()\n\n### Home-made implemenantation of a SGD. \nv=0\ngamma=0.5\nalpha0=1e-2\nalpha=alpha0\nminibatch=30\nidd=0\nfor i in range(1):\n    r=nm.random.permutation(m)\n    for j in range(m%minibatch):\n        tt=r[j*minibatch:(j+1)*minibatch]\n        J, grad=costFunction(theta, images[tt,:,:], y[tt], all_units)\n    \n        v=gamma*v +alpha * grad\n        theta=theta - v \n        print('i=', i,'j=', j, 'J=',J, 'alpha=', alpha)\n    alpha =alpha/2.\n\n    \n### one the training is done compute the training error:\n#p=costFunction(theta, images, y, all_units, pred=True)\n\n#print(nm.mean(p==y))\n"}], "nbformat": 4}