{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":3004,"databundleVersionId":861823,"sourceType":"competition"}],"dockerImageVersionId":30035,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Simple MNIST NN from scratch\n\nIn this notebook, I implemented a simple two-layer neural network and trained it on the MNIST digit recognizer dataset. It's meant to be an instructional example, through which you can understand the underlying math of neural networks better.\n\nHere's a video I made explaining all the math and showing my progress as I coded the network: https://youtu.be/w8yWXqWQYmU","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom matplotlib import pyplot as plt\nimport matplotlib as mpl","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-06-12T15:31:55.208082Z","iopub.execute_input":"2024-06-12T15:31:55.208456Z","iopub.status.idle":"2024-06-12T15:31:55.213758Z","shell.execute_reply.started":"2024-06-12T15:31:55.208419Z","shell.execute_reply":"2024-06-12T15:31:55.212856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nfrom matplotlib import pyplot as plt\nimport pandas as pd\n\ndata = pd.read_csv('/kaggle/input/digit-recognizer/train.csv')\n  \ndata = np.array(data)\nm, n = data.shape\nnp.random.shuffle(data) # shuffle before splitting into dev and training sets\n\ndata_dev = data[0:1000].T\nY_dev = data_dev[0]\nX_dev = data_dev[1:n]\nX_dev = X_dev / 255.\n\ndata_train = data[1000:m].T\nY_train = data_train[0]\nX_train = data_train[1:n]\nX_train = X_train / 255.\n_,m_train = X_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:31:55.220008Z","iopub.execute_input":"2024-06-12T15:31:55.220400Z","iopub.status.idle":"2024-06-12T15:31:59.586941Z","shell.execute_reply.started":"2024-06-12T15:31:55.220364Z","shell.execute_reply":"2024-06-12T15:31:59.585824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our NN will have a simple two-layer architecture. Input layer $a^{[0]}$ will have 784 units corresponding to the 784 pixels in each 28x28 input image. A hidden layer $a^{[1]}$ will have 10 units with ReLU activation, and finally our output layer $a^{[2]}$ will have 10 units corresponding to the ten digit classes with softmax activation.\n\n**Forward propagation**\n\n$$Z^{[1]} = W^{[1]} X + b^{[1]}$$\n$$A^{[1]} = g_{\\text{ReLU}}(Z^{[1]}))$$\n$$Z^{[2]} = W^{[2]} A^{[1]} + b^{[2]}$$\n$$A^{[2]} = g_{\\text{softmax}}(Z^{[2]})$$\n\n**Backward propagation**\n\n$$dZ^{[2]} = A^{[2]} - Y$$\n$$dW^{[2]} = \\frac{1}{m} dZ^{[2]} A^{[1]T}$$\n$$dB^{[2]} = \\frac{1}{m} \\Sigma {dZ^{[2]}}$$\n$$dZ^{[1]} = W^{[2]T} dZ^{[2]} .* g^{[1]\\prime} (z^{[1]})$$\n$$dW^{[1]} = \\frac{1}{m} dZ^{[1]} A^{[0]T}$$\n$$dB^{[1]} = \\frac{1}{m} \\Sigma {dZ^{[1]}}$$\n\n**Parameter updates**\n\n$$W^{[2]} := W^{[2]} - \\alpha dW^{[2]}$$\n$$b^{[2]} := b^{[2]} - \\alpha db^{[2]}$$\n$$W^{[1]} := W^{[1]} - \\alpha dW^{[1]}$$\n$$b^{[1]} := b^{[1]} - \\alpha db^{[1]}$$\n\n**Vars and shapes**\n\nForward prop\n\n- $A^{[0]} = X$: 784 x m\n- $Z^{[1]} \\sim A^{[1]}$: 10 x m\n- $W^{[1]}$: 10 x 784 (as $W^{[1]} A^{[0]} \\sim Z^{[1]}$)\n- $B^{[1]}$: 10 x 1\n- $Z^{[2]} \\sim A^{[2]}$: 10 x m\n- $W^{[1]}$: 10 x 10 (as $W^{[2]} A^{[1]} \\sim Z^{[2]}$)\n- $B^{[2]}$: 10 x 1\n\nBackprop\n\n- $dZ^{[2]}$: 10 x m ($~A^{[2]}$)\n- $dW^{[2]}$: 10 x 10\n- $dB^{[2]}$: 10 x 1\n- $dZ^{[1]}$: 10 x m ($~A^{[1]}$)\n- $dW^{[1]}$: 10 x 10\n- $dB^{[1]}$: 10 x 1","metadata":{}},{"cell_type":"code","source":"def init_params():\n    W1 = np.random.rand(10, 784) - 0.5\n    b1 = np.random.rand(10, 1) - 0.5\n    W2 = np.random.rand(10, 10) - 0.5\n    b2 = np.random.rand(10, 1) - 0.5\n    return W1, b1, W2, b2\n\n\ndef softmax(Z):\n    Z -= np.max(Z, axis=0)  # Subtract max value for numerical stability\n    A = np.exp(Z) / np.sum(np.exp(Z), axis=0)\n    return A\n\ndef forward_prop(W1, b1, W2, b2, X):\n    Z1 = W1.dot(X) + b1\n    A1 = ReLU(Z1)\n    Z2 = W2.dot(A1) + b2\n    A2 = softmax(Z2)\n    return Z1, A1, Z2, A2\n\n\ndef ReLU(Z):\n    return np.maximum(Z, 0)\n\ndef ReLU_deriv(Z):\n    return Z > 0\n\n#def ReLU(Z):\n#    return 1/(1+np.exp(-Z))\n\n#def ReLU_deriv(Z):\n#    return ReLU(Z)*(1-ReLU(Z))\n\ndef one_hot(Y):\n    one_hot_Y = np.zeros((Y.size, Y.max() + 1))\n    one_hot_Y[np.arange(Y.size), Y] = 1\n    one_hot_Y = one_hot_Y.T\n    return one_hot_Y\n\ndef backward_prop(Z1, A1, Z2, A2, W1, W2, X, Y):\n    one_hot_Y = one_hot(Y)\n    dZ2 = A2 - one_hot_Y\n    dW2 = 1 / m * dZ2.dot(A1.T)\n    db2 = 1 / m * np.sum(dZ2)\n    dZ1 = W2.T.dot(dZ2) * ReLU_deriv(Z1)\n    dW1 = 1 / m * dZ1.dot(X.T)\n    db1 = 1 / m * np.sum(dZ1)\n    return dW1, db1, dW2, db2\n\ndef update_params(W1, b1, W2, b2, dW1, db1, dW2, db2, alpha):\n    W1 = W1 - alpha * dW1\n    b1 = b1 - alpha * db1    \n    W2 = W2 - alpha * dW2  \n    b2 = b2 - alpha * db2    \n    return W1, b1, W2, b2\n\ndef show_2d(ar):\n    for i in ar:\n        print(\"\\n\")\n        for j in i:\n            print(j,end=\" \")","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:31:59.589446Z","iopub.execute_input":"2024-06-12T15:31:59.589902Z","iopub.status.idle":"2024-06-12T15:31:59.608752Z","shell.execute_reply.started":"2024-06-12T15:31:59.589854Z","shell.execute_reply":"2024-06-12T15:31:59.607488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_predictions(A2):\n    return np.argmax(A2, 0)\n\ndef get_accuracy(predictions, Y):\n    #print(predictions, Y)\n    return np.sum(predictions == Y) / Y.size\n\ndef gradient_descent(X, Y, alpha, iterations):\n    W1, b1, W2, b2 = init_params()\n    for i in range(iterations):\n        Z1, A1, Z2, A2 = forward_prop(W1, b1, W2, b2, X)\n        dW1, db1, dW2, db2 = backward_prop(Z1, A1, Z2, A2, W1, W2, X, Y)\n        W1, b1, W2, b2 = update_params(W1, b1, W2, b2, dW1, db1, dW2, db2, alpha)\n        if i % 100 == 0:\n            print(\"Iteration: \", i)\n            predictions = get_predictions(A2)\n            print(get_accuracy(predictions, Y))\n    return W1, b1, W2, b2\n\n\n","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:31:59.610569Z","iopub.execute_input":"2024-06-12T15:31:59.610952Z","iopub.status.idle":"2024-06-12T15:31:59.623249Z","shell.execute_reply.started":"2024-06-12T15:31:59.610887Z","shell.execute_reply":"2024-06-12T15:31:59.622046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"W1, b1, W2, b2 = gradient_descent(X_train, Y_train, 0.10, 1500)","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:31:59.625069Z","iopub.execute_input":"2024-06-12T15:31:59.625507Z","iopub.status.idle":"2024-06-12T15:34:08.446188Z","shell.execute_reply.started":"2024-06-12T15:31:59.625463Z","shell.execute_reply":"2024-06-12T15:34:08.445152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"~85% accuracy on training set.","metadata":{}},{"cell_type":"code","source":"def make_predictions(X, W1, b1, W2, b2):\n    _, _, _, A2 = forward_prop(W1, b1, W2, b2, X)\n    predictions = get_predictions(A2)\n    return predictions\n\ndef test_prediction(index, W1, b1, W2, b2):\n    current_image = X_dev[:, index, None]\n    prediction = make_predictions(X_dev[:, index, None], W1, b1, W2, b2)\n    label = Y_dev[index]\n    print(\"Prediction: \", prediction)\n    print(\"Label: \", label)","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:34:08.449371Z","iopub.execute_input":"2024-06-12T15:34:08.449849Z","iopub.status.idle":"2024-06-12T15:34:08.458083Z","shell.execute_reply.started":"2024-06-12T15:34:08.449801Z","shell.execute_reply":"2024-06-12T15:34:08.456921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"w1_copy = W1.copy()\nb1_copy = b1.copy()\nw2_copy = W2.copy()\nb2_copy = b2.copy()","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:34:08.459521Z","iopub.execute_input":"2024-06-12T15:34:08.459836Z","iopub.status.idle":"2024-06-12T15:34:08.470124Z","shell.execute_reply.started":"2024-06-12T15:34:08.459804Z","shell.execute_reply":"2024-06-12T15:34:08.469108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"normal\")\nplt.colorbar(mpl.cm.ScalarMappable(norm=mpl.colors.Normalize(-1, 1), cmap='RdBu_r'),\n             orientation='vertical', label='weight values')\nplt.pcolormesh( W2, vmin=-1., vmax=1., cmap='RdBu')\nplt.show()\nprint(\"fixed\")\nplt.colorbar(mpl.cm.ScalarMappable(norm=mpl.colors.Normalize(-1, 1), cmap='RdBu_r'),\n             orientation='vertical', label='weight values')\nplt.pcolormesh( W2fixed, vmin=-1., vmax=1., cmap='RdBu')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:34:08.471649Z","iopub.execute_input":"2024-06-12T15:34:08.472089Z","iopub.status.idle":"2024-06-12T15:34:08.957138Z","shell.execute_reply.started":"2024-06-12T15:34:08.472042Z","shell.execute_reply":"2024-06-12T15:34:08.955911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def best_matches(W,b,data):\n    a = []\n    for i, X in enumerate(data):\n        Z = W.dot(X) + b\n        A = ReLU(Z)\n        a.append((float(A),i))\n    return a\n\ndef look_at_hidden_node(ind,W1):\n    rank = best_matches(W1[ind],b1[ind],X_dev.T.copy())\n    srank = sorted(rank,key = lambda x: x[0], reverse = True)\n    arr = []\n    plt.gca().set_aspect('equal')\n    plt.pcolormesh( np.array(list(i*3 for i in W1[ind])).reshape((28,28)), vmin=-1., vmax=1., cmap='RdBu')\n    #plt.imshow(W1[ind].reshape((28,28)),cmap = \"bwr_r\")\n    plt.show()\n    for i in range(3):\n        temp = np.transpose(X_dev).copy()\n        temp = temp[srank[i][1]]\n        arr.append(temp)\n        for j in range(784):\n            arr[i][j] *= W1[ind][j]\n        plt.gca().set_aspect('equal')\n        plt.pcolormesh( arr[i].reshape((28,28)), vmin=-1., vmax=1., cmap='RdBu')\n        #plt.imshow( arr[i].reshape((28,28)), cmap='RdBu')\n        plt.show()\n\nlook_at_hidden_node(6,W1_blur_blur)","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:47:43.013862Z","iopub.execute_input":"2024-06-12T15:47:43.014360Z","iopub.status.idle":"2024-06-12T15:47:43.595400Z","shell.execute_reply.started":"2024-06-12T15:47:43.014312Z","shell.execute_reply":"2024-06-12T15:47:43.594390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, let's find the accuracy on the dev set:","metadata":{}},{"cell_type":"code","source":" ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"W2fixed = W2.copy()\nfor i in range(10):\n    for j in range(10):\n        if abs(W2fixed[i][j]) < 0.2: W2fixed[i][j] = 0\n\n    \ndev_predictions = make_predictions(X_dev, W1, b1, W2, b2)\nprint(f'normal W2 accuracy {get_accuracy(dev_predictions, Y_dev)*100}')\n\ndev_predictions = make_predictions(X_dev, W1, b1, W2fixed, b2)\nprint(f'fixed W2 accuracy {get_accuracy(dev_predictions, Y_dev)*100}')","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:34:09.266588Z","iopub.execute_input":"2024-06-12T15:34:09.267043Z","iopub.status.idle":"2024-06-12T15:34:09.280228Z","shell.execute_reply.started":"2024-06-12T15:34:09.266977Z","shell.execute_reply":"2024-06-12T15:34:09.278827Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def blur_weights(W1):\n    W1_blur = []\n    for ind in range(10):\n        te = W1[ind].reshape((28,28)).copy()\n        for i in range(28):\n            for j in range(28):\n                te[i][j] = te[i][j]/4 + te[(i+1)%28][j]/8 + te[i%28][(j+1)%28]/8 + te[(i-1)%28][j]/8 + te[i][(j-1)%28]/8 + te[(i+1)%28][(j-1)%28]/16 + te[(i-1)%28][(j-1)%28]/16+ te[(i-1)%28][(j+1)%28]/16+ te[(i+1)%28][(j+1)%28]/16\n        W1_blur.append(te.reshape((784)))\n    return np.array(W1_blur)\n\n# plt.imshow(te, cmap='RdBu')\n# plt.show()\nW1_blur = blur_weights(W1)\nW1_blur = np.array(W1_blur)\ndev_predictions = make_predictions(X_dev, W1, b1, W2, b2)\nprint(f'normal accuracy {get_accuracy(dev_predictions, Y_dev)*100}')\ndev_predictions = make_predictions(X_dev, W1_blur, b1, W2, b2)\nprint(f'blurred accuracy {get_accuracy(dev_predictions, Y_dev)*100}')\nW1_blur_blur = np.array(blur_weights(blur_weights(W1)))\ndev_predictions = make_predictions(X_dev, W1_blur_blur, b1, W2, b2)\nprint(f'blurred blurred accuracy {get_accuracy(dev_predictions, Y_dev)*100}')\n\n# for i in range(10):\n#     W1_blur_blur = blur_weights(W1_blur_blur)\n#     dev_predictions = make_predictions(X_dev, W1_blur_blur, b1, W2, b2)\n#     print(f'blurred {i+2} accuracy {get_accuracy(dev_predictions, Y_dev)*100}')\n\n    \n\nfor i in blur_weights(W1_blur_blur):\n    plt.imshow(i.reshape((28,28)), cmap = \"RdBu\")\n    plt.show()\n    ","metadata":{"execution":{"iopub.status.busy":"2024-06-12T16:05:51.088300Z","iopub.execute_input":"2024-06-12T16:05:51.088692Z","iopub.status.idle":"2024-06-12T16:05:52.990592Z","shell.execute_reply.started":"2024-06-12T16:05:51.088655Z","shell.execute_reply":"2024-06-12T16:05:52.989594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X_dev_COPY = X_dev.T.copy()\n\ninverted = []\nfor i in X_dev_COPY:\n    te = []\n    for j in i:\n        te.append(1-j)\n    inverted.append(te)\n\ndev_predictions = make_predictions(X_dev, W1, b1, W2fixed, b2)\nprint(f'dev set accuracy {get_accuracy(dev_predictions, Y_dev)*100}')\n\na = [list([0 for i in range(10)]) for i in range(10)]\nfor i in range(len(dev_predictions)):\n    a[dev_predictions[i]][Y_dev[i]] += 1\nplt.imshow(a, cmap = \"plasma\", origin = \"lower\")\ntab = []\nfor i in range(10):\n    s = 0\n    for j in range(10):\n        s += a[i][j]\n    print(f'kai speja {i}, acc {round(a[i][i]/s,2)}')\n    tab.append((i,round(a[i][i]/s,2)))\ntab.sort(key = lambda x: x[1], reverse = True)\nprint(tab)\n#show_2d(a)","metadata":{"execution":{"iopub.status.busy":"2024-06-12T15:34:11.036497Z","iopub.execute_input":"2024-06-12T15:34:11.036812Z","iopub.status.idle":"2024-06-12T15:34:11.668805Z","shell.execute_reply.started":"2024-06-12T15:34:11.036781Z","shell.execute_reply":"2024-06-12T15:34:11.667864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Still 84% accuracy, so our model generalized from the training data pretty well.","metadata":{}}]}