{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Shameless Plug\nI recently published a free course on Neural Networks, on YouTube, called [Neural Networks For Your Dog](https://www.youtube.com/playlist?list=PL9oKUrtC4VP5N3VtTTjhTfiHoFXmnrgPW) - So easy your dog could learn them ;)\n\nIn my course, I derive and implement the neural network model used below. So, if you're wondering how I came up with the below code, it's all explained in my videos. \n\nI could not have put this together without the help of Kaggle and its awesome community. So, thanks!","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the data\ntrain = pd.read_csv(\"/kaggle/input/digit-recognizer/train.csv\")\ntest = pd.read_csv(\"/kaggle/input/digit-recognizer/test.csv\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Inspect the data\nprint(f'train.shape: {train.shape}')\nprint(f'test.shape: {test.shape}')\n\ntrain.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Split train into (85% train_train, 15% train_validation)\ngen = np.random.default_rng(seed = 1234)\nval_ids, train_ids = np.split(gen.choice(len(train), len(train), replace=False), [int(0.15 * len(train))])\ntrain_train = train.iloc[train_ids]\ntrain_val = train.iloc[val_ids]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Helper Functions","metadata":{}},{"cell_type":"code","source":"def softmax(x):\n    \"\"\"\n    Calculate row-wise softmax\n    :param x: 2d array where (i,j) gives the jth input value for the ith sample\n    :return: 2d array with the same shape as the input, with softmax applied to each row-vector\n             As a result, the elements in each row can be interpretted as probabilities that sum to one\n    \"\"\"\n\n    return np.exp(x)/np.sum(np.exp(x), axis=1)[:, None]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def one_hot(x):\n    \"\"\"\n    One-hot-encode an array\n    :param x: 1d array where element (i) gives the true label for sample i\n    :return: tuple of (onehot, classes) where:\n             - onehot is a NxK array where N = len(x), K = len(np.unique(x)) and\n               element (i,j) = 1 if string_arr[i] == np.unique(x)[j], 0 otherwise\n             - classes is a 1d array of classes corresponding to the columns of onehot\n    \"\"\"\n\n    classes, inverse = np.unique(x, return_inverse=True)\n    onehot = np.eye(classes.shape[0], dtype='int64')[inverse]\n    return (onehot, classes)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def cross_entropy(Yhat, Y):\n    \"\"\"\n    Calculate row-wise cross entropy\n    :param Yhat: NxK array where (i,j) gives the predicted probability of class j for sample i\n    :param Y: either:\n              1) NxK array where (i,j) gives the true probability of class j for sample i or\n              2) a 1-D array where element i gives the index of the true class for sample i\n    :return: 1-D array with N elements, where element i gives the cross entropy for the ith sample\n    \"\"\"\n\n    if Y.ndim == 1:\n        ce = -np.log(Yhat[np.arange(len(Y)), Y])\n    else:\n        ce = -np.sum(Y * np.log(Yhat), axis=1)\n\n    return ce","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def logistic(x):\n    \"\"\"\n    standard logistic function\n    uses the identity: 1/(1 + e^(-x)) = e^x/(e^x + 1)\n    to prevent double precision issues when x is a big negative number\n    :param x: numpy array\n    :return: 1/(1 + e^(-x))\n    \"\"\"\n\n    mask = x > 0\n    y = np.full(shape=x.shape, fill_value=np.nan)\n    y[mask] = 1 / (1 + np.exp(-(x[mask])))\n    y[~mask] = np.exp(x[~mask]) / (np.exp(x[~mask]) + 1)\n    return y","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Neural Network","metadata":{}},{"cell_type":"code","source":"class NNet():\n    \"\"\"\n    NNet with stochastic gradient descent and validation loss monitoring\n    \"\"\"\n\n    def __init__(self, Ws=None, y_classes=None):\n        \"\"\"\n        Initialization\n        :param Ws: optional list of weight matrices (list of 2-D numpy arrays)\n        :param y_classes: optional array of y_classes (1-D numpy array with >= 2 elements)\n        \"\"\"\n\n        self.Ws = Ws\n        self.y_classes = y_classes\n\n    def fit(self, X, y, hiddenNodes, Xval=None, yval=None, stepSize=0.01, ITERS=100, batchSize=None, seed=None):\n        \"\"\"\n        Find the best weights via stochastic gradient descent\n        :param X: training features\n        :param y: training labels. 1-d array with >= 2 classes\n        :param hiddenNodes: list indicating how many nodes to use in each hidden layer, excluding bias nodes\n        :param Xval: optional validation features\n        :param yval: optional validation labels. 1-d array with >= 2 classes\n        :param stepSize: AKA \"learning rate\" AKA \"alpha\" used in gradient descent\n        :param ITERS: How many gradient descent steps to make?\n        :param batchSize: How many samples to user per batch? If None, use all samples\n        :return: None. Update self.y_classes, self.W1, self.W2\n        \"\"\"\n\n        # Validate X dimensionality\n        if X.ndim != 2:\n            raise AssertionError(f\"X should have 2 dimensions but it has {X.ndim}\")\n\n        # Validate Ws type\n        if not isinstance(hiddenNodes, list):\n            AssertionError(\"hiddenNodes should be a list of integers\")\n\n        # Determine unique y classes\n        y01, y_classes = one_hot(y)\n        if len(y_classes) < 2:\n            AssertionError(f\"y should have at least 2 distinct classes, but instead it has {len(y_classes)}\")\n\n        if yval is not None:\n            y01_val, yval_classes = one_hot(yval)\n\n        # Initialization (note Ws is a list of weight matrices)\n        gen = np.random.default_rng(seed)\n        X1 = np.insert(X / 255, obj=X.shape[1], values=1, axis=1)\n        Ws = [None] * (len(hiddenNodes) + 1)\n        Ws[0] = gen.uniform(low=-1, high=1, size=(X1.shape[1], hiddenNodes[0]))\n        for i in range(1, len(hiddenNodes)):\n            Ws[i] = gen.uniform(low=-1, high=1, size=(hiddenNodes[i - 1] + 1, hiddenNodes[i]))\n        Ws[i + 1] = gen.uniform(low=-1, high=1, size=(hiddenNodes[i] + 1, len(y_classes)))\n\n        # Initialize lists to store Xs, Zs, and gradients\n        Zs = [None] * len(Ws)\n        Xs = [None] * len(Ws)\n        gradWs = [None] * len(Ws)\n\n        # Determine number of batches\n        if batchSize is None:\n            Nbatches = 1\n        else:\n            Nbatches = np.ceil(X1.shape[0]/batchSize).astype('int64')\n\n        # Initialize lists to store performance stats\n        CE_train_list = []\n        Accuracy_train_list = []\n        CE_validation_list = []\n        Accuracy_validation_list = []\n\n        # Train\n        for i in range(ITERS):\n\n            # mini batches\n            idxs = gen.choice(X1.shape[0], size=X1.shape[0], replace=False)\n            batches = np.array_split(idxs, Nbatches)\n\n            # Loop over batches\n            for b in range(Nbatches):\n                batch_idxs = batches[b]\n                Xs[0] = X1[batch_idxs]\n\n                # Make predictions (forward pass)\n                for j in range(len(Ws)):\n                    Zs[j] = Xs[j] @ Ws[j]\n                    if j + 1 < len(Xs):\n                        Xs[j + 1] = np.insert(logistic(Zs[j]), obj=Zs[j].shape[1], values=1, axis=1)\n                yhat_probs = softmax(Zs[-1])\n\n                if b == Nbatches - 1:\n                    # Calculate training cross entropy loss, accuracy\n                    yhat_probs_train = self.predict(X, type='probs', Ws=Ws, y_classes=y_classes)\n                    yhat_classes_train = y_classes[np.argmax(yhat_probs_train, axis=1)]\n                    ce_train = cross_entropy(yhat_probs_train, y01)\n                    CE_train = np.mean(ce_train)\n                    accuracy_train = np.mean(yhat_classes_train == y)\n\n                    CE_train_list.append(CE_train)\n                    Accuracy_train_list.append(accuracy_train)\n\n                    if Xval is None or yval is None:\n                        print(f'iteration: {i}, train cross entropy loss: {CE_train}, train accuracy: {accuracy_train}')\n                    else:\n                        yhat_probs_val = self.predict(Xval, type='probs', Ws = Ws, y_classes = y_classes)\n                        yhat_classes_val = y_classes[np.argmax(yhat_probs_val, axis=1)]\n                        ce_val = cross_entropy(yhat_probs_val, y01_val)\n                        CE_val = np.mean(ce_val)\n                        accuracy_val = np.mean(yhat_classes_val == yval)\n\n                        CE_validation_list.append(CE_val)\n                        Accuracy_validation_list.append(accuracy_val)\n\n                        print(f'iter: {i}, '\n                              f'train CE: {np.round(CE_train, 3)}, val CE: {np.round(CE_val, 3)}, '\n                              f'train acc: {np.round(accuracy_train, 3)}, val acc: {np.round(accuracy_val, 3)}')\n\n\n                # Calculate gradients (backward pass)\n                gradZ = (yhat_probs - y01[batch_idxs])[:, None, :]\n                for j in range(len(Ws) - 1, -1, -1):\n                    gradWs[j] = np.transpose(Xs[j][:, None, :], axes=[0, 2, 1]) @ gradZ\n                    gradWs[j] = gradWs[j].mean(axis=0)\n                    gradX = (gradZ @ np.transpose(Ws[j]))[:, :, :-1]\n                    gradZ = gradX * (Xs[j] * (1 - Xs[j]))[:, None, :-1]\n\n                # Update weights (gradient step)\n                for j in range(len(Ws)):\n                    Ws[j] -= gradWs[j] * stepSize\n\n        # Update class vars\n        self.y_classes = y_classes\n        self.Ws = Ws\n        self.CE_train_list = CE_train_list\n        self.Accuracy_train_list = Accuracy_train_list\n        self.CE_validation_list = CE_validation_list\n        self.Accuracy_validation_list = Accuracy_validation_list\n\n    def predict(self, X, type='classes', Ws = None, y_classes = None):\n        \"\"\"\n        Predict on X\n        :param X: 2-D array with >= 1 column of real-valued features\n        :param type: If 'classes', predicted classes, else if 'probs', predicted class probabilities\n        :param Ws: list of 2-D arrays (weight matrices). If None, use self.Ws\n        :param y_classes: numpy array of y classes. If None, use self.y_classes\n        :return: if type = 'probs' then probabilities else if type = 'classes' then classes\n        \"\"\"\n\n        if Ws is None:\n            Ws = self.Ws\n        if y_classes is None:\n            y_classes = self.y_classes\n\n        if Ws is None:\n            raise AssertionError(f\"Need to fit() before predict()\")\n        if X.ndim != 2:\n            raise AssertionError(f\"X should have 2 dimensions but it has {X.ndim}\")\n        if X.shape[1] != len(Ws[0]) - 1:\n            raise AssertionError(\n                f\"Perceptron was fit on X with {len(Ws[0]) - 1} columns but this X has {X.shape[1]} columns\")\n\n        # Make predictions (forward pass)\n        X1 = np.insert(X / 255, obj=X.shape[1], values=1, axis=1)\n        for j in range(len(Ws)):\n            Z = X1 @ Ws[j]\n            if j < len(Ws) - 1:\n                X1 = np.insert(logistic(Z), obj=Z.shape[1], values=1, axis=1)\n        yhat_probs = softmax(Z)\n\n        if type == 'probs':\n            return yhat_probs\n        elif type == 'classes':\n            yhat_classes = y_classes[np.argmax(yhat_probs, axis=1)]\n            return yhat_classes","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train the model","metadata":{}},{"cell_type":"code","source":"# Initialize & fit neural network\nnn = NNet()\nnn.fit(\n    X = train_train.drop(columns='label').to_numpy(),\n    y = train_train.label.to_numpy(),\n    Xval = train_val.drop(columns='label').to_numpy(),\n    yval = train_val.label.to_numpy(),\n    hiddenNodes = [100, 50, 20],\n    stepSize = 0.1,\n    batchSize = 100,\n    ITERS = 250,\n    seed = 0\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Identify the best number of epochs (where the validation loss is lowest)\nbest_epochs = np.argmax(nn.Accuracy_validation_list)\nprint(f'best number of epochs to run: {best_epochs}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Refit the neural network on the full training data with this many epochs\nnn2 = NNet()\nnn2.fit(\n    X = train.drop(columns='label').to_numpy(),\n    y = train.label.to_numpy(),\n    hiddenNodes = [100, 50, 20],\n    stepSize = 0.1,\n    batchSize = 100,\n    ITERS = best_epochs,\n    seed = 0\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Make Predictions","metadata":{}},{"cell_type":"code","source":"# Build submissions dataframe\npred_labels = nn2.predict(test.to_numpy())\nsubmission = pd.DataFrame({\n    'ImageId': range(1,28001),\n    'Label': pred_labels\n})\nsubmission.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Save submissions to CSV\nsubmission.to_csv(\"nnets_for_your_dog.csv\", index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}