{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"b55a88ba-cc45-e017-8115-3592bcab614c"},"source":"This script is for digit recognition with 2-convolutions-CNN\n\nHowever, it takes long time to learn with full data set. Now, 'train_cnt' is only 100, and 'test_cnt' is only 10\n\nUsing full data, my score is 0.87586. Parameter may need to be tuned."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1048f06c-5f72-bda2-0eb2-1dc53284a40b"},"outputs":[],"source":"import numpy as np\nimport math\nfrom sklearn.preprocessing import LabelBinarizer\nimport time\nimport csv\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd\n%matplotlib inline"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0080e46e-966e-3149-256e-f458fb8a8e16"},"outputs":[],"source":"#Prepare function\ndef gaborc(H,k): #H is size of filter, l is direction, k is coeff for sigma\n    buf = np.zeros([H, H])\n    l = 2*math.pi*math.cos(k)/H\n    m = 2*math.pi*math.sin(k)/H\n    n = 1/math.sqrt(2*math.pi*H*H)\n    for i in range(H):\n        for j in range(H):\n            buf[i,j] = math.cos((l*(i-H/2)+m*(H/2-j))) * n * math.exp(-((i-H/2)**2+(H/2-j)**2)/(2*H*H))\n    return buf\n\ndef gabors(H,k): #H is size of filter, l is direction, k is coeff for sigma\n    buf = np.zeros([H, H])\n    l = 2*math.pi*math.cos(k)/H\n    m = 2*math.pi*math.sin(k)/H\n    n = 1/math.sqrt(2*math.pi*H)\n    for i in range(H):\n        for j in range(H):\n            buf[i,j] = math.sin((l*(i-H/2)+m*(H/2-j))) * n * math.exp(-((i-H/2)**2+(H/2-j)**2)/(2*H*H))\n    return buf\n\ndef zpadding(x, n, m, p): # x is (nxm)x1 array, p is addtional zeros\n    y = x.reshape([n,m])    \n    ab = np.zeros([p,m])\n    side = np.zeros([n+2*p, p])\n    y = np.r_[ab,y,ab]\n    y = np.c_[side,y,side]\n    return y.reshape([1,(n+2*p)*(m+2*p)])\n    \ndef Wconv(h, n, m): # h is filter matrix. nxm is conv-ed matrix\n    nh = len(h)    \n    nc = n-nh+1\n    mc = m-nh+1\n    buf = np.zeros([nc*mc,n*m])\n    for i in range(mc):\n        for j in range(nc):\n            for p in range(nh):\n                for q in range(nh): \n                    buf[i+j*mc,i+j*m+p+q*m] = h[p,q]\n    return buf\n    \ndef pooling(x,n,m,p,s): # p is size of pooling, s is stride, maximum pooling\n    nc = int((n-p)/s)+1\n    mc = int((m-p)/s)+1\n    buf = np.zeros([nc,mc])\n    x = x.reshape([n,m]) \n    w = np.zeros([nc*mc,n*m])\n    for i in range(nc):\n        for j in range(mc):\n            maxr = 0\n            maxq = 0\n            for r in range(p):\n                for q in range(p):\n                    if (buf[i,j] < x[i*s+r,j*s+q]):\n                        maxr = r #update \n                        maxq = q #update                 \n            buf[i,j] = x[i*s+maxr,j*s+maxq]\n            w[i*mc+j,(i*s+maxr)*m+j*s+maxq] = 1 # w is pick up matrix\n    buf = buf.reshape([nc*mc])\n    return (buf, w) \n    \ndef relu(x): # x is 1-dim array\n    n = x.size\n    buf = np.zeros([n])    \n    pos = np.zeros([n])\n    for i in range(n):\n        if (x[0,i]>0):\n            buf[i] = x[0,i]\n            pos[i] = 1 # pos is activated position\n    return (buf, pos) \n\ndef connorm(x):\n    c = 1e-3\n    return (x - np.mean(x))/math.sqrt((c**2+ np.std(x**2)))\n\ndef softmax(y):\n    buf = np.exp(y)\n    return buf/np.sum(buf) "},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"08340ec3-a829-6cfc-eaea-8a02c1e20fa1"},"outputs":[],"source":"#It takes long time to learn with full data set\n#Now, 'train_cnt' is only 100, and 'test_cnt' is only 10\ndef main ():\n# load mnist data from kaggle dataset\n    data = np.loadtxt(\"../input/train.csv\",delimiter=\",\", skiprows=1)\n    T_n, X_train = np.split(data, [1], axis=1) #data: [label, px1, px2,,,]\n    X_train /= X_train.max() #adjust maximum as 1\n    T_train = LabelBinarizer().fit_transform(T_n) # transoform to 1-of-k form   \n    N = len(T_train)\n    D = len(X_train[0,:]) # dimension of image\n\n# parameters\n    cl = 10 # num of class\n    d = int(math.sqrt(D))\n    pd = 4 # zero padding size\n    st = 2 # stride\n    pl = 3 # pooling size\n    eta = 0.01 # learning coeff\n    train_cnt = 100#N # train data (by one image)\n\n# first convolution filter [dh1 x dh1]matrix x nh1\n    dh1 = 5 # dim of filter\n    nh1 = 6 # number of filter\n    H1 = np.zeros([nh1,dh1,dh1])\n    for i in range(int(nh1/2)):\n        H1[i,:,:] = gaborc(dh1,i*math.pi/nh1*2)\n        H1[i+int(nh1/2),:,:] = gabors(dh1,i*math.pi/nh1*2)\n        \n# second convolution filter [dh2 x dh2]matrix x nh2\n    dh2 = 3 # dim of filter\n    nh2 = 4 # number of filter\n    H2 = np.zeros([nh2,dh2,dh2])\n    for i in range(int(nh2/2)):\n        H2[i,:,:] = gaborc(dh2,i*math.pi/nh2*2)\n        H2[i+int(nh2/2),:,:] = gabors(dh2,i*math.pi/nh2*2)\n        \n# dimensions for each propagation\n    d_pd = d+2*pd # 1st padding 28 + 2 * 4\n    d_cv = d_pd - dh1 + 1 # 1st convolution 36 - 11 +1\n    d_pl = int((d_cv - pl)/st) + 1 # 1st pooling\n    \n    d_pd2 = d_pl+2*pd # 2nd padding 28 + 2 * 4\n    d_cv2 = d_pd2 - dh2 + 1 # 2nd convolution 36 - 11 +1\n    d_pl2 = int((d_cv2 - pl)/st) + 1 # 2nd pooling\n                                \n# learning\n# in/out of layer is 1-dim array\n    w2 = 0.01 * np.random.random(nh2*nh1*d_pl2*d_pl2*cl).reshape([cl,nh1*nh2*d_pl2*d_pl2]) # all connection 10 <- 12x12x12  \n    bo = 0.01 * np.random.random(cl) #bias for all connection\n    bc1 = 0.01 * np.random.random(nh1) # bias for convolution\n    bc2 = 0.01 * np.random.random(nh2) # bias for convolution\n    correct = 0\n\n# prepare matrix - 1st convolution layer                                \n    wc1 = np.zeros([nh1, d_cv*d_cv, d_pd*d_pd]) # convolution matrix 36x36 -> 26x26\n    conv1 = np.zeros([nh1,d_cv*d_cv])\n    pos1 = np.zeros([nh1,d_cv*d_cv])\n    B1 = np.zeros([nh1,d_cv*d_cv])\n    pool1 = np.zeros([nh1,d_pl*d_pl])\n    wp1 = np.zeros([nh1,d_pl*d_pl,d_cv*d_cv]) #pooling matrix\n    conm1 = np.zeros([nh1, d_pl*d_pl])\n                   \n# second convolution layer\n    zpd2 = np.zeros([nh1,d_pd2*d_pd2])\n    wc2 = np.zeros([nh2, nh1, d_cv2*d_cv2, d_pd2*d_pd2]) # convolution matrix 36x36 -> 26x26\n    conv2 = np.zeros([nh2, nh1, d_cv2*d_cv2])\n    pos2 = np.zeros([nh2, nh1, d_cv2*d_cv2])\n    B2 = np.zeros([nh2, d_cv2*d_cv2])\n    pool2 = np.zeros([nh2, nh1, d_pl2*d_pl2])\n    wp2 = np.zeros([nh2, nh1, d_pl2*d_pl2, d_cv2*d_cv2]) #pooling matrix\n    conm2 = np.zeros([nh2, nh1, d_pl2*d_pl2])        \n\n# start learning -- learn one by one (not batch learning)\n    start = time.time()\n    Y_train = []\n\n# forward propergation  \n    for a in range(train_cnt):\n        m = int(a) #pick up parameter\n        x = X_train[m]\n        zpd1 = zpadding(x,d,d,pd)\n        for i in range(nh1):\n            for j in range(d_cv*d_cv):\n                B1[i,j] = bc1[i]\n        for i in range(nh2):\n            for j in range(d_cv2*d_cv2):\n                B2[i,j] = bc2[i]\n        for k in range(nh1): # 1st convolution layers\n            wc1[k,:,:] = Wconv(H1[k,:,:],d_pd,d_pd)\n            conv_dum1 = wc1[k,:,:].dot(zpd1[0,:]).reshape([1,d_cv*d_cv]) + B1[k,:]\n            conv1[k,:], pos1[k,:] = relu(conv_dum1)  # 1st convoluition\n            pool1[k,:], wp1[k,:,:] = pooling(conv1[k,:],d_cv,d_cv, pl, st) # 1st pooling \n            conm1[k,:]= connorm(pool1[k,:])\n            zpd2[k,:] = zpadding(conm1[k,:],d_pl,d_pl,pd)\n            for l in range(nh2): # 2nd convolution layers\n                wc2[l,k,:,:] = Wconv(H2[l,:,:],d_pd2,d_pd2)\n                conv_dum2 = wc2[l,k,:,:].dot(zpd2[0,:]).reshape([1,d_cv2*d_cv2]) + B2[l,:] #2nd Relu\n                conv2[l,k,:], pos2[l,k,:] = relu(conv_dum2)  # 2nd convoluition\n                pool2[l,k,:], wp2[l,k,:,:] = pooling(conv2[l,k,:],d_cv2,d_cv2, pl, st) # 2nd pooling       \n                conm2[l,k,:]= connorm(pool2[l,k,:])\n                \n        conm_rs = conm2.reshape([nh2*nh1*d_pl2*d_pl2])\n        nw = conm_rs.size\n        out_dum = w2.dot(conm_rs[:])+bo[:]\n        y = softmax(out_dum) # 10-dim\n        \n#pickup simulated label\n        maxc = 0\n        buf = 0\n        for i in range(cl):\n            if y[i] > buf:\n                buf = y[i]\n                maxc = i\n        Y_train.append(maxc)    \n        if maxc == T_n[a]:\n            correct += 1\n            \n# back propagation\n# back propagation for 2nd layer\n        d1 = y-T_train[m,:] #10 dim, diff of output\n        d2 = w2.T.dot(d1).reshape([nh2,nh1,d_pl2*d_pl2]) #diff of 2nd pooling  \n        d3 = np.zeros([nh2,nh1,d_cv2*d_cv2]) # diff of 2nd Relu\n        for k in range(nh1):\n            for l in range(nh2):\n                d3[l,k,:] = wp2[l,k,:,:].T.dot(d2[l,k,:]) \n        d4 = np.zeros([nh2,nh1,d_cv2*d_cv2]) #diff of 2nd convolution\n        for k in range(nh1):\n            for l in range(nh2):\n                d4[l,k,:] = pos2[l,k,:]*d3[l,k,:] \n        d4_rs = d4.reshape([nh2,nh1,d_cv2,d_cv2])   \n        d5 = np.zeros([nh2,dh2,dh2]) #diff of 2nd filters\n\n        dw2 = np.zeros([cl,nw])\n        for i in range(cl):\n            bo[i] -= eta*d1[i]\n            for j in range(nw):\n                dw2[i,j] = d1[i]*conm_rs[j]\n        w2 -= eta*dw2\n        \n        zpd2_rs = zpd2.reshape([nh1,d_pd2,d_pd2])\n        for k in range(nh1):\n            for l in range(nh2):\n                for i in range(d_cv2):\n                    for j in range(d_cv2):\n                        bc2[l] -= eta * d4_rs[l,k,i,j]\n                        for s in range(dh2):\n                            for t in range(dh2):\n                                d5[l,s,t] += d4_rs[l,k,i,j]*zpd2_rs[k,i+s,j+t]\n        H2 -= eta*d5\n        \n#back propagation of 1st layer\n        d6_dum = np.zeros([nh1,d_pd2,d_pd2]) # diff of 2nd padding\n        d6 = np.zeros([nh1,d_pl,d_pl]) #diff of 1st pooling\n        for k in range(nh1):\n            for l in range(nh2):\n                d6_dum[k,:,:] += wc2[l,k,:,:].T.dot(d4[l,k,:]).reshape([d_pd2,d_pd2])\n            for i in range(d_pl):\n                for j in range(d_pl):\n                    d6[k,i,j] = d6_dum[k,i+pd,j+pd]\n        d6_rs = d6.reshape([nh1,d_pl*d_pl])\n        d7 = np.zeros([nh1,d_cv*d_cv]) # diff of 1st Relu\n        for k in range(nh1):\n            d7[k,:] = wp1[k,:,:].T.dot(d6_rs[k,:]) \n        d8 = np.zeros([nh1, d_cv*d_cv]) # diff of 1st convolution\n        for k in range(nh1):        \n            d8[k,:] = pos1[k,:]*d7[k,:] \n        d8_rs = d8.reshape([nh1,d_cv,d_cv]) # 26x26 dim, d of Convolution   \n        d9 = np.zeros([nh1,dh1,dh1]) # diff of 1st filters\n\n        zpd1_rs = zpd1.reshape([d_pd,d_pd])\n        for k in range(nh1):\n            for i in range(d_cv):\n                for j in range(d_cv):\n                    bc1[k] -= eta * d8_rs[k,i,j]\n                    for s in range(dh1):\n                        for t in range(dh1):\n                            d9[k,s,t] +=  d8_rs[k,i,j]*zpd1_rs[i+s,j+t]\n        H1 -= eta*d9\n \n    elapsed_time = time.time() - start\n    print(\"time\", elapsed_time)\n    \n    print('Y_train',Y_train[(train_cnt-10):train_cnt])\n    print('T_train',T_n[(train_cnt-10):train_cnt].T)\n    print('correct',correct/train_cnt)\n\n#check with test data\n    X_test = np.loadtxt(\"../input/test.csv\",delimiter=\",\", skiprows=1)\n    n_test = len(X_test[:,0]) # the amount of Y_test\n    X_test /= X_test.max() #adjust maximum as 1 \n    test_cnt =  10#n_test\n\n    Y_test = np.zeros(test_cnt)\n    for a in range(test_cnt):\n        m = int(a)\n        x = X_test[m]\n        zpd1 = zpadding(x,d,d,pd)   # zero padding 28x28 -> 36x36  \n        for k in range(nh1): # 1st convolution layers\n            wc1[k,:,:] = Wconv(H1[k,:,:],d_pd,d_pd)\n            conv_dum1 = wc1[k,:,:].dot(zpd1[0,:]).reshape([1,d_cv*d_cv]) + B1[k,:]\n            conv1[k,:], pos1[k,:] = relu(conv_dum1)  # 1st convoluition\n            pool1[k,:], wp1[k,:,:] = pooling(conv1[k,:],d_cv,d_cv, pl, st) # 1st pooling \n            conm1[k,:]= connorm(pool1[k,:])\n            zpd2[k,:] = zpadding(conm1[k,:],d_pl,d_pl,pd)\n            for l in range(nh2): # 2nd convolution layers\n                wc2[l,k,:,:] = Wconv(H2[l,:,:],d_pd2,d_pd2)\n                conv_dum2 = wc2[l,k,:,:].dot(zpd2[0,:]).reshape([1,d_cv2*d_cv2]) + B2[l,:] #2nd Relu\n                conv2[l,k,:], pos2[l,k,:] = relu(conv_dum2)  # 2nd convoluition\n                pool2[l,k,:], wp2[l,k,:,:] = pooling(conv2[l,k,:],d_cv2,d_cv2, pl, st) # 2nd pooling       \n                conm2[l,k,:]= connorm(pool2[l,k,:])\n                \n        conm_rs = conm2.reshape([nh2*nh1*d_pl2*d_pl2])\n        nw = conm_rs.size\n        out_dum = w2.dot(conm_rs[:])+bo[:]\n        y = softmax(out_dum) # 10-dim\n        \n        buf = 0\n        for i in range(cl):\n            if y[i] > buf:\n                buf = y[i]\n                maxc = i\n        Y_test[a] = maxc\n        \n    print('out put label for test data, last 10')\n    print(Y_test[(test_cnt-10):(test_cnt)])\n    \n    f = open('out.csv', 'w')\n    writer = csv.writer(f)\n    writer.writerows([['ImageId','Label']])\n    for i in range(test_cnt):\n        writer.writerows([[i+1,int(Y_test[i])]])\n    f.close()\n\n    print('last one of X_test is ',Y_test[test_cnt-1], '. Is it correct??')\n    sns.heatmap(X_test[test_cnt-1,:].reshape(28,28))"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0578d8dc-92e2-c417-c5aa-8b1a74730045"},"outputs":[],"source":"if __name__ == '__main__':\n    main()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b330f431-7ecb-4314-d26b-f74423b5eaab"},"outputs":[],"source":""}],"metadata":{"_change_revision":0,"_is_fork":false,"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":0}