{"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":"The idea of this kernel is to extract embedded features from the `test_set` and `training_set` by using autoencoders. Here I've implemented autoencoders by using convolutional layers, but with the same idea it can be used with LSTM or it can be also extended to Variational Autoencoders. The nice property of the autoencoders is that it reduces the dimensionality by removing the noise. At that stage I am not certain doe, how well the implemented autoencoding scheme is performing in terms of noise removal and preserving the information from the time series, but it will be tested soon.  :)\n\nThe extracter features then can be used as an input into some machiene learning forecaster. Here for the first test I am using random forest, which is fast, but in our unbalanced labels case it is not the best choice.","metadata":{"_uuid":"e0062ce149146256a8e88555824f351c5e94309f"}},{"cell_type":"markdown","source":"Helpers for data pre-processing","metadata":{"_uuid":"4881a3a281c1c49f6513737369db75852c625688"}},{"cell_type":"code","source":"import codecs\nimport csv\nimport bisect\nimport sys\nimport numpy as np\nimport gc\n\nFIRST_TIME_SERIE_DATE = 59580 # Wednesday November 17, 1858 is the zero date\nPB_COUNT = 6\nVAR_COUNT = 5\nVAR_META_COUNT = 10\nTIME_SERIES_DIM = 1095 # 3 years of history(3 * 365)\nTEST_N = 3490000\nTRAIN_N = 7848\nTRAIN_PATH = \"../input/training_set.csv\"\nTEST_PATH = \"../input/test_set.csv\"\nTRAIN_METADATA_PATH = \"../input/training_set_metadata.csv\"\nTEST_METADATA_PATH = \"../input/test_set_metadata.csv\"\n\nENCODING = 'utf-8'\n\ndef rescale_time_serie(time_serie):\n    '''\n    Standardize the time-serie for flux and flux error for each passband and object.\n    ''' \n    '''\n    scaling_time_serie = time_serie[:,0:2,:,:]\n    mask = time_serie[:,4:5,:,:]\n    aggregation_axis = 2\n    n = mask.sum(axis=aggregation_axis, keepdims=True)\n    mean = scaling_time_serie.sum(axis=aggregation_axis, keepdims=True) / n\n    std = np.sum((scaling_time_serie - mean) * (scaling_time_serie - mean),\\\n                                  axis=aggregation_axis, keepdims=True) / n\n    std = np.sqrt(std)\n    \n    std[std==0] = 1\n    \n    return ((scaling_time_serie - mean) / std, mean, std)\n    '''\n    \n    scaling_time_serie = time_serie[:,0:2,:,:]\n    mask = time_serie[:,4:5,:,:]\n    aggregation_axis = 2\n    max = scaling_time_serie.max(axis=aggregation_axis, keepdims=True)\n    min = scaling_time_serie.max(axis=aggregation_axis, keepdims=True)\n    \n    diff = min - max\n    diff[diff==0] = 1\n    \n    return ((scaling_time_serie - min) / diff, min, diff)\n    \n    \ndef parse_line_to_timeserie(time_serie,line):\n    \n    mjd = float(line[1])\n    day = int(mjd)\n    hour = mjd % 1\n    \n    passband = int(line[2])\n    flux = float(line[3])\n    flux_err = float(line[4])\n    detected = int(line[5])\n    \n    time_serie[0,day - FIRST_TIME_SERIE_DATE,passband] = flux\n    time_serie[1,day - FIRST_TIME_SERIE_DATE,passband] = flux_err\n    time_serie[2,day - FIRST_TIME_SERIE_DATE,passband] = detected\n    time_serie[3,day - FIRST_TIME_SERIE_DATE,passband] = hour\n    time_serie[4,day - FIRST_TIME_SERIE_DATE,passband] = 1\n    \n    return time_serie\n\ndef create_time_serie():\n    # TODO: sparse np format?\n    # One extra for the mask\n    time_serie = np.zeros((VAR_COUNT,TIME_SERIES_DIM,PB_COUNT))\n    return time_serie\n\ndef initialize_time_serie(line):\n    time_serie = create_time_serie()\n    return parse_line_to_timeserie(time_serie,line)\n\ndef iter_time_serie_by_id(file, chunk_size = 10000,max_n=sys.maxsize):\n    \"\"\"\n    The method loops over the csv file and yields\n    all the time series in a chunks with the specified\n    chunk size. The function assumes that the lines of\n    the csv files are grouped by object id. The chunk\n    size represents the count of time series within a chunk.\n    The time serie is the collection of all chronologically\n    ordered observations for an object id and it consists of\n    the dimensions for each feature value of interest.\n    \"\"\"\n    reader = csv.reader(file)\n    chunk = list()\n    ids = list()\n    headers = next(reader)\n    first_line = next(reader)\n    time_serie = initialize_time_serie(first_line)\n    id = first_line[0]\n\n    counter = 0\n    for l in reader:\n        if id == l[0]:\n            id = l[0]\n            time_serie = parse_line_to_timeserie(time_serie,l)\n            continue\n\n        chunk.append(time_serie)\n        ids.append(id)\n        counter += chunk_size\n        if counter > max_n:\n            break\n        if len(chunk) >= chunk_size:\n            yield (np.array(chunk), ids)\n            chunk = list()\n            ids = list()\n            \n        id = l[0]\n        time_serie = initialize_time_serie(l)\n        \n    if counter < max_n:\n        yield (np.array(chunk), ids)\n    \ndef write_input(file, max_n=sys.maxsize): \n    '''\n    The function assumes that the lines of\n    the csv files are grouped by object id.\n    '''\n    flux_list = list()\n    id_to_idx = dict()\n\n    reader = csv.reader(file)\n    headers = next(reader)\n\n    for l in reader:\n        object_id = l[0]\n        if object_id in id_to_idx:\n            flux_list[id_to_idx[object_id]] = parse_line_to_timeserie(time_serie,l)\n        else:\n           # TODO: we have to use some shuffeling insted of taking the first n samples\n            if len(flux_list) > max_n:\n                break\n\n            id_to_idx[object_id] = len(flux_list)\n            time_serie = initialize_time_serie(l)      \n            flux_list.append(time_serie)\n                \n    return (np.array(flux_list), id_to_idx)\n        ","metadata":{"_uuid":"da8b131b5c511bc44ef003301fca0663f3fff8ef","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We prepare the times-series input for autoencoding.  Ideally we have to train on the whole data(`test` + `train`).  That is done later by using the `fit_generator`.","metadata":{"_uuid":"0a7688057f86a299238c24913dd230c465eefb5e"}},{"cell_type":"code","source":"\nVAL_COUNT = int(TRAIN_N * 0.4)\n'''with codecs.open(TRAIN_PATH, \"r\", ENCODING) as fp: \n    train_ts, train_id_to_idx = write_input(fp)\n\nval_ts = train_ts[(N-VAL_COUNT):(N-1),]\n\n(val_scaled, val_add_factor, val_mult_factor) = rescale_time_serie(val_ts)\n(train_scaled, train_add_factor, train_mult_factor) = rescale_time_serie(train_ts)\n'''","metadata":{"_uuid":"f55a2b8b6cb2418268d9055d76003eb7e9ac8288","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The design of the autoencoder is composed by convolutional layers extracting the weekly, monthly and yearly information from the time series. We distinguish two types of input for the Convolutional Layers. One is with input the the values of interest which we want to reconstruct: `flux` and `flux_error` and the extra variables from the input which are just for helping reconstructing better `non-censored`, `detected` and `hours`. We have two parallel convolutions for both group of inputs, but it is good idea also input them stacked into one convolution since they are all scaled between 0 and 1. We have extra two inputs namely `add_factor` and `mult_factor` which are basically helping the network to reconstruct the real scale of the time-series . In the last output layer we scale-up the values by multiplying by the `mult_factor` and adding `add_factor` so that the output can be compared in the loss function to the real time-series values. We are learning parameters on the scaled-down values first so that we improve the numerical stability of the learning part. As a loss function we use the `mse` criterion. We put the same weights for each one of the two losses(`flux_loss` and ` flux_err_loss`). By computing the losses it is crucial to mask also the input, because for each time-serie we have considerable amount of the values censored. We multiply by zero the errors from these censored spots so that they don't end up into the loss function.","metadata":{"_uuid":"0cd9425f50f23a0b0a4c95932fb26283d8045a3f"}},{"cell_type":"code","source":"from keras.layers import Input, Dense, Dropout, Convolution2D, Reshape, Flatten,\\\n                         concatenate, BatchNormalization, Multiply, Add\nfrom keras.models import Model\nfrom keras import backend as K, objectives\nfrom keras.optimizers import Adam\nfrom keras import regularizers\n\nENCODING_DIM = 32\nDROPOUT_RATE = 0.03\nLR = 0.0001\n\nFACTOR_INPUT_SHAPE = (2, 1, PB_COUNT)\n\n# The input shape of the scaled time-series\nSCALED_INPUT_SHAPE = (2, TIME_SERIES_DIM, PB_COUNT)\n\n# The output shape of the scaled time-series\nSCALED_OUTPUT_DIM = 2 * TIME_SERIES_DIM * PB_COUNT\nSCALED_OUTPUT_SHAPE = (2, TIME_SERIES_DIM, PB_COUNT)\n\n# The shape of the inputs between 0 and 1\nBIN_INPUT_SHAPE = (3, TIME_SERIES_DIM, PB_COUNT)\n\ndef tailored_loss(true_output, loss_input):\n    flux_loss = objectives.mean_squared_error(loss_input[:,0,:,:] * true_output[:,2,:,:], true_output[:,0,:,:])\n    flux_err_loss = objectives.mean_squared_error(loss_input[:,1,:,:] * true_output[:,2,:,:], true_output[:,1,:,:])\n\n    return 0.5 * flux_loss + 0.5 * flux_err_loss\n\nscaled_ts_input = Input(shape = SCALED_INPUT_SHAPE)\nbin_ts_input = Input(shape = BIN_INPUT_SHAPE)\nadd_factor = Input(shape = FACTOR_INPUT_SHAPE)\nmult_factor = Input(shape = FACTOR_INPUT_SHAPE)\n\nxf = concatenate([add_factor, mult_factor])\nxf = Dense(2 * ENCODING_DIM, activation='relu')(xf)\nxf = Dropout(DROPOUT_RATE)(xf)\nxf = Flatten()(xf)\n\n# Weekly pattern convolution\nxs = Convolution2D(32, (7, 6), strides=(7, 6), activation='relu', data_format=\"channels_first\")(scaled_ts_input)\nxs = BatchNormalization()(xs)\nxb = Convolution2D(32, (7, 6), strides=(7, 6), activation='relu', data_format=\"channels_first\")(bin_ts_input)\nprint(\"Output dimensions of the weekly pattern convolution\", xs.shape)\n\n# Monthly pattern convolution\nxs = Convolution2D(16, (4, 1), strides=(4, 1), activation='relu', data_format=\"channels_first\")(xs)\nxs = BatchNormalization()(xs)\nxb = Convolution2D(16, (4, 1), strides=(4, 1), activation='relu', data_format=\"channels_first\")(xb)\nprint(\"Output dimensions of the monthly pattern convolution\", xs.shape)\n\n# Annual pattern convolution\nxs = Convolution2D(16, (12, 1), strides=(12, 1), activation='relu', data_format=\"channels_first\")(xs)\nxs = BatchNormalization()(xs)\nxb = Convolution2D(16, (12, 1), strides=(12, 1), activation='relu', data_format=\"channels_first\")(xb)\nprint(\"Output dimensions of the annual pattern convolution\", xs.shape)\n\nx = concatenate([xs, xb])\nx = Dense(4 * ENCODING_DIM, activation='relu')(x)\nx = Dropout(DROPOUT_RATE)(x)\nx = Flatten()(x)\nx = concatenate([x, xf])\n\nprint(\"Input dimensions of the embedding layer.\", x.shape)\nencoded = Dense(4 * ENCODING_DIM, activation='relu')(x)\nencoded = Dropout(DROPOUT_RATE)(encoded)\nencoded = Dense(2 * ENCODING_DIM, activation='relu')(encoded)\nencoded = Dropout(DROPOUT_RATE)(encoded)\nencoded = Dense(ENCODING_DIM, activation='relu')(encoded)\nprint(\"Output dimensions of the embedding layer.\", encoded.shape)\n\nx_mid = Dropout(DROPOUT_RATE)(encoded)\nx_mid = Dense(2 * ENCODING_DIM, activation='relu')(x_mid)\nx_mid = Dropout(DROPOUT_RATE)(x_mid)\nx_mid = Dense(4 * ENCODING_DIM, activation='relu')(x_mid)\nx_mid = Dropout(DROPOUT_RATE)(x_mid)\n\ndecoded = Dense(SCALED_OUTPUT_DIM, activation='sigmoid')(x_mid)\ndecoded = Reshape(SCALED_OUTPUT_SHAPE, input_shape=(SCALED_OUTPUT_DIM,))(decoded)\n\ndecoded = Multiply()([decoded, mult_factor]) # Scaling up in the output layer\ndecoded = Add()([decoded, add_factor]) # Scaling up in the output layer\n\nprint(\"Decoded input dimensions.\", decoded.shape)\n\nautoencoder = Model([scaled_ts_input, bin_ts_input, add_factor, mult_factor], decoded)\n\noptimizer = Adam(lr = LR)\nautoencoder.compile(optimizer = optimizer, loss = tailored_loss)\n\nencoder = Model([scaled_ts_input, bin_ts_input, add_factor, mult_factor], encoded)","metadata":{"_uuid":"b1dff87531d76d18faa32556fd3fbe9a65afdb81","scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In order to use for training of the autoencoder the whole data we feed it in chunks by using `fit_generator()`. Here the first loop is over the training set and the second loop is over the whole test set.  The batch size chosen to be 1000, but it might very well be that bigger or smaller batch size provides better results in terms of computational time and accuracy.","metadata":{"_uuid":"7981ad829431800eb5ae78eb6dd3d461b49eb476"}},{"cell_type":"code","source":"from keras.callbacks import EarlyStopping\n\ndef time_series_generator(train_path, test_path, chunk_size,\\\n                          steps_per_epoch, train_n=sys.maxsize, test_n=sys.maxsize):\n    gc.enable()\n    while 1:\n        with codecs.open(train_path, \"r\", ENCODING) as fp:\n            time_series_redaer =iter_time_serie_by_id(fp, chunk_size)\n            for (time_series, ids) in time_series_redaer:\n                (scaled_ts, add_factor, mult_factor) = rescale_time_serie(time_series)\n                yield ([scaled_ts, time_series[:,2:5,], add_factor, mult_factor], time_series[:,[0,1,4],])\n                gc.collect() \n        with codecs.open(test_path, \"r\", ENCODING) as fp:\n            time_series_redaer =iter_time_serie_by_id(fp, chunk_size)\n            for (time_series, ids) in time_series_redaer:\n                (scaled_ts, add_factor, mult_factor) = rescale_time_serie(time_series)\n                yield ([scaled_ts, time_series[:,2:5,], add_factor, mult_factor], time_series[:,[0,1,4],])\n                gc.collect()\nearlystop = EarlyStopping(monitor='loss', min_delta=0.0001, patience=3, \\\n                          verbose=1, mode='auto')\ncallbacks_list = [earlystop]\n\n# We comment out the whole number of items to be used for training, so that \n# the computation for the commit is faster\n#TRAIN_N = (N-VAL_COUNT-1)\n#TEST_N = 10000\n\nCHUNK_SIZE = 2 * 1024\nSTEPS_PER_EPOCH = 2 * 853\n\nautoencoder.fit_generator(\n    time_series_generator(TRAIN_PATH, TEST_PATH, CHUNK_SIZE, STEPS_PER_EPOCH, TRAIN_N, TEST_N),\n    steps_per_epoch = STEPS_PER_EPOCH,\n    epochs=4,\n    max_queue_size = 3,\n    shuffle=True,\n    #validation_data=([val_scaled, val_ts[:,2:5,], val_add_factor, val_mult_factor], val_ts[:,[0,1,4],]),\n    callbacks=callbacks_list)\n\n# serialize model to JSON\nencoder_json = encoder.to_json()\nautoencoder_json = autoencoder.to_json()\nwith open(\"encoder_bs-{}_loss-mse_el-{}.json\".format(CHUNK_SIZE, ENCODING_DIM), \"w\") as json_file:\n    json_file.write(encoder_json)\nwith open(\"autoencoder_bs-{}_loss-mse_el-{}.json\".format(CHUNK_SIZE, ENCODING_DIM), \"w\") as json_file:\n    json_file.write(autoencoder_json)\n# serialize weights to HDF5\nencoder.save_weights(\"encoder_bs-{}_loss-mse_el-{}.h5\".format(CHUNK_SIZE, ENCODING_DIM))\nautoencoder.save_weights(\"autoencoder_bs-{}_loss-mse_el-{}.h5\".format(CHUNK_SIZE, ENCODING_DIM))\nprint(\"Saved model to disk\")","metadata":{"_uuid":"30eee6a692106317ab103d8381f159df5ff16f95","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Below is demonstration how to use the `encoder` for extracting features.","metadata":{"_uuid":"e90a6cc3f96bfd12eb3436f3832f650827514dd8"}},{"cell_type":"code","source":"#embedded_features = encoder.predict([val_scaled, val_ts[:,2:5,], val_add_factor, val_mult_factor])","metadata":{"_uuid":"b83a2c5d098825cfd31de11641896e221ce6243a","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We visually inspect how the autoencoder is performing with a reconstruction of the time-series.  `train_ts` has time-series which have been included into the training, whereas `val_ts` have been just used for validation and we can see how the `encoder` reconstructs unseen time-series. As we can see the reconstruction is good in terms of pattern, but the scale is shifted. Analysis will be done to see how to reconstruct the exact scale.  If you dear kaggelers spot the reason for the wrong scale, shout =D","metadata":{"_uuid":"2138f8c2bd65d6259fd019acd0a72fad847548ad"}},{"cell_type":"code","source":"'''\nimport matplotlib.pyplot as plt\nval_obj_id = 0\ntrain_obj_id = 0\npassband = 1\nfeature =  1\n\n# The indicator of observed values in the time series is used as a mask\nval_mask = val_ts[val_obj_id,feature,:,passband]\nrec_val_ts = autoencoder.predict([val_scaled, val_ts[:,2:5,], val_add_factor, val_mult_factor])\nplt.plot(val_ts[val_obj_id,feature,:,passband])\nplt.show()\n# masking the unobserved part of the time-serie, which we don't care about\nplt.plot(rec_val_ts[val_obj_id,feature,:,passband] * val_mask)\nplt.show()\n\n# The indicator of observed values in the time series is used as a mask\ntrain_mask = train_ts[train_obj_id,feature,:,passband]\nrec_train_ts = autoencoder.predict([train_scaled, train_ts[:,2:5,], train_add_factor, train_mult_factor])\n\nplt.plot(train_ts[train_obj_id,feature,:,passband])\nplt.show()\n# masking the unobserved part of the time-serie, which we don't care about\nplt.plot(rec_train_ts[train_obj_id,feature,:,passband] * train_mask)\nplt.show()\n'''","metadata":{"_uuid":"8e8b4eddbf1a7268ffdcafc310dbebf0a4d08bf1","trusted":true},"execution_count":null,"outputs":[]}]}