{"cells":[{"metadata":{},"cell_type":"markdown","source":"Proces modelowania zaczniemy tak jak w analizie od wcześniejszego wczytania wszystkich danych treningowych i testowych, a następnie zdefiniujemy funkcje potrzebne do wyznaczania feature'ów \n\nDalej zgodnie z wnioskami z Analizy i Eksploracji Danych wybierzemy feature'y z korelacją > 0.4, a także obliczymy dla nich współczynnik Pearsona, gdzie wartości poniżej 0.05, zostaną dodane do modelu\n\nNie będziemy sprowadzać średniej do zera, pomimo iż wydawało się to słuszną koncepcją podczas analizy wstepne wyniki naszych modeli pokazały, że nie przynosiło to jednak żadnych rezultatów jeśli chodzi o sam proces predykcji czasu trzęsień"},{"metadata":{"trusted":true},"cell_type":"code","source":"import seaborn as sns\nimport matplotlib.pyplot as plt\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import GridSearchCV, KFold, train_test_split\nfrom sklearn.svm import SVR, NuSVR\nfrom sklearn.kernel_ridge import KernelRidge\nimport pandas as pd\nimport numpy as np\nimport os\nimport gc\nimport warnings\nwarnings.filterwarnings(\"ignore\")\n\nDATA_DIR = \"../input\"\nTEST_DIR = r'../input/test'\n\nld = os.listdir(TEST_DIR)\nsizes = np.zeros(len(ld))\n\nfrom scipy.signal import hilbert\nfrom scipy.signal import hann\nfrom scipy.signal import convolve\nfrom scipy.stats import pearsonr\nfrom scipy import stats\nfrom sklearn.kernel_ridge import KernelRidge\n\nimport lightgbm as lgb\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.metrics import mean_absolute_error\n\nfrom tsfresh.feature_extraction import feature_calculators\nfrom tqdm import tqdm\n\n%matplotlib inline\nsns.set_style('darkgrid')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def classic_sta_lta(x, length_sta, length_lta):\n    \n    sta = np.cumsum(x ** 2)\n\n    # Zamiana na float\n    sta = np.require(sta, dtype=np.float)\n\n    # Kopia dla LTA\n    lta = sta.copy()\n\n    # Obliczanie STA i LTA\n    sta[length_sta:] = sta[length_sta:] - sta[:-length_sta]\n    sta /= length_sta\n    lta[length_lta:] = lta[length_lta:] - lta[:-length_lta]\n    lta /= length_lta\n\n    # Uzupełnienie zerami\n    sta[:length_lta - 1] = 0\n\n    # Aby nie dzielić przez 0 ustawiamy 0 na małe liczby typu float\n    dtiny = np.finfo(0.0).tiny\n    idx = lta < dtiny\n    lta[idx] = dtiny\n\n    return sta / lta","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def calc_change_rate(x):\n    change = (np.diff(x) / x[:-1]).values\n    change = change[np.nonzero(change)[0]]\n    change = change[~np.isnan(change)]\n    change = change[change != -np.inf]\n    change = change[change != np.inf]\n    return np.mean(change)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"percentiles = [1, 5, 10, 20, 25, 30, 40, 50, 60, 70, 75, 80, 90, 95, 99]\nhann_windows = [50, 150, 1500, 15000]\nspans = [300, 3000, 30000, 50000]\nwindows = [10, 50, 100, 500, 1000, 10000]\nborders = list(range(-4000, 4001, 1000))\npeaks = [10, 20, 50, 100]\ncoefs = [1, 5, 10, 50, 100]\nlags = [10, 100, 1000, 10000]\nautocorr_lags = [5, 10, 50, 100, 500, 1000, 5000, 10000]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def gen_features(x, zero_mean=False):\n    if zero_mean==True:\n        x = x-x.mean()\n    strain = {}\n    strain['mean'] = x.mean()\n    strain['std']=x.std()\n    strain['max']=x.max()\n    strain['kurtosis']=x.kurtosis()\n    strain['skew']=x.skew()\n    zc = np.fft.fft(x)\n    realFFT = np.real(zc)\n    imagFFT = np.imag(zc)\n    strain['min']=x.min()\n    strain['sum']=x.sum()\n    strain['mad']=x.mad()\n    strain['median']=x.median()\n    \n    strain['mean_change_abs'] = np.mean(np.diff(x))\n    strain['mean_change_rate'] = np.mean(np.nonzero((np.diff(x) / x[:-1]))[0])\n    strain['abs_max'] = np.abs(x).max()\n    strain['abs_min'] = np.abs(x).min()\n    \n    strain['avg_first_50000'] = x[:50000].mean()\n    strain['avg_last_50000'] = x[-50000:].mean()\n    strain['avg_first_10000'] = x[:10000].mean()\n    strain['avg_last_10000'] = x[-10000:].mean()\n    \n    strain['min_first_50000'] = x[:50000].min()\n    strain['min_last_50000'] = x[-50000:].min()\n    strain['min_first_10000'] = x[:10000].min()\n    strain['min_last_10000'] = x[-10000:].min()\n    \n    strain['max_first_50000'] = x[:50000].max()\n    strain['max_last_50000'] = x[-50000:].max()\n    strain['max_first_10000'] = x[:10000].max()\n    strain['max_last_10000'] = x[-10000:].max()\n    \n    strain['max_to_min'] = x.max() / np.abs(x.min())\n    strain['max_to_min_diff'] = x.max() - np.abs(x.min())\n    strain['count_big'] = len(x[np.abs(x) > 500])\n           \n    strain['mean_change_rate_first_50000'] = calc_change_rate(x[:50000])\n    strain['mean_change_rate_last_50000'] = calc_change_rate(x[-50000:])\n    strain['mean_change_rate_first_10000'] = calc_change_rate(x[:10000])\n    strain['mean_change_rate_last_10000'] = calc_change_rate(x[-10000:])\n    \n    strain['q95'] = np.quantile(x, 0.95)\n    strain['q99'] = np.quantile(x, 0.99)\n    strain['q05'] = np.quantile(x, 0.05)\n    strain['q01'] = np.quantile(x, 0.01)\n    \n    strain['abs_q95'] = np.quantile(np.abs(x), 0.95)\n    strain['abs_q99'] = np.quantile(np.abs(x), 0.99)\n    strain['abs_q05'] = np.quantile(np.abs(x), 0.05)\n    strain['abs_q01'] = np.quantile(np.abs(x), 0.01)\n    \n    for autocorr_lag in autocorr_lags:\n        strain['autocorrelation_' + str(autocorr_lag)] = feature_calculators.autocorrelation(x, autocorr_lag)\n    \n    # percentiles on original and absolute values\n    for p in percentiles:\n        strain['percentile_'+str(p)] = np.percentile(x, p)\n        strain['abs_percentile_'+str(p)] = np.percentile(np.abs(x), p)\n    \n    strain['abs_mean'] = np.abs(x).mean()\n    strain['abs_std'] = np.abs(x).std()\n    \n    strain['quantile_0.95']=np.quantile(x, 0.95)\n    strain['quantile_0.99']=np.quantile(x, 0.99)\n    strain['quantile_0.05']=np.quantile(x, 0.05)\n    strain['realFFT_mean']=realFFT.mean()\n    strain['realFFT_std']=realFFT.std()\n    strain['realFFT_max']=realFFT.max()\n    strain['realFFT_min']=realFFT.min()\n    strain['imagFFT_mean']=imagFFT.mean()\n    strain['imagFFT_std']=realFFT.std()\n    strain['imagFFT_max']=realFFT.max()\n    strain['imaglFFT_min']=realFFT.min()\n    \n    strain['std_first_50000']=x[:50000].std()\n    strain['std_last_50000']=x[-50000:].std()\n    strain['std_first_25000']=x[:25000].std()\n    strain['std_last_25000']=x[-25000:].std()\n    strain['std_first_10000']=x[:10000].std()\n    strain['std_last_10000']=x[-10000:].std()\n    strain['std_first_5000']=x[:5000].std()\n    strain['std_last_5000']=x[-5000:].std()\n        \n    strain['Hilbert_mean'] = np.abs(hilbert(x)).mean()\n    strain['Hann_window_mean'] = (convolve(x, hann(150), mode='same') / sum(hann(150))).mean()\n    strain['classic_sta_lta1_mean'] = classic_sta_lta(x, 500, 10000).mean()\n    strain['classic_sta_lta2_mean'] = classic_sta_lta(x, 5000, 100000).mean()\n    strain['classic_sta_lta3_mean'] = classic_sta_lta(x, 3333, 6666).mean()\n    strain['classic_sta_lta4_mean'] = classic_sta_lta(x, 10000, 25000).mean()\n    strain['classic_sta_lta6_mean'] = classic_sta_lta(x, 100, 5000).mean()\n    strain['classic_sta_lta8_mean'] = classic_sta_lta(x, 4000, 10000).mean()\n    strain['Moving_average_700_mean'] = x.rolling(window=700).mean().mean(skipna=True)\n    moving_average_700_mean = x.rolling(window=700).mean().mean(skipna=True)\n    ewma = pd.Series.ewm\n    strain['exp_Moving_average_300_mean'] = (ewma(x, span=300).mean()).mean(skipna=True)\n    strain['exp_Moving_average_3000_mean'] = ewma(x, span=3000).mean().mean(skipna=True)\n    strain['exp_Moving_average_30000_mean'] = ewma(x, span=30000).mean().mean(skipna=True)\n    no_of_std = 3\n    strain['MA_700MA_std_mean'] = x.rolling(window=700).std().mean()\n    strain['MA_1000MA_std_mean'] = x.rolling(window=1000).std().mean()\n    \n    strain['iqr'] = np.subtract(*np.percentile(x, [75, 25]))\n    strain['q999'] = np.quantile(x,0.999)\n    strain['q001'] = np.quantile(x,0.001)\n    strain['ave10'] = stats.trim_mean(x, 0.1)\n        \n    for window in windows:\n        x_roll_std = x.rolling(window).std().dropna().values\n        x_roll_mean = x.rolling(window).mean().dropna().values\n        \n        strain['ave_roll_std_' + str(window)] = x_roll_std.mean()\n        strain['std_roll_std_' + str(window)] = x_roll_std.std()\n        strain['max_roll_std_' + str(window)] = x_roll_std.max()\n        strain['min_roll_std_' + str(window)] = x_roll_std.min()\n        strain['q01_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.01)\n        strain['q05_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.05)\n        strain['q95_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.95)\n        strain['q99_roll_std_' + str(window)] = np.quantile(x_roll_std, 0.99)\n        strain['av_change_abs_roll_std_' + str(window)] = np.mean(np.diff(x_roll_std))\n        strain['av_change_rate_roll_std_' + str(window)] = np.mean(np.nonzero((np.diff(x_roll_std) / x_roll_std[:-1]))[0])\n        strain['abs_max_roll_std_' + str(window)] = np.abs(x_roll_std).max()\n        \n        for p in percentiles:\n            strain['percentile_roll_std_' + str(p) + '_window_' + str(window)] = np.percentile(x_roll_std, p)\n            strain['percentile_roll_mean_' + str(p) + '_window_' + str(window)] = np.percentile(x_roll_mean, p)\n        \n        strain['ave_roll_mean_' + str(window)] = x_roll_mean.mean()\n        strain['std_roll_mean_' + str(window)] = x_roll_mean.std()\n        strain['max_roll_mean_' + str(window)] = x_roll_mean.max()\n        strain['min_roll_mean_' + str(window)] = x_roll_mean.min()\n        strain['q01_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.01)\n        strain['q05_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.05)\n        strain['q95_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.95)\n        strain['q99_roll_mean_' + str(window)] = np.quantile(x_roll_mean, 0.99)\n        strain['av_change_abs_roll_mean_' + str(window)] = np.mean(np.diff(x_roll_mean))\n        strain['av_change_rate_roll_mean_' + str(window)] = np.mean(np.nonzero((np.diff(x_roll_mean) / x_roll_mean[:-1]))[0])\n        strain['abs_max_roll_mean_' + str(window)] = np.abs(x_roll_mean).max()\n        \n        \n    return pd.Series(strain)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'), iterator=True, chunksize=150_000, dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32})\nX_train = pd.DataFrame()\ny_train = pd.Series()\n\nfor df in tqdm(train_df):\n    features = gen_features(df['acoustic_data'])\n    X_train = X_train.append(features, ignore_index=True)\n    y_train = y_train.append(pd.Series(df['time_to_failure'].values[-1]), ignore_index=True)\n\nX_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train_df\nX_test = pd.DataFrame()\n\nfor i, f in tqdm(enumerate(ld)):\n    df = pd.read_csv(os.path.join(TEST_DIR, f))\n    features = gen_features(df['acoustic_data'])\n    X_test = X_test.append(features, ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Postanowiliśmy więc zbadać korelacje między wszystkimi feature'ami a wartością zmiennej wyjściowej\n\nW tym kroku postanowiliśmy wyznaczyć poszczególne korelacje pomiędzy cechami tak, aby zobaczyć, które mają najlepsze wyniki. Pod uwagę braliśmy korelacje o wartości bezwzględnej większej lub równej 0.45. Jak widać w wyniku poniższego kawałku kodu takich korelacji mamy: X. Najlepsze z nich mają wartości powyżej 0.6, to na nich powinniśmy skupić swoją największą uwagę i to one powinny mieć największy wpływ na nasz model."},{"metadata":{"trusted":true},"cell_type":"code","source":"corelations = np.abs(X_train.corrwith(y_train)).sort_values(ascending=False)\ncorelations_df = pd.DataFrame(data=corelations, columns=['corr'])\nprint(\"Number of high corelated values: \",corelations_df[corelations_df['corr']>=0.55]['corr'].count())\n\nhigh_corr = corelations_df[corelations_df['corr']>=0.55]\nprint(high_corr)\nhigh_corr_labels = high_corr.reset_index()['index'].values\n#print(high_corr_labels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X_train_high_corr = X_train[high_corr_labels]\nX_test_high_corr = X_test[high_corr_labels]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Poniżej wykonujemy standaryzację danych. StandardScaler odpowiedzialny jest za przekształcenie naszych danych w taki sposób, że wartość mean będzie wynosiła 0, a odchylenie standardowe 1. Spowoduje to, że dla przekazanych danych od każdej wartości w zbiorze danych zostanie odjęta wartość średnia próbki, a następnie podzielona zostanie przez odchylenie standardowe całego zbioru danych. Dzięki temu nasz zbiór danych zostanie znormalizowany. "},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.preprocessing import MinMaxScaler\n\nscaler = StandardScaler()\nscaler.fit(X_train_high_corr)\nX_train_scaled = pd.DataFrame(scaler.transform(X_train_high_corr), columns=X_train_high_corr.columns)\nX_test_scaled = pd.DataFrame(scaler.transform(X_test_high_corr), columns=X_test_high_corr.columns)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Poniżej liczyliśmy dla przeskalowanego zbioru treningowego korelację Pearsona. Mierzy ona liniową zależność pomiędzy dwoma zbiorami danych. Korelację musimy liczyć na danych znormalizowanych - co też uczyniliśmy w poprzednich krokach, dodatkowo wymagane jest, aby przekazywane zbiory danych miały rozkład normalny. \n\nPodobnie jak w przypadku zwykłej korelacji, korelacja Pearson'a zwraca wartości między -1, a 1, lecz wartość 0 określa - brak korelacji pomiędzy danymi. Skrajne wartości reprezentują dokładną liniową zależność. \n\nJednak w tym wypadku zależy nam na tzw p-wartości zwracanej przez funkcję __pearsonr__. Wartość symbolizuje prawdopodobieństwo wyprodukowania nieskorelowanych zbiorów danych przez system, które mają korelację przynajmniej tak wysoką jak ta powstałą z tego zbioru danych. Wszystkie wartości z p-wartością < 0.05 zostają przez nas przyjęto do modelu"},{"metadata":{"trusted":true},"cell_type":"code","source":"p_columns = []\np_corr = []\np_values = []\n\nfor col in X_train_scaled.columns:\n    p_columns.append(col)\n    p_corr.append(abs(pearsonr(X_train_scaled[col], y_train.values)[0]))\n    p_values.append(abs(pearsonr(X_train_scaled[col], y_train.values)[1]))\n\ndf = pd.DataFrame(data={'column': p_columns, 'corr': p_corr, 'p_value': p_values}, index=range(len(p_columns)))\ndf.sort_values(by=['corr', 'p_value'], inplace=True)\ndf.dropna(inplace=True)\ndf = df.loc[df['p_value'] <= 0.05]\n\ndrop_cols = []\n\nfor col in X_train_scaled.columns:\n    if col not in df['column'].tolist():\n        drop_cols.append(col)\n\nprint(drop_cols)\nprint('--------------------')\nprint(X_train_high_corr.columns.values)\n        \nX_train_scaled = X_train_scaled.drop(labels=drop_cols, axis=1)\nX_test_scaled = X_test_scaled.drop(labels=drop_cols, axis=1)\n\nX_train_scaled_minmax = X_train_scaled_minmax.drop(labels=drop_cols, axis=1)\nX_test_scaled_minmax = X_test_scaled_minmax.drop(labels=drop_cols, axis=1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## DNN model"},{"metadata":{},"cell_type":"markdown","source":"Do stworzenia głębokiej sieci neuronowej użyliśmy frameworku Keras. Model Sequential pozwala tworzyć sieć warstwa po warstwie. Warstwy typu Dense oznaczają, że każdy neuron jest połączony do każdego neuronu w kolejnej warstwie. W przypadku regresji jako funkcję aktywacji w warstwach ukrytych użyliśmy relu - Rectified Linear Unit. Po dodaniu każdej wartwy występuje Dropout(), zapobiegający przetrenowaniu sieci. Działanie sprowadza się do losowego usunięcia połowy neutronów. \nNastępnie ustalamy sieci optymalizator, który kontroluje współczynnik uczenia."},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.callbacks import ModelCheckpoint\nfrom keras.models import Sequential\nfrom keras.layers import Dense, Activation, Flatten, Dropout\nfrom sklearn.model_selection import train_test_split\nNN_model = Sequential()\n\n# The Input Layer :\nNN_model.add(Dense(128, kernel_initializer='RandomUniform',input_dim = X_train_scaled.shape[1], activation='relu'))\nNN_model.add(Dropout(0.5))\n# The Hidden Layers :\nNN_model.add(Dense(256, kernel_initializer='RandomUniform',activation='relu'))\nNN_model.add(Dropout(0.5))\nNN_model.add(Dense(256, kernel_initializer='RandomUniform',activation='relu'))\nNN_model.add(Dropout(0.5))\nNN_model.add(Dense(128, kernel_initializer='RandomUniform',activation='relu'))\n# The Output Layer :\nNN_model.add(Dense(1, kernel_initializer='RandomUniform',activation='linear'))\n\n# Compile the network :\nNN_model.compile(loss='mean_absolute_error', optimizer='adam', metrics=['mean_absolute_error'])\nNN_model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"checkpoint_name = 'Weights-{epoch:03d}--{val_loss:.5f}.hdf5' \ncheckpoint = ModelCheckpoint(checkpoint_name, monitor='val_loss', verbose = 1, save_best_only = True, mode ='auto')\ncallbacks_list = [checkpoint]\nNN_model.fit(X_train_scaled, y_train, epochs=100, batch_size=32, validation_split = 0.2, callbacks=callbacks_list)\npredictions_DNN = NN_model.predict(X_test_scaled)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission_DNN = pd.read_csv(os.path.join(DATA_DIR, 'sample_submission.csv'), dtype={\n    'acoustic_data': np.int16, 'time_to_failure': np.float32})\nsubmission_DNN['time_to_failure'] = predictions_DNN\nsubmission_DNN.to_csv('result_DNN.csv', index=False)","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.7"}},"nbformat":4,"nbformat_minor":1}