{"cells":[{"metadata":{},"cell_type":"markdown","source":"# EDA for Volcanic Eruption Competition\n\n### Goals\n\n* For competition https://www.kaggle.com/c/predict-volcanic-eruptions-ingv-oe/data\n\n### Comments"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom pathlib import Path\nimport sqlite3\nfrom tqdm.notebook import tqdm\ntqdm.pandas()\nimport missingno as msno\nfrom multiprocessing import Pool\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nsns.set_context(\"talk\", font_scale=1.4)\nsns.set_style('whitegrid')\n\n\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nimport os\n\n# You can write up to 5GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!apt-get install tree","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step I: Business Goal\n\n* To allow for an early warning of vulcanic erruption the goal is to predict when a volcano's next eruption will occur. \n* You'll analyze a large geophysical dataset collected by sensors deployed on active volcanoes\n* Identify signatures in seismic waveforms that characterize the development of an eruption. \n\n* Metrics: MAE for the time until erruption\n* Final submission deadline: December 30, 2020\n\n* Hypotheses\n    * Sensor signatures indicate the development towards an erruption\n    * The sensor signatures for eruptions are similar across different vulcanos\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Step II: Data Extraction\n\n* readings are from several seismic sensors around a volcano and challenges you to estimate how long it will be until the next eruption. \n* The data represent a classic signal processing setup that has resisted traditional methods.\n* Each file contains 10 minutes of logs from 10 different sensors arrayed around a volcano.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"! tree /kaggle/input/predict-volcanic-eruptions-ingv-oe -L 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"erruptions = pd.read_csv(\"/kaggle/input/predict-volcanic-eruptions-ingv-oe/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"erruptions.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_path_train = []\nsensors_path_test = []\n\n\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        path = os.path.join(dirname, filename)\n        if 'train/' in path:\n            sensors_path_train.append(Path(path))\n        if 'test/' in path:\n            sensors_path_test.append(Path(path))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train = pd.Series(sensors_path_train).to_frame('path')\n\nsensors_test = pd.Series(sensors_path_test).to_frame('path')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train['segment_id'] = sensors_train['path'].apply(lambda path: int(path.stem))\nsensors_test['segment_id'] = sensors_test['path'].apply(lambda path: int(path.stem))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt = pd.merge(sensors_train, erruptions, on ='segment_id', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Read one 10-min file/sensor time series for exploration "},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample = pd.read_csv(sensors_train['path'].iloc[0])\nsensor_sample.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"f'Total number of sensor rows: {sensors_train.shape[0] * sensor_sample.shape[0]} x 10 sensor data columns'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Dataset too large to load into the Kaggle Machine memory of 16GB (cpu) at once."},{"metadata":{"trusted":true},"cell_type":"code","source":"storage_sql=False\n\nif storage_sql:\n\n    create_db=False\n    if create_db:\n        conn = sqlite3.connect('train2.db')\n        print(\"Opened database successfully\")\n\n        conn.execute('''CREATE TABLE sensors\n                 (SERIE         INT    NOT NULL,\n                 S1           REAL    NOT NULL,\n                 S2           REAL    NOT NULL,\n                 S3           REAL    NOT NULL,\n                 S4           REAL    NOT NULL,\n                 S5           REAL    NOT NULL,\n                 S6           REAL    NOT NULL,\n                 S7          REAL    NOT NULL,\n                 S8           REAL    NOT NULL,\n                 S9           REAL    NOT NULL,\n                 S10           REAL    NOT NULL,\n                 ERRUPT          REAL    NOT NULL);''')\n        print(\"Table created successfully\")\n        conn.close()\n    \n    def write_to_sql(bunch):\n        \n        assert type(bunch[0]) == tuple\n        # [(serie, S1, S2, S3, S4, S5... , ERRUPT)]\n        #examples = [(2, \"def\"), (3, \"ghi\"), (4, \"jkl\")]\n        cur.executemany(\"INSERT INTO sensors VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)\", bunch)\n    \n    def ingest_data(row):\n        sensor_series = pd.read_csv(row['path'])\n\n        sensor_series['time_to_eruption'] = row['time_to_eruption']\n        sensor_series['segment_id'] = str(row['segment_id'])\n        #sensor_series['adsfb'] = 'abc'\n        #print(sensor_series.dtypes)\n\n        sensor_series = sensor_series[['segment_id'] + sensor_series.columns.to_list()[:-1]]\n\n        to_ingest = [tuple(x) for i, x in sensor_series.iterrows()]\n\n\n\n\n        #to_ingest = []\n        #for i, row_sens in sensor_test.iterrows():\n        #    to_ingest.append(tuple([row['segment_id']] + row_sens.to_list() + [row['time_to_eruption']]))\n\n        #print(to_ingest)\n        write_to_sql(to_ingest)   \n    \n    conn = sqlite3.connect('train2.db')\n    print(\"Opened database successfully\")\n    \n    _ = sensors_train_erupt.head(1).progress_apply(ingest_data, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# II. Meet and Greet the Data\n\nIn the following I analyze the high-level sensor data as well as 1 sample cycle file.\n\n\n* Missing values identified via pandas: \n    * the sensor sample cycle file has sensor_9 with about 1/6 of missing data.\n    * All train cycle series have a time-to-eruption, no missing values\n\n* According to sample file sensor values are varying, can be positive and negative.\n\n* Data types in proper format, tested sensor cycle sample\n    * segment_id: discrete numeric identifier\n    * time to erruption: discrete numeric. Unit is not defined.\n    * sensor data: \n        * documentation suggests \"you may find that you still need to load the data as float32 due to the presence of some nulls\"\n \n* Assumptions:\n    *  The readings have been normalized within each segment, in part to ensure that the readings fall within the range of int16 values.\n    \n"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.sample(10, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.sample(10, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.info()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Missing Values? Missigness not randomly distributed but occurs in certain time intervals."},{"metadata":{"trusted":true},"cell_type":"code","source":"msno.matrix(sensor_sample)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Working with all sensor cycle files\n\n* Check for number of steps per cycle: Is it always 60001 as in sample above? - yes always.\n* Always same data types for cycle files? - yes.\n* how much missing values are there in other cycle files? Missingness is common, 10% ahve missing values, 5% of sensors have no sensor values at all. Best to have algorithm to deal with missing values."},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_aggregates(sensor_cycle):\n        \n    # series sensor 1 to 10\n    missing_values = sensor_cycle.isna().sum().values\n    \n    cycle_steps = sensor_cycle.shape[0]\n    \n    dtypes = [str(val) for val in sensor_cycle.dtypes.values]\n    \n    return [cycle_steps, dtypes, missing_values]\n\ndef process_cycle(path):\n    \n    sensor_cycle = pd.read_csv(path)\n    \n    return extract_aggregates(sensor_cycle)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"paths_train = [str(path) for path in sensors_train_erupt['path'].values]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with Pool(4) as p:\n  aggs = list(tqdm(p.imap(process_cycle, paths_train), total=len(paths_train)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(aggs)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats = pd.DataFrame(aggs, columns=['num_rows', 'dtypes', 'missingness'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats = sensors_train_erupt.join(aggs_stats)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats['num_rows'].unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats['dtypes'].apply(lambda dtypes_cycle:  dtypes_cycle== aggs_stats['dtypes'].iloc[0]).all()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Missing Values"},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_cycles = pd.DataFrame(np.array(aggs_stats['missingness'].to_list()), columns=sensor_sample.columns)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_cycles.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If missing values for sensor, are there all sensor values for the whole cycle missing?"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_missing_values=  missing_cycles.applymap(lambda x: x>0).sum().sum()\nprint(f'{round(100*sensors_missing_values/(missing_cycles.shape[0]*missing_cycles.shape[1]), 2)}% of sensors have missing data')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_missing_values=  missing_cycles.applymap(lambda x: x>0).sum().sum()\nsensors_all_missing_values=  missing_cycles.applymap(lambda x: x==60001).sum().sum()\nprint(f'{round(100*sensors_all_missing_values/sensors_missing_values, 2)}% of sensors with missing data have not sensor values at all.')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Missingess Analysis:\n* Different cycles have different sensor values missing.\n* We expect that there will be data in test set that have no sensor values at all. More difficult to impute, if it becomes necessary."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(20,8))\nsns.heatmap(missing_cycles, ax=ax)\nax.set_ylabel('sensor cycles')\nax.set_title('Number of missing values accross sensor cycles')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax = missing_cycles.applymap(lambda val: val>0).sum().div(missing_cycles.shape[0]).plot(kind='bar')\nax.set_title('fraction cycles with missing values')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Some sensors seem to be more prone towards having missing values. Particular `sensor_9` has in 30% of all cycles missing values."},{"metadata":{"trusted":true},"cell_type":"code","source":"missing_cycles.applymap(lambda val: val>0).sum(axis=1).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats[missing_cycles.applymap(lambda val: val>0).sum(axis=1)==10]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A some cycles below have more than 4+ sensors with signficant missing values. Missingness is not a big problem for other cycles."},{"metadata":{"trusted":true},"cell_type":"code","source":"aggs_stats[missing_cycles.applymap(lambda val: val>0).sum(axis=1)>3][['segment_id', 'missingness']]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del aggs_stats","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del aggs","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step IV: Univariate Analysis\n\n* Sensor cycle sample insights\n    * no obvious correlations but appear non-linear affects\n    * "},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt[['time_to_eruption']].describe().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.axes_style(\"ticks\"):\n    fig, ax = plt.subplots(figsize=(16,4))\n    sns.distplot(np.log10(sensors_train_erupt['time_to_eruption']), ax=ax, kde=True, hist=True)\n    #ax.set_xscale('log')\n    ax.set_xlabel('log time to eruption')\n    plt.yticks(np.arange(0,2,0.25))\n    plt.grid()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(16,4))\nsns.boxplot(np.log10(sensors_train_erupt['time_to_eruption']), ax=ax)\nax.set_xlabel('log time to eruption')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Outliers using IQR for short times until erruption."},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt[['path', 'segment_id']].astype('O').describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Investigating Sensor Cycle Sample\n\n* All data appears normalized to median equal to zero (reference normalization)\n* "},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.describe().T","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Is there any correlation between sensor data?"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(sensor_sample)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> No linear correlation between sensor values but there appears to be some non-linear correlations (shapes peculiar)."},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.plotting_context(\"talk\", font_scale=1):\n    fig, ax = plt.subplots(figsize=(16,8))\n    sns.heatmap(sensor_sample.corr(), annot=True, fmt=\".2\", ax=ax)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"to_plot = sensor_sample.unstack() \nto_plot = to_plot.droplevel(1).reset_index()\n\ng = sns.FacetGrid(to_plot, col='index', col_wrap=3, sharex=False, sharey=False, aspect=1.2, height=6)\ng.map(sns.distplot, 0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Time Series Analysis\n\nstep size: 60x10x100 = 10min > 10ms intervals between sensor readings"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.plotting_context(\"talk\", font_scale=1):\n    fig, axes = plt.subplots(len(sensor_sample.columns), figsize=(30, 20), sharex=True)\n    axes = axes.flatten()\n    for i, col in enumerate(sensor_sample.columns):\n        sns.lineplot(data=sensor_sample, x=sensor_sample.index, y=col, ax=axes[i])\n        axes[i].set_ylabel(col); axes[i].set_xlabel('')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Normalization effect visible. \n    * There appear to be signals (e.g. sensor_5 and sensor 10 have larger variances for certain time periods) and there does not appear to be only simple white noise. We measure something! \n    * Also, magnitude of change differs and larger variance differences, e.g. sensor1 vs sensor2. \n    * Also, no seasonal effects expected so do not do any decomposition\n    * Also, no investigation into autoregression parts and other things needed due to the nature of the problem.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.isna().sum()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rolmean = sensor_sample.rolling(window=100).mean()# 1s averages","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(16,8))\nplt.plot(sensor_sample['sensor_1'], color='blue', label='Actual Series')\nplt.plot(rolmean['sensor_1'], color='red', label='Actual Series')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ADF test\ndef adf_test(series):\n    from statsmodels.tsa.stattools import adfuller\n    test = adfuller(series)\n    output = pd.Series(test[0:4], index=['Test Statistic','p-value','Lags Used','No. of Observations'])\n    for key,value in test[4].items():\n        output['Critical Value (%s)'%key] = value\n    print(output)\nadf_test(sensor_sample['sensor_1'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Test reveals stationarity while our mean estimation  over the 1-s window did not show stationarity."},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.date_range(\"00:00\", \"01:00\", freq=\"10ms\").time.shape","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Extraction Sensor Features and comparing with predictor\n\n* Using quantile binning to extract features from time series/summarise the time series for each cycle. Maybe choose other quantiles in the future.\n* Investigate the pearson correlation between sensor values for each cycle."},{"metadata":{},"cell_type":"markdown","source":"#### 1. Extraction"},{"metadata":{"trusted":true},"cell_type":"code","source":"def extract_stats(sensor_cycle, cycle_ID):\n        \n    # extract sensors stats for cycles \n    description = sensor_cycle.describe().T\n    description['segment_id'] = cycle_ID\n    \n    # identify cycles with linear correlation\n    corr = sensor_cycle.corr().fillna(0).values\n    np.fill_diagonal(corr, 0)\n    abs_threshold = 0.4\n    corr_signal = (np.abs(corr)>abs_threshold).sum().sum() /2. # factor 2 due to symmetric matrix\n    \n    \n    return [description, corr_signal]\n\ndef process_cycle_stats(path):\n    \n    sensor_cycle = pd.read_csv(path)\n    \n    return extract_stats(sensor_cycle, path.stem)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with Pool(4) as p:\n  aggs = list(tqdm(p.imap(process_cycle_stats, sensors_train_erupt['path']), total=len(paths_train)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"summary_stats_cycles = pd.concat([val[0] for val in aggs])\nsummary_stats_cycles.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"summary_stats_cycles['segment_id'] = summary_stats_cycles['segment_id'].astype(int)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### 2. Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"summary_stats_cycles.head().drop(columns=['segment_id', 'count']).stack().to_frame().T","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"stats_by_cycle = summary_stats_cycles.groupby('segment_id').apply(lambda cycle: cycle.drop(columns=['segment_id', 'count']).stack().to_frame().T)\n\ncols = [col[0]+'_'+col[1] for col in stats_by_cycle.columns]\nstats_by_cycle.columns = cols\nstats_by_cycle.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Distribution of maximum and minimum values per cycle"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_max = [col for col in stats_by_cycle.columns if 'max' in col]\nto_plot = stats_by_cycle[sensor_max].stack().reset_index().drop(columns=['level_1', 'segment_id'])\n\nwith sns.plotting_context(\"talk\", font_scale=0.7):\n    g = sns.FacetGrid(to_plot, col=\"level_2\", col_wrap=3, sharex=False, sharey=False,  aspect=1.5, height=4)\n    g.map(sns.distplot, 0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> normalish shaped with long tail towards high values. Outliers indicated. It could be beneficial to do log transform or general power-transform box-cox"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_max = [col for col in stats_by_cycle.columns if 'max' in col]\nto_plot = stats_by_cycle[sensor_max].stack().reset_index().drop(columns=['level_1', 'segment_id'])\n\nwith sns.plotting_context(\"talk\", font_scale=0.7):\n    g = sns.FacetGrid(to_plot, col=\"level_2\", col_wrap=3, sharex=False, sharey=False,  aspect=1.5, height=4)\n    g.map(sns.boxplot, 0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_min = [col for col in stats_by_cycle.columns if 'min' in col]\nto_plot = stats_by_cycle[sensor_min].stack().reset_index().drop(columns=['level_1', 'segment_id'])\n\nwith sns.plotting_context(\"talk\", font_scale=0.7):\n    g = sns.FacetGrid(to_plot, col=\"level_2\", col_wrap=3, sharex=False, sharey=False,  aspect=1.5, height=4)\n    g.map(sns.distplot, 0)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Again, clearly there are strong outliers for some sensor series with some very large and very small values. This could be badly calibrated sensors or missfunctional sensors. Should I remove those? These outliers could also be related to the time to erruption as closer erruption time leads to more and higher signals?"},{"metadata":{},"cell_type":"markdown","source":"#### Outlier sensors"},{"metadata":{"trusted":true},"cell_type":"code","source":"def outlier_iqr(series):\n\n    q75=series.quantile(q=0.75)\n    q25=series.quantile(q=0.25)\n    IQR = q75-q25\n    low_IQR = q25 -1.5*IQR\n    high_IQR = q75+1.5*IQR\n    #print(low_IQR, high_IQR)\n    to_keep = stats_by_cycle[(series>low_IQR) & (series<high_IQR)]\n    print(f'drop {series.shape[0]-to_keep.shape[0]} out of {series.shape[0]}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"outlier_iqr(stats_by_cycle['sensor_1_max'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"quantile of 75% might be better suited to identify vastly different cycles/sensor setups:"},{"metadata":{"trusted":true},"cell_type":"code","source":"outlier_iqr(stats_by_cycle['sensor_1_75%'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_75 = [col for col in stats_by_cycle.columns if '75%' in col]\nto_plot = stats_by_cycle[sensor_75].stack().reset_index().drop(columns=['level_1', 'segment_id'])\n\nwith sns.plotting_context(\"talk\", font_scale=0.7):\n    g = sns.FacetGrid(to_plot, col=\"level_2\", col_wrap=3, sharex=False, sharey=False,  aspect=1.5, height=4)\n    g.map(sns.boxplot, 0)\n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step V: Multivariate Analysis\n\n* Test how extracted features depend on erruption time.\n\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt_stats = pd.merge(sensors_train_erupt, stats_by_cycle.reset_index(), \n                                     on='segment_id', how='left').drop(columns=['path', 'level_1'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### 2. Analyis of extracted features"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt_stats.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cols_plot = [col for col in sensors_train_erupt_stats.columns if '25%' in col]\ncols_plot","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt_stats.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.plotting_context(\"talk\", font_scale=0.9):\n    fig, ax = plt.subplots(figsize=(20,4))\n    sensors_train_erupt_stats.corr()['time_to_eruption'].plot(kind='bar', ax=ax)\n    #ax.xaxis.set_visible(False)\n    plt.yticks([-0.5, 0, 0.5])\n    ax.set_ylabel('person correlation')\n    ax.set_title('correlation of sensor features with time to erruption')\n    ax.set_ylim(-0.5, 0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> There are not strong correlations but weak correlation of some features. This weak correlation varies between sensors, e.g. std, 25, 75, max. Median is always zero due to normalization and mean is also very small. "},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt_stats['time_to_eruption_log'] = np.log10(sensors_train_erupt_stats['time_to_eruption'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Investigating the maximum values per series. Could investigate more sensor stats."},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.plotting_context(\"talk\", font_scale=0.5):\n\n    col_max_features = [col for col in sensors_train_erupt_stats.columns if 'max' in col]\n    fig, axes = plt.subplots(len(col_max_features), figsize=(20, 14))\n    axes = axes.flatten()\n    for i, col in enumerate(col_max_features):\n        sns.scatterplot(data=sensors_train_erupt_stats, x='time_to_eruption_log', y=col, ax=axes[i])\n    ax.set_title('Maximum sensor values for eruption times')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.plotting_context(\"talk\", font_scale=0.5):\n\n    col_max_features = [col for col in sensors_train_erupt_stats.columns if 'mean' in col]\n    fig, axes = plt.subplots(len(col_max_features), figsize=(20, 14))\n    axes = axes.flatten()\n    for i, col in enumerate(col_max_features):\n        sns.scatterplot(data=sensors_train_erupt_stats, x='time_to_eruption_log', y=col, ax=axes[i])\n    ax.set_title('Average sensor values for eruption times')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### One could investigate how the correlation of sensor values relates to the time until eruption.\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt['corr_sensors'] = [val[1] for val in aggs]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Most cycle sensors have no linear correlation. Less than 2% have and are hence not relevant here."},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt['corr_sensors'].value_counts(normalize=True).head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Extraction of Time series for certain times until erruption."},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"with sns.axes_style(\"ticks\"):\n    fig, ax = plt.subplots(figsize=(16,4))\n    sns.distplot(np.log10(sensors_train_erupt['time_to_eruption']), ax=ax, kde=True, hist=True)\n    #ax.set_xscale('log')\n    ax.set_xlabel('log time to eruption')\n    plt.yticks(np.arange(0,2,0.25))\n    plt.grid()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(16,4))\nsns.boxplot(np.log10(sensors_train_erupt['time_to_eruption']), ax=ax)\nax.set_xlabel('log time to eruption')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt['time_to_eruption_log'] = np.log10(sensors_train_erupt['time_to_eruption'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Choose quantiles to represent log distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.qcut(sensors_train_erupt['time_to_eruption_log'], q=[1e-8, 1e-7, 1e-4, 1e-3, 1e-2, 0.1, 1]).cat.categories","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.qcut(sensors_train_erupt['time_to_eruption_log'], q=[1e-4, 1e-3, 1e-2, 0.1, 1]).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt['q_time_to_eruption'] = pd.qcut(sensors_train_erupt['time_to_eruption_log'], q=[1e-4, 1e-3, 1e-2, 0.1, 1], \n                                                    labels=[1e-3, 1e-2, 0.1, 1]).astype(float)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt['q_time_to_eruption'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.sort_values('time_to_eruption_log').head(1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.loc[sensors_train_erupt['segment_id'] == 601524801, 'q_time_to_eruption'] = 1e-7","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sample cycle for each time_to_eruption bin:"},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples = sensors_train_erupt.groupby('q_time_to_eruption').apply(lambda x: x.sample(1, random_state=42)).drop(\n    columns='q_time_to_eruption').reset_index().sort_values('q_time_to_eruption').drop(columns='level_1')\ncycle_samples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Get all sensor data for the samples"},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensors = []\nfor i, row in cycle_samples.iterrows():\n    sens_tmp = pd.read_csv(row['path'])\n    sens_tmp['segment_id'] = row['segment_id']\n    sens_tmp['step'] = sens_tmp.index\n    cycle_samples_sensors.append(sens_tmp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensors = pd.merge(pd.concat(cycle_samples_sensors), cycle_samples, on='segment_id', how= 'left').drop(columns='path')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensors.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensors.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_cols = [col for col in cycle_samples_sensors.columns if 'sensor_' in col]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_cols_plot = cycle_samples_sensors.melt(id_vars=['q_time_to_eruption', 'time_to_eruption_log', 'step'], \n        var_name=\"sensors\", value_vars=sensor_cols,\n        value_name=\"Value\")\nsensor_cols_plot.head(2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_cols_plot['time_to_eruption_log'] = sensor_cols_plot['time_to_eruption_log'].apply(lambda x: round(x, 2))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensors.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\n\nj=-1\nwith sns.plotting_context(\"talk\", font_scale=0.8):\n    \n    fig, axes = plt.subplots(10, 1, figsize=(16,40))\n    axes = axes.flatten()\n    for sensor_type in sensor_cols:\n        j+=1\n    #sensor_type = 'sensor_1'\n        plot_subset = sensor_cols_plot[(sensor_cols_plot['sensors'] == sensor_type)]        \n        colors = ['y', 'r', 'b', 'g', 'k']\n        i=-1\n        for label, df in plot_subset.groupby('time_to_eruption_log'):\n            i+=1\n            df.reset_index(drop=True).Value.plot(ax=axes[j], label=label, color=colors[i])\n            ax=axes[j].set_title(sensor_type)\n            #df.Value.plot(kind=\"kde\", ax=axes[1], label=label, color=colors[i])\n            ax=axes[j].legend()\n            #axes[1].legend()\n            #axes[1].set_xlim(-2500, 2500)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Signals at different times until erruption do NOT differe greatly. One can spot some time ranges with signal patterns attributed to shorter times to erruption, at 3.8 and 4.46. But they are not clearly distinct from patterns at other times."},{"metadata":{},"cell_type":"markdown","source":"### Looking at 10 signals per binned time-to-erruption"},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt['q_time_to_eruption'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def draw_sample(df):\n    if df.shape[0] <40:\n        return df\n    else:\n        return df.sample(10, random_state=42)\n    \ncycle_samples = sensors_train_erupt.groupby('q_time_to_eruption').apply(draw_sample).drop(\n    columns='q_time_to_eruption').reset_index().sort_values('q_time_to_eruption').drop(columns='level_1')\ncycle_samples","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_sensor_samples(sensor):\n    \n    cycle_samples_sensors = []\n\n    for i, row in cycle_samples.iterrows():\n        sens_tmp = pd.read_csv(row['path'])\n        sens_tmp['segment_id'] = row['segment_id']\n        sens_tmp['step'] = sens_tmp.index\n        cycle_samples_sensors.append(sens_tmp[[sensor, 'step', 'segment_id']])\n\n    cycle_samples_sensors = pd.merge(pd.concat(cycle_samples_sensors), cycle_samples, on='segment_id', how= 'left').drop(columns='path')\n    return cycle_samples_sensors","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensor = get_sensor_samples('sensor_1')\ncycle_samples_sensor.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cycle_samples_sensor.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looking just at sensor 1."},{"metadata":{"trusted":true},"cell_type":"code","source":"#colors = ['y', 'r', 'b', 'g', 'k']\n\nwith sns.plotting_context(\"talk\", font_scale=0.8):\n\n\n    fig, axes = plt.subplots(5, 1, figsize=(16,10), sharex=True, sharey=True)\n    fig.suptitle(\"sensor_1\")\n    \n    i=-1\n    for label, df in cycle_samples_sensor.groupby('q_time_to_eruption'):\n        i+=1\n        df.reset_index(drop=True).groupby('segment_id').plot('step','sensor_1',ax=axes[i])\n        \n        #sns.lineplot(data=df.reset_index(drop=True), x='step', y='sensor_1', hue='segment_id', ax=axes[i])\n        \n        ax=axes[i].set_title(f\"log time erruption: {round(df['time_to_eruption_log'].mean(), 2)}\")\n        #df.Value.plot(kind=\"kde\", ax=axes[1], label=label, color=colors[i])\n        ax=axes[i].legend().set_visible(False)\n        #axes[i].legend()\n        #axes[1].set_xlim(-2500, 2500)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> there are some cycles shown here where the sensor has extremely large values. This is related to my outlier analysis before. Again I ask: Are this simply differently calibrated sensors or malfunction sensors?\n\n> One can also see that the red seensor for erruptoin 6.35 seem to hit a threshold."},{"metadata":{"trusted":true},"cell_type":"code","source":"filter_extreme = cycle_samples_sensor.groupby('segment_id')['sensor_1'].describe()[\n    (cycle_samples_sensor.groupby('segment_id')['sensor_1'].describe()['75%'] > 1000)].index.to_list()\nfilter_extreme","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#colors = ['y', 'r', 'b', 'g', 'k']\n\nwith sns.plotting_context(\"talk\", font_scale=0.8):\n\n\n    fig, axes = plt.subplots(5, 1, figsize=(16,10), sharex=True, sharey=True)\n    fig.suptitle(\"sensor_1\")\n    \n    i=-1\n    for label, df in cycle_samples_sensor[~cycle_samples_sensor['segment_id'].isin(filter_extreme)].groupby('q_time_to_eruption'):\n        i+=1\n        df.reset_index(drop=True).groupby('segment_id').plot('step','sensor_1',ax=axes[i])\n        \n        #sns.lineplot(data=df.reset_index(drop=True), x='step', y='sensor_1', hue='segment_id', ax=axes[i])\n        \n        ax=axes[i].set_title(f\"log time erruption: {round(df['time_to_eruption_log'].mean(), 2)}\")\n        #df.Value.plot(kind=\"kde\", ax=axes[1], label=label, color=colors[i])\n        ax=axes[i].legend().set_visible(False)\n        #axes[i].legend()\n        axes[i].set_ylim(-2500, 2500)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Step VI. Downsample Series\n\n* Downsample series to 1-s intervals, which is a factor of 100 less in data.\n* write all downsampled series into one csv file.\n* Which kind of aggregation used for series?"},{"metadata":{"trusted":true},"cell_type":"code","source":"rolmean = sensor_sample.rolling(window=100).mean()# 1s averages\nplt.plot(sensor_sample['sensor_1'], color='blue', label='Actual Series')\nplt.plot(rolmean['sensor_1'], color='red', label='Actual Series')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.shape, rolmean.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.plot(rolmean['sensor_1'], color='red', label='Actual Series')\nplt.plot(rolmean[0:-1:100]['sensor_1'], color='green', label='Actual Series')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> A lot of high signal gets removed when creating a rolling average over 1s.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"rolmax = sensor_sample.rolling(window=100).max()# 1s averages\nplt.plot(sensor_sample['sensor_1'], color='blue', label='Actual Series')\nplt.plot(rolmax['sensor_1'], color='red', label='Actual Series')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rolsample = sensor_sample[0:-1:100]\nplt.plot(sensor_sample['sensor_1'], color='blue', label='Actual Series')\nplt.plot(rolsample['sensor_1'], color='red', label='Actual Series')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Above would correspond to a measure every 1s. Also notice the difference to the mean avg window.which does not preserve the max signal values.."},{"metadata":{"trusted":true},"cell_type":"code","source":"import csv","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_sample.cols","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"writer = csv.writer(open(\"low_freq_train.cv\", \"w\"))\nwriter.writerow([])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_train_erupt.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensor_cols = sensor_sample.columns.to_list()\nsensor_cols","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fname_train = \"train_low_freq.csv\"\nfor i, row in tqdm(sensors_train_erupt.iterrows(), total=sensors_train_erupt.shape[0]):\n    \n    sens_tmp = pd.read_csv(row['path'])\n    \n    sens_tmp = sens_tmp.rolling(window=100).mean()\n    sens_tmp = sens_tmp[0:-1:100]\n    sens_tmp['step'] = sens_tmp.index\n    \n    #sens_tmp.dropna(subset=sensor_cols, inplace=True, how='all') # drop row only if all sensors are missing\n    sens_tmp.reset_index(drop=True, inplace=True)\n    \n    sens_tmp['segment_id'] = row['segment_id']\n    sens_tmp['time_to_eruption_log'] = row['time_to_eruption_log']\n\n    \n    if i==0:\n        sens_tmp.to_csv(fname_train, index=False, mode='w', header=True)\n    else:\n        sens_tmp.to_csv(fname_train, index=False, mode='a', header=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sensors_test.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fname_test = \"test_low_freq.csv\"\nfor i, row in tqdm(sensors_test.iterrows(), total=sensors_test.shape[0]):\n    \n    sens_tmp = pd.read_csv(row['path'])\n    \n    sens_tmp = sens_tmp.rolling(window=100).mean()\n    sens_tmp = sens_tmp[0:-1:100]\n    sens_tmp['step'] = sens_tmp.index\n    \n    #sens_tmp.dropna(subset=sensor_cols, inplace=True, how='all') # drop row only if all sensors are missing\n    sens_tmp.reset_index(drop=True, inplace=True)\n    \n    sens_tmp['segment_id'] = row['segment_id']\n    \n    if i==0:\n        sens_tmp.to_csv(fname_test, index=False, mode='w', header=True)\n    else:\n        sens_tmp.to_csv(fname_test, index=False, mode='a', header=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_check = pd.read_csv(fname_train)\ntrain_check.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_check['segment_id'].nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_check = pd.read_csv(fname_test)\ntest_check.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_check['segment_id'].nunique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}