{"cells":[{"metadata":{"_uuid":"f6ae395b51319ea53b4689ab2e00bb9056987e44"},"cell_type":"markdown","source":"Hey everyone! I've seen that people have started using NNs for a while now, but I didn't see people using CNNs. Indeed, it is pretty unusual, at least to me, to use a kind of network that performs so well on image data to treat time series. However, I'm completely in love with this kind of NN and I started looking for ways to use it here and I found [this paper](https://arxiv.org/pdf/1710.00886.pdf). There are several issues that we have to deal with along the preprocessing and evaluation steps but I believe that this could be a good starting point for those who want to try using CNNs. I prepared this notebook as an introduction in the topic and necessary preprocessing in case someone wants to investigate it deeper. IMHO, the greatest advantage of this method - and other DL techniques - is that we don't need worry about feature extraction.\n\nFirst of all, we'll need to load the data. "},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport cv2 #deal with images\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom tqdm import tqdm\nnp.random.seed(42)\n\ntraining = pd.read_csv('../input/training_set.csv')\nmeta_training = pd.read_csv(\"../input/training_set_metadata.csv\")\nmerged = training.merge(meta_training, on = \"object_id\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"80bf41426357c30094c7b4b9932634dd3b1be95a"},"cell_type":"markdown","source":"The basic idea described in the paper is that we can transform a 1-dimensional time series into a 2-dimensional image by calculating its recurrences, that is, revealing in which points the series returns to a previous value.  This is done by using the equation:\n\n$$ R_{i,j} = \\theta(\\epsilon - || s_i - s_j ||)$$\n\nwhere a $s$ are the states (time series values), ||.|| is the norm  and $\\epsilon$ is a threshold distance and $\\theta$ is the Heaviside (step) function. \n\nBelow are the functions I used to create a dictionary of such R matrices, where each key is a different object. You'll notice that I preferred to use a slightly different approach by using a sigmoid rather than Heaviside step function. I did so because since we are dealing with several different objects, it turns out that each one has a different \"optimal\" $\\epsilon$ in order to create a proper image. So for simplicity I opted to use sigmoid because it doesn't need tunning for each image.\n\nAfter each R matrix for each object and each passband is made, there is another issue that we need to worry about.  Since each object has 6 passbands, we need a way to put all of those images - one for each passband -  into a single image. My solution was to stack them up , in some kind of resemblance to a RGB image that has 3 channels. Here we will have 6 channels instead. The problem is: it is not guaranteed that each passband has the same number of observations, resulting in images of different sizes. If we want to stack them, they need to have the same shape. To solve that, I cropped all of them to the minimum size possible for that object - and here is where we begin to lose information. For instance, if an object has R matrices of shapes (30x30), (50x50) and (48x48), I'd shape them all to (30x30) so they can be stacked. This is done by `crop_obj_plots` function."},{"metadata":{"trusted":true,"_uuid":"9275ed05e8576b5362be46c4425c1f252e392769","scrolled":false},"cell_type":"code","source":"###recurrent plot\n\ndef sigmoid(x):\n    '''\n    Returns the sigmoid of a value\n    '''\n    return 1/(1+np.exp(-x))\n\ndef R_matrix(signal, eps):\n    '''\n    Given a time series (signal) and an epsilon,\n    return the Recurrent Plot matrix\n    '''\n    R = np.zeros((signal.shape[0], signal.shape[0]))\n    for i in range(R.shape[0]):\n        for j in range(R.shape[1]):\n            R[i][j] = np.heaviside((eps - abs(signal[i] - signal[j])),1)\n    return R\n\n#using sigmoid rather than heaviside\n#because in this dataset the epsilon parameter needs to\n#change from object to object and therefore should be learned as well\ndef R_matrix_modified(signal):\n    '''\n    Given a time series (signal) and an epsilon,\n    return the modified Recurrent Plot matrix\n    using sigmoid rather than heaviside\n    '''\n    R = np.zeros((signal.shape[0], signal.shape[0]))\n    for i in range(R.shape[0]):\n        for j in range(R.shape[1]):\n            R[i][j] = sigmoid((abs(signal[i] - signal[j])))\n    return R\n\ndef create_objects_dict(merged_dataset):\n    '''\n    Input: dataset containing both training data and metadata\n    Creates a dictionary using each object as keys and\n    one R matrix for each passband in that object\n    '''\n    objects = {}\n    for obj in tqdm(np.unique(merged_dataset.object_id)):\n        R_passbands = []\n        for passband in np.unique(merged_dataset.passband):\n            obj_flux = merged_dataset[(merged_dataset.object_id == obj) & (merged_dataset.passband == passband)].flux.values\n            R_passbands.append(R_matrix_modified(obj_flux))\n        objects[obj] = (np.asarray(R_passbands), max(merged_dataset[merged_dataset.object_id == obj].target))\n    return objects\n\ndef get_minmax_shapes(obj_R_matrices):\n    '''\n    Given an R matrix, get the min and max width \n    to be used to crop and let all images from a given\n    object be of the same size so they can be concatenated\n    '''\n    min_length = 0\n    max_length = 0\n    for passband in np.unique(merged.passband):\n        if passband == 0:\n            length = len(obj_R_matrices[passband])\n            min_length = length\n            max_length = length\n        else:\n            length = len(obj_R_matrices[passband])\n            min_length = min(min_length, length)\n            max_length = max(max_length, length)\n    return (min_length, max_length)\n\ndef crop_obj_plots(objects):\n    '''\n    Accepts a dictionary where each key is a different object\n    and each value is a tuple - one slot with a list of R matrices and \n    the other with the target value (object class)\n    '''\n    for obj in tqdm(objects.keys()):\n        min_len, max_len = get_minmax_shapes(objects[obj][0])\n        for passband in np.unique(merged.passband):\n            objects[obj][0][passband] = objects[obj][0][passband][:min_len, :min_len]\n    return objects\n\nobjects = create_objects_dict(merged)\ncropped_objects = crop_obj_plots(objects)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b8ee58b710371d85ea845dbca519633d9aea9e7f"},"cell_type":"markdown","source":"Now each object has a 6-channel image. They all look like the image below:"},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"6c23ed7e28a7a2211e202c8d39b9f5f52e465376"},"cell_type":"code","source":"cropped_objects[730][0][3].shape\nplt.imshow(cropped_objects[730][0][0])","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"58b56fe7a6e14e1dad21956c3c1e9162fc3726b0"},"cell_type":"markdown","source":"Now let's the distribution of shapes across all images."},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"523cf3426033520cd9256f5e890098a113157d60"},"cell_type":"code","source":"from collections import Counter\n\nshapes = []\nfor key in tqdm(cropped_objects.keys()):\n    shapes.append(cropped_objects[key][0][0].shape[0])\nplt.hist(shapes, bins = 50)\nCounter(shapes)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"20e6ec82ae97d4677a37c1bb5940753d65d7e6cf"},"cell_type":"markdown","source":"Now we found the second problem and the key source of information loss: we have lots of different shapes, most around 10-12. The most complete images (longer observations) are just a small portion of the dataset. In order to properly work with the dataset, I will try to get most of them into the same shape. From my experience, the very small images are not very good predictors, so I'll just use the larger ones. Now I'll put all images with shapes (50,50) or larger  into shape (57,57). "},{"metadata":{"trusted":true,"_uuid":"22aab92ed0c8d11faa24b664b894e4715655dacd"},"cell_type":"code","source":"import math\ncropped_2 = np.copy(cropped_objects).item()\nfor key in tqdm(cropped_2.keys()):\n    shape = cropped_2[key][0][0].shape[0]\n    if shape < 11:\n        for passband in np.unique(merged.passband):\n            #how much we will increase the border\n            increaseBorder = abs(shape-11)/2\n            cropped_2[key][0][passband] = cv2.copyMakeBorder(src = cropped_2[key][0][passband],\n                                                             top = math.ceil(increaseBorder), \n                                                             left = math.ceil(increaseBorder),\n                                                             bottom = round(increaseBorder),\n                                                             right = round(increaseBorder),\n                                                             borderType = cv2.BORDER_REFLECT)\n    elif shape>11 and shape < 25:\n        for passband in np.unique(merged.passband):\n            cropped_2[key][0][passband] = cropped_2[key][0][passband][:-(shape-11), :-(shape-11)]\n            \n    \n    elif shape >= 50 and shape < 57:\n        for passband in np.unique(merged.passband):\n            increaseBorder57 = abs(shape-57)/2\n            cropped_2[key][0][passband] = cv2.copyMakeBorder(src = cropped_2[key][0][passband],\n                                                             top = math.ceil(increaseBorder57), \n                                                             left = math.ceil(increaseBorder57),\n                                                             bottom = round(increaseBorder57),\n                                                             right = round(increaseBorder57),\n                                                             borderType = cv2.BORDER_REFLECT)\n    else:\n        continue","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"320469dd93d13fb806e7055fee6fa1c4d73dc1ee"},"cell_type":"markdown","source":"Now let's check how the shapes look like."},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"e069117d6f3bb7363bd2c1534a413b5913aad388"},"cell_type":"code","source":"shapes = []\nfor key in tqdm(cropped_2.keys()):\n    shapes.append(cropped_2[key][0][0].shape[0])\nplt.hist(shapes, bins = 50)\nCounter(shapes)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"92b86d1a47cb0aad827b53dc7bc61d0849a580bb"},"cell_type":"markdown","source":"Now we can stack all the images with the same shape to form our input."},{"metadata":{"trusted":true,"_uuid":"7e51c39a5426602ea2bff88332f94c6699999945","scrolled":true},"cell_type":"code","source":"objects = list(cropped_2.keys())\ninput_images = list()\nlabels = list()\nfor key in tqdm(objects):\n    if cropped_2[key][0][0].shape[0] == 57:\n        img = np.stack((cropped_2[key][0][0],\n                        cropped_2[key][0][1],\n                        cropped_2[key][0][2],\n                        cropped_2[key][0][3],\n                        cropped_2[key][0][4]), axis = -1)  \n                                           \n        input_images.append(np.expand_dims(img, axis = 0))\n        labels.append(cropped_2[key][1])                                                        \ninput_images = np.vstack(input_images)     \ninput_images.shape","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0f863d607075c4cc9d2ef4a46068e52e165e5257"},"cell_type":"markdown","source":"We have lots of images with the same shape and therefore we can use them in our CNN. Here I'll just use images with shape 57x57x5.  \n\nThe following cell is to split between train, validation and test samples and also binarize the labels. "},{"metadata":{"trusted":true,"_uuid":"aaf22c66ead98a0f7a20c51c4e5c79651102eb6b"},"cell_type":"code","source":"#LabelBinarizer and train-test split\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelBinarizer\nfrom keras.utils import np_utils\n\n\ntrain_fraction = 0.8\n\nencoder = LabelBinarizer()\ny = encoder.fit_transform(labels)\nx = input_images\n\ntrain_tensors, test_tensors, train_targets, test_targets =\\\n    train_test_split(x, y, train_size = train_fraction, random_state = 42)\n\nval_size = int(0.5*len(test_tensors))\n\nval_tensors = test_tensors[:val_size]\nval_targets = test_targets[:val_size]\ntest_tensors = test_tensors[val_size:]\ntest_targets = test_targets[val_size:]\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a8157e13d43eceb23d5618359a89493502ccefc1"},"cell_type":"markdown","source":"We then train the model. The input_shape needs to be (None, None, 5) rather than (57, 57, 5) - that we know will be our training input shapes - because we want to test the model in images of different sizes as well (remember the huge amount of images of size 11x11?). If we don't set input_shape's width and height to None, our convolutional layers will throw an error. "},{"metadata":{"trusted":true,"_uuid":"9b5bfe7796eadd4e97dc51b9d1ab86d687de421b","scrolled":false},"cell_type":"code","source":"from keras.layers import Conv2D, MaxPooling2D, GlobalAveragePooling2D, GlobalMaxPooling2D\nfrom keras.layers import Dropout, Flatten, Dense, LeakyReLU\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom keras.models import Sequential\nfrom tensorflow import set_random_seed\n\nset_random_seed(42)\n\nearly_stopping = EarlyStopping(monitor = 'val_loss', patience = 10)\ncheckpointer = ModelCheckpoint(filepath='weights.hdf5', \n                               verbose=1, save_best_only=True)\n\n\nmodel = Sequential()\nmodel.add(Conv2D(filters = 128, kernel_size = (3,3), padding = 'same', activation = 'elu', input_shape = (None, None,5)))\nmodel.add(Conv2D(filters = 128, kernel_size = (3,3), padding = 'same', activation = 'elu'))\nmodel.add(Dropout(0.3))\nmodel.add(MaxPooling2D(pool_size = 2)) \n\nmodel.add(Dense(256, activation = 'relu'))\nmodel.add(Dense(5, activation = 'softmax'))\n\n\nmodel.compile(optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'])\nepochs = 10\nmodel.fit(train_tensors, train_targets, \n          validation_data=(val_tensors, val_targets),\n          epochs=epochs, batch_size=128, verbose=1, callbacks = [early_stopping, checkpointer])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"a366577a1c8de9820517a7d3c67e94c264f267d6","scrolled":false},"cell_type":"code","source":"model.load_weights('weights.hdf5')\n\ncell_predictions =  [np.argmax(model.predict(np.expand_dims(tensor, axis=0))) for tensor in test_tensors]\n\ntest_accuracy = 100*np.sum(np.array(cell_predictions)==np.argmax(test_targets, axis=1))/len(cell_predictions)\nprint('Test accuracy: %.4f%%' % test_accuracy)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"155d5f3c15d1aec65fc71f015d97ba65def6aea1"},"cell_type":"markdown","source":"Note that I'm not considering any weights for any particular class here. The loss function that gets the best results so far is categorical crossentropy. By doing so, I'm getting an accuracy of around 60% in the classification of the larger images (objects that have observations for a longer period of time). Making the predictions on the images of shape (11x11), we got:"},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"cfdbe5f0e0bd51ad1884e00f13d57323d05f0de8"},"cell_type":"code","source":"#manipulations to get our images into the proper dimensions (done before with shape 57)\nobjects = list(cropped_2.keys())\nsmall_labels = list()\nsmall_input_images = list()\nfor key in tqdm(objects):\n    if cropped_2[key][0][0].shape[0] == 11:\n        img = np.stack((cropped_2[key][0][0],\n                        cropped_2[key][0][1],\n                        cropped_2[key][0][2],\n                        cropped_2[key][0][3],\n                        cropped_2[key][0][4]), axis = -1)  \n                                           \n        small_input_images.append(np.expand_dims(img, axis = 0))\n        small_labels.append(cropped_2[key][1])                                                        \nsmall_input_images = np.vstack(small_input_images)     \nsmall_input_images.shape\n\n#predictions\ncell_predictions =  [np.argmax(model.predict(np.expand_dims(small_image, axis=0))) for small_image in small_input_images]\n\nsmall_encoder = LabelBinarizer()\nsmall_labels_encoded = encoder.fit_transform(small_labels)\ntest_accuracy = 100*np.sum(np.array(cell_predictions)==np.argmax(small_labels_encoded, axis=1))/len(cell_predictions)\nprint('Accuracy in small images: %.4f%%' % test_accuracy)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"f1d9d1072153b028711717c748c72bd9810b32d8"},"cell_type":"markdown","source":"The model's accuracy for a much smaller test set is awfully low, around 13%-15%. There are at least two possible ways out of this problem (I'm not sure if any would actually work. Just some hypothesis here):\n- Train the model using just the smaller images and try to generalize to the larger ones\n- Train the model using both large and small images\n\nI'll leave this open for more experienced data scientists to discuss if they want to :) Meanwhile I'll keep working on how to solve this generalization issue and update here if I have an Eureka moment."}],"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":4,"nbformat_minor":4}