{"cells":[{"metadata":{},"cell_type":"markdown","source":"# General information\nIn this notebook, we explore the statistical distribution of the reference signal (`time_to_failure`) and the implications of chopping the signals in blocks of 150.000 samples."},{"metadata":{},"cell_type":"markdown","source":"# Preliminaries\nLet's import everything we need:"},{"metadata":{"trusted":false},"cell_type":"code","source":"import gc\nimport os\nimport time\nimport random\nimport datetime\nimport warnings\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom sklearn.model_selection import StratifiedKFold, KFold, RepeatedKFold, GridSearchCV, cross_val_score","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's set some configurations"},{"metadata":{"trusted":false},"cell_type":"code","source":"%matplotlib inline\npd.options.display.precision = 15\nrandom.seed(6) #totally random seed (got from a dice)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training data\nWe load the raw input data as is."},{"metadata":{"trusted":false},"cell_type":"code","source":"%%time\ntrain = pd.read_csv('../input/train.csv', dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\nfs = 4000000 #sampling frequency of the sensor signal","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"print(f'Train: rows:{train.shape[0]} cols:{train.shape[1]}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's visualize the data"},{"metadata":{"trusted":false},"cell_type":"code","source":"train_acoustic_data_small = train['acoustic_data'].values[::50]\ntrain_time_to_failure_small = train['time_to_failure'].values[::50]\n\nfig, ax1 = plt.subplots(figsize=(16, 8))\nplt.title('Acoustic_data and time_to_failure (sampled)')\nplt.plot(train_acoustic_data_small, color='b')\nax1.set_ylabel('acoustic_data', color='b')\nplt.legend(['acoustic_data'])\nax2 = ax1.twinx()\nplt.plot(train_time_to_failure_small, color='g')\nax2.set_ylabel('time_to_failure', color='g')\nplt.legend(['time_to_failure'], loc=(0.875, 0.9))\nplt.grid(False)\n\ndel train_acoustic_data_small\ndel train_time_to_failure_small","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It seems that the reference data (green line) won't be uniformly distributed. We can visualize it "},{"metadata":{"scrolled":false,"trusted":false},"cell_type":"code","source":"plt.figure()\nplt.hist(train['time_to_failure'].values[::50], bins='auto', density=True)  # arguments are passed to np.histogram\nplt.title('Histogram of time_to_failure')\nplt.xlabel('time_to_failure')\nplt.ylabel('count')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The time to training is not uniforally distribured on the training set. Resampling could be interesting."},{"metadata":{},"cell_type":"markdown","source":"# Sampling the training set\nSince the testign set is build of sequences of 150.000 samples, let's chop the training set on comparable chunks. \n\nSeveral strategies of how to select the chunks are possible:\n- Deterministic\n- Random\n\nMoreover, the *jumps* on the `time_to_failure` could be an issue (to be seen the impact). In this situation, we can\n- Include the *jumps* on the training set\n- Exlcude the *jumps* on the training set\n\nLet's see what happens to the histogram in some of these situations:"},{"metadata":{"trusted":false},"cell_type":"code","source":"# Create a training file with simple derived features\nsegment_size = 150000","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"def generate_segment_start_ids(sampling_method):\n    \"\"\" Generates the indeces where the segments for the training data start \"\"\"\n    if sampling_method == 'uniform':\n        # With this approach we obtain 4194 segments\n        num_segments_training = int(np.floor(train.shape[0] / segment_size))\n        segment_start_ids = [i * segment_size for i in range(num_segments_training)]\n    elif sampling_method == 'uniform_no_jump':\n        # With this approach we obtain 4178 segments (99.5% of 'uniform')\n        already_sampled = np.full(train.shape[0], False)\n        num_segments_training = int(np.floor(train.shape[0] / segment_size))\n        time_to_failure_jumps = np.diff(train['time_to_failure'].values)\n        num_good_segments_found = 0\n        segment_start_ids = []\n        for i in range(num_segments_training):\n            idx = i * segment_size\n            # Detect if there is a discontinuity on the time_to_failure signal within the segment\n            max_jump = np.max(time_to_failure_jumps[idx:idx + segment_size])\n            if max_jump < 5:\n                segment_start_ids.append(idx)\n                num_good_segments_found += 1\n            else:\n                print(f'Rejected candidate segment since max_jump={max_jump}')\n        segment_start_ids.sort()\n    elif sampling_method == 'random_no_jump':\n        # With this approach we obtain 4194 segments\n        num_segments_training = int(np.floor(train.shape[0] / segment_size)) #arbitrary choice\n        time_to_failure_jumps = np.diff(train['time_to_failure'].values)\n        num_good_segments_found = 0\n        segment_start_ids = []\n        while num_segments_training != num_good_segments_found:\n            # Generate a random sampling position\n            idx = random.randint(0, train.shape[0] - segment_size - 1)\n            # Detect if there is a discontinuity on the time_to_failure signal within the segment\n            max_jump = np.max(time_to_failure_jumps[idx:idx + segment_size])\n            if max_jump < 5:\n                segment_start_ids.append(idx)\n                num_good_segments_found += 1\n            else:\n                print(f'Rejected candidate segment since max_jump={max_jump}')\n        segment_start_ids.sort()\n    else:\n        raise NameError('Method does not exist')\n    return segment_start_ids","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Currently, we have three strategies implemented.\n- Uniform sampling: Here the training set is built from consecutive chunks of data. No special care is taken into what's the content of each segment.\n- Uniform sampling with rejection: Here the training set is built from consecutive chunks of data. There is a control to avoid having in a segment the jump of the `time_to_failure` signal from zero to a high value. *This is very dangerous when doing the splitting between training and validation*.\n- Random sampling with rejection: Here the training set is built from sampling randomly the data. There is a control to avoid having in a segment the jump of the `time_to_failure` signal from zero to a high value. The segments will most likely overlap."},{"metadata":{},"cell_type":"markdown","source":"Let's visualize what happens to the distribution of the `time_to_failure` on teh different samplign strategies:"},{"metadata":{"scrolled":false,"trusted":false},"cell_type":"code","source":"print(f'Generating uniformly sampled training set')\nsegment_start_ids_uniform = generate_segment_start_ids('uniform')\n\nprint(f'Generating uniformly sampled training set excluding discontinuities in time_to_failure.')\nsegment_start_ids_uniform_no_jump = generate_segment_start_ids('uniform_no_jump')\n\nprint(f'Generating randomly sampled training set excluding discontinuities in time_to_failure.')\nprint(f'This method may yield overlaping segments')\nsegment_start_ids_random_no_jump = generate_segment_start_ids('random_no_jump')\n\n\ny_tr_samples_uniform = train['time_to_failure'].values[np.array(\n    segment_start_ids_uniform) + segment_size - 1]\ny_tr_samples_uniform_no_jump = train['time_to_failure'].values[\n    np.array(segment_start_ids_uniform_no_jump) + segment_size - 1]\ny_tr_samples_random_no_jump = train['time_to_failure'].values[\n    np.array(segment_start_ids_random_no_jump) + segment_size - 1]\n\nplt.subplots(figsize=(16, 5))\nplt.subplot(1, 3, 1)\nplt.hist(y_tr_samples_uniform, bins='auto', alpha=0.5, density=True)\nplt.hist(train['time_to_failure'].values[::50], bins='auto', alpha=0.5, density=True)\nplt.title('With discontinuities (contiguous)')\nplt.legend(['Sampled', 'All data'])\n\nplt.subplot(1, 3, 2)\nplt.hist(y_tr_samples_uniform_no_jump, bins='auto', alpha=0.5, density=True)\nplt.hist(train['time_to_failure'].values[::50], bins='auto', alpha=0.5, density=True)\nplt.title('Discarding discontinuities (contiguous)')\nplt.legend(['Sampled', 'All data'])\n\nplt.subplot(1, 3, 3)\nplt.hist(y_tr_samples_random_no_jump, bins='auto', alpha=0.5, density=True)\nplt.title('Discarding discontinuities (rand)')\nplt.hist(train['time_to_failure'].values[::50], bins='auto', alpha=0.5, density=True)\nplt.legend(['Sampled', 'All data'])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we conclude that the histogram of the `time_to_failure` does not change significantly when excluding *all* jumps. However, the histogram starts to be significantly distorded when randon sampling is applied."},{"metadata":{},"cell_type":"markdown","source":"[@konradb](https://www.kaggle.com/konradb), on a [post](https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/85108#496570) on March 22ds, he pointed out that \"By adjusting the median of my prediction to match the target, I got a boost both in cv (2.0702 -> 2.0538) and lb (1.512 -> 1.489).\"\nLet's explore what happens to the distribution when making a k-fold."},{"metadata":{"trusted":false},"cell_type":"code","source":"n_fold = 5\nfolds = KFold(n_splits=n_fold, shuffle=True, random_state=11)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"trusted":false},"cell_type":"code","source":"n_fold = folds.get_n_splits()\n\nfor fold_n, (train_index, valid_index) in enumerate(folds.split(y_tr_samples_uniform_no_jump)):\n    print('Fold', fold_n, 'started at', time.ctime())\n    y_train, y_valid = y_tr_samples_uniform_no_jump[train_index], y_tr_samples_uniform_no_jump[valid_index]\n\n    plt.figure()\n    plt.hist(y_train, bins='auto', alpha=0.5, density=True)\n    plt.hist(y_valid, bins='auto', alpha=0.5, density=True)\n    plt.title(f\"Histogram with fold {fold_n}\")\n    plt.legend([f'Train (median={np.median(y_train):.4f})', f'Validation (median={np.median(y_valid):.4f})'])\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we see here, the distributions of the training and validation datasets are not identical. It is to be seen if a more clever split can be done. "},{"metadata":{},"cell_type":"markdown","source":"Next questions to expore:\n- For training, should we augment the dataset to have more equalized histogram?"},{"metadata":{"trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"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.8"}},"nbformat":4,"nbformat_minor":1}