<jupyter_start><jupyter_text>Component Wear Anomalies  OCME Vega Shrinker**Let's consider the Vega skinwrapper family of packaging machines by [OCME](https://www.ocme.com/en?_ptc=4)*** They work by wrapping products (bottles) in a _plastic film_* ...Which is _cut and heated_, so that the film shrinks and stabilizes the content  OCME Vega Shrinker**A public dataset about one of their machines is [publicly available from Kaggle](https://www.kaggle.com/inIT-OWL/one-year-industrial-component-degradation)*** The dataset contains a _run-to-failure experiment_* I.e. the machine was left running until one of it components became unserviceable  - Specificallly, it was the blade for cutting the film**This is an example of anomaly _due to component wear_*** It's a common type of anomaly* ...And run-to-failure experiments are a typical way to investigate them**All problems in this class share a few _properties_*** The behavior becomes _more and more distant_ from normal over time* There is a critical anomaly _at the end of the experiment_**We will try to tackle the problem using the techniques we know**  The Dataset**Let's have a first look at the dataset**<jupyter_code>print(f'Number of examples: {len(vs)}, number of inputs: {len(vs_in)}')
vs.head()<jupyter_output>Number of examples: 1062912, number of inputs: 8<jupyter_text>* There aren't many columns, but there are _many examples_!* The data refers to different measurement intervals (or "segments")* Each segment contains data sampled _every 4ms_  The Dataset**Let's check some statistics**<jupyter_code>vs.describe()<jupyter_output><empty_output><jupyter_text>* The data is neither normalized nor standardized  The Dataset**Let's check for missing values**<jupyter_code>vs[vs_in].isnull().any()<jupyter_output><empty_output><jupyter_text>* There are none  The Dataset**And let's check the length of each segment**<jupyter_code>vs.groupby('segment')['mode'].count().describe()<jupyter_output><empty_output><jupyter_text>* There are 519 segments overall* ...Each with 2048 samples  The Dataset**The machine has multiple operating modes**<jupyter_code>nn.plot_series(vs['mode'], figsize=figsize)<jupyter_output><empty_output><jupyter_text>* The mode is _a controlled parameter_ and does not change in the middle of a segment* Intuitively, the mode _has an impact_ on the machine behavior  The Dataset**Let's have a look at all the sensor readings**<jupyter_code>vss = vs.copy()
vss[vs_in] = (vss[vs_in] - vss[vs_in].mean()) / vss[vs_in].std()
nn.plot_dataframe(vss[vs_in], figsize=figsize)<jupyter_output><empty_output><jupyter_text>Data Preparation  Binning**This dataset contain _high-frequency data_ (4ms sampling period)*** In this situation, feeding the raw data to a model does not usually make sense* Instead, we reduce the frequency of the data via a process called _binning_**A binning approach typically works as follows**We apply a sliding window, but so that its consecutive applications _do not overlap_* Each window application is called a _bin_* ...From which we extract one or more _features_**The result is series that contains a _smaller number of samples_**...But typically a _larger number of features_  Binning**There are two broad classes of features that are usually extracted*** _Time-domain_ features (e.g. mean, standard deviation)* _Frequency-domain_ features (e.g. specific FFT amplitudes)For this case, we will focus on time-domain features**As a first step, we defined _which features_ we are going to extract**<jupyter_code>functions = ['mean', 'std', 'skew', lambda x: x.kurtosis()]
aggmap = {a: functions for a in vs_in}
aggmap['mode'] = 'first'
str(aggmap)<jupyter_output><empty_output><jupyter_text>Binning**Then we define bin numbers, and extract the features via a `groupby` operation**We need to take care so that _no bin crosses between different segments_<jupyter_code>%%time
binsize = 128
bins = []
for sname, sdata in vs.groupby('segment'):
    # Build the bin numbers
    sdata['bin'] = sdata.index // binsize
    # Apply the aggregation functions
    tmp = sdata.groupby('bin').agg(aggmap)
    bins.append(tmp)
vsb = pd.concat(bins)<jupyter_output>CPU times: user 32 s, sys: 23.9 ms, total: 32 s
Wall time: 32.1 s<jupyter_text>* This can be a relatively slow operation* Bin numbers are usually easy to define using positional indexes and an integer division  Binning**If we choose the bin size correctly, we can speed up the operation*** In particular, if all segments have the same length...* ...And we choose bin size that is a submultiple of the segment length...Then we can avoid processing each segment separately:<jupyter_code>%%time
# Build the bin numbers
binsize = 128
vsb = vs.copy()
vsb['bin'] = vs.index // binsize
vsb = vsb.groupby('bin').agg(aggmap)<jupyter_output>CPU times: user 16.5 s, sys: 46 ms, total: 16.5 s
Wall time: 16.5 s<jupyter_text>* This kind of approach is significantly faster  Binning**The resulting dataframe has a hierarchical column index**<jupyter_code>vsb.iloc[:1]<jupyter_output><empty_output><jupyter_text>**It may be worth _flattening_, so as to simplify access:**<jupyter_code>if isinstance(vsb.columns, pd.MultiIndex):
    vsb.columns = ['::'.join(c) for c in vsb.columns]
vsb.iloc[:1]<jupyter_output><empty_output><jupyter_text>Standardization**Before we can train any model, we need some preparation**We will _standardize sensor inputs_ (all except the mode) using the first third of the series<jupyter_code>sep = int(np.round(len(vsb) * 0.34))
vsb_in = vsb.columns[:-1]
vsbs = vsb.copy()
tmp = vsbs[vsb_in].iloc[:sep]
vsbs[vsb_in] = (vsbs[vsb_in] - tmp.mean()) / tmp.std()
vsbs.iloc[:3]<jupyter_output><empty_output><jupyter_text>Categorical Mode**We will also adopt a _categorical encoding_ for the operating mode**This is critical for neural network approaches in particular<jupyter_code>from tensorflow.keras.utils import to_categorical
cmode = to_categorical(vsbs['mode::first'])
cols = [f'm{i}' for i in range(cmode.shape[1])]
cmode = pd.DataFrame(index=vsbs.index, data=cmode, columns=cols)
vsbs[cols] = cmode
vsbs.head()<jupyter_output><empty_output><jupyter_text>KDE Approach  KDE Approach**Now, let's try anomaly detection via KDE**First we estimate the optimal bandwidth:<jupyter_code>%%time
vsbs_tr = vsbs.iloc[:sep]

params = {'bandwidth': np.linspace(0.2, 0.8, 10)}
opt = GridSearchCV(KernelDensity(kernel='gaussian'), params, cv=5)
opt.fit(vsbs_tr)
best_params = pd.Series(index=opt.best_params_.keys(), data=opt.best_params_.values())
print(best_params)<jupyter_output>bandwidth    0.333333
dtype: float64
CPU times: user 8.4 s, sys: 0 ns, total: 8.4 s
Wall time: 8.41 s<jupyter_text>KDE Approach**Then we can train an estimator**<jupyter_code>h = opt.best_params_['bandwidth']
kde = KernelDensity(bandwidth=h)
kde.fit(vsbs_tr)<jupyter_output><empty_output><jupyter_text>**...And we can generate the alarm signal**<jupyter_code>%%time
ldens = kde.score_samples(vsbs)
signal_kde = pd.Series(index=vsbs.index, data=-ldens)<jupyter_output>CPU times: user 3.24 s, sys: 6.64 ms, total: 3.25 s
Wall time: 3.25 s<jupyter_text>KDE Approach**Let's plot the signal**<jupyter_code>nn.plot_signal(signal_kde, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* A peak at mid run makes the signal difficult to read  KDE Approach**Since the anomaly arises slowly, it makes sense to _clip_ and _smooth_ the signal**<jupyter_code>nn.plot_signal(signal_kde.clip(upper=25).rolling(512).mean(), figsize=figsize)<jupyter_output><empty_output><jupyter_text>* The smoothed signal has a tendency to grow* There is a plateau in the middle that is difficult to explain  Autoencoder Approach  Autoencoder Approach**Let's repeat our analysis using an autoencoder**First we define its structure:<jupyter_code>input_shape = vsbs.shape[1]
output_shape = len(vsb_in)
ae_x = keras.Input(shape=input_shape, dtype='float32')
ae_z = Dense(16, activation='relu')(ae_x)
ae_y = Dense(output_shape, activation='linear')(ae_z)
ae = keras.Model(ae_x, ae_y)<jupyter_output><empty_output><jupyter_text>* The input includes the operating mode* ...But _the output does not_!We have _no interest in reconstructing controlled parameters_  Autoencoder Approach**Now we can perform training**<jupyter_code>%%time
ae.compile(optimizer='RMSProp', loss='mse')
callbacks = [EarlyStopping(monitor='val_loss', patience=20, restore_best_weights=True)]
history = ae.fit(vsbs_tr, vsbs_tr[vsb_in], validation_split=0.2, callbacks=callbacks,
                 batch_size=32, epochs=100, verbose=2)<jupyter_output>Epoch 1/100
71/71 - 0s - loss: 1.1143 - val_loss: 0.7524
Epoch 2/100
71/71 - 0s - loss: 0.8494 - val_loss: 0.5817
Epoch 3/100
71/71 - 0s - loss: 0.6537 - val_loss: 0.4514
Epoch 4/100
71/71 - 0s - loss: 0.5175 - val_loss: 0.3687
Epoch 5/100
71/71 - 0s - loss: 0.4367 - val_loss: 0.3198
Epoch 6/100
71/71 - 0s - loss: 0.3797 - val_loss: 0.2831
Epoch 7/100
71/71 - 0s - loss: 0.3361 - val_loss: 0.2562
Epoch 8/100
71/71 - 0s - loss: 0.3025 - val_loss: 0.2329
Epoch 9/100
71/71 - 0s - loss: 0.2757 - val_loss: 0.2133
Epoch 10/100
71/71 - 0s - loss: 0.2545 - val_loss: 0.1975
Epoch 11/100
71/71 - 0s - loss: 0.2369 - val_loss: 0.1827
Epoch 12/100
71/71 - 0s - loss: 0.2218 - val_loss: 0.1687
Epoch 13/100
71/71 - 0s - loss: 0.2086 - val_loss: 0.1577
Epoch 14/100
71/71 - 0s - loss: 0.1966 - val_loss: 0.1481
Epoch 15/100
71/71 - 0s - loss: 0.1861 - val_loss: 0.1376
Epoch 16/100
71/71 - 0s - loss: 0.1774 - val_loss: 0.1295
Epoch 17/100
71/71 - 0s - loss: 0.1691 - val_loss: 0.1244
Epoch 18/100
71/71 - 0s[...]<jupyter_text>Autoencoder Approach**Let's see the loss evolution over time**<jupyter_code>nn.plot_training_history(history, figsize=figsize)<jupyter_output><empty_output><jupyter_text>Autoencoder Approach**Then we obtain our predictions**<jupyter_code>%%time
preds = ae.predict(vsbs)
preds = pd.DataFrame(index=vsbs.index, columns=vsb_in, data=preds)<jupyter_output>CPU times: user 182 ms, sys: 10.5 ms, total: 192 ms
Wall time: 157 ms<jupyter_text>...And we generate the alarm signal:<jupyter_code>se = np.square(preds - vsbs[vsb_in])
sse = np.sum(se, axis=1)
signal_ae = pd.Series(index=vsbs.index, data=sse)<jupyter_output><empty_output><jupyter_text>* We stored also the individual errors `se` for a later analysis  Autoencoder Approach**We can now plot the alarm signal**<jupyter_code>nn.plot_signal(signal_ae, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* Once again, the signal is difficult to read due to the peak  Autoencoder Approach**Let's apply clipping and smoothing**<jupyter_code>nn.plot_signal(signal_ae.clip(upper=10).rolling(512).mean(), figsize=figsize)<jupyter_output><empty_output><jupyter_text>* Again, a tendency to grow (weaker than before)* ...But this time _we can try  to explain the peak_!  Mutiple Signal Analysis**Let's investigate the situation**<jupyter_code>signals_ae = pd.DataFrame(index=vsbs.index, columns=vsb_in, data=se)
nn.plot_dataframe(signals_ae, vmin=-1, vmax=1, figsize=figsize)<jupyter_output><empty_output><jupyter_text>* As expected, errors are concentrated on a few features  Multiple Signal Analysis**Let's focus on the peak at mid run and check the largest errors**<jupyter_code>tmp = se.iloc[3500:4000].mean().sort_values(ascending=False)[:10]
nn.plot_bars(tmp, figsize=figsize, tick_gap=-1)
print(f'The largest error is on {tmp.index[0]}')<jupyter_output>The largest error is on pSvolFilm::CTRL_Position_controller::Actual_position::<lambda_0><jupyter_text>Multiple Signal Analysis**Let's what was going on with the original series**<jupyter_code>nn.plot_series(vs['pCut::CTRL_Position_controller::Actual_position'], figsize=figsize)<jupyter_output><empty_output><jupyter_text>* Indeed, _there is an unusual oscillation_! A domain expert may make sense of that  RNVP Approach  RNVP Approach**Let's make an attempt with Real NVP**<jupyter_code>%%time
vs_rnvp = nn.RealNVP(input_shape=vsbs.shape[1],
                num_coupling=6, units_coupling=32, depth_coupling=1, reg_coupling=0.01)
vs_rnvp.compile(optimizer='Adam')
X = vsbs_tr.astype(np.float32).values
cb = [EarlyStopping(monitor='loss', patience=10, min_delta=0.001, restore_best_weights=True)]
history = vs_rnvp.fit(X, batch_size=256, epochs=300, verbose=1, callbacks=cb)<jupyter_output>Epoch 1/300
12/12 [==============================] - 2s 3ms/step - loss: 207.9114
Epoch 2/300
12/12 [==============================] - 0s 3ms/step - loss: 92.4916
Epoch 3/300
12/12 [==============================] - 0s 4ms/step - loss: 75.9836
Epoch 4/300
12/12 [==============================] - 0s 4ms/step - loss: 70.8679
Epoch 5/300
12/12 [==============================] - 0s 3ms/step - loss: 69.0279
Epoch 6/300
12/12 [==============================] - 0s 3ms/step - loss: 65.2163
Epoch 7/300
12/12 [==============================] - 0s 3ms/step - loss: 62.8103
Epoch 8/300
12/12 [==============================] - 0s 3ms/step - loss: 63.1269
Epoch 9/300
12/12 [==============================] - 0s 3ms/step - loss: 58.0721
Epoch 10/300
12/12 [==============================] - 0s 3ms/step - loss: 55.8594
Epoch 11/300
12/12 [==============================] - 0s 3ms/step - loss: 54.1732
Epoch 12/300
12/12 [==============================] - 0s 3ms/step - loss: 52.3561
Epoch 13/300
12/12 [====[...]<jupyter_text>RNVP Approach**Let's have a look at the loss evolution over time**<jupyter_code>nn.plot_training_history(history, figsize=figsize)<jupyter_output><empty_output><jupyter_text>RNVP Approach**Now we can generate and plot the signal**<jupyter_code>ldens = vs_rnvp.score_samples(vsbs.astype(np.float32).values)
signal_vs = pd.Series(index=vsbs.index, data=-ldens)
nn.plot_signal(signal_vs, figsize=figsize)<jupyter_output><empty_output><jupyter_text>RNVP Approach**Finally, we apply clipping and smoothing**<jupyter_code>nn.plot_signal(signal_vs.clip(upper=10).rolling(512).mean(), figsize=figsize)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14247660187185024, 'fraction_numerical': 0.05169186465082793, 'mean_word_length': 5.900645802285147, 'pattern_counts': {'":': 0, '<': 121, '<?xml version=': 0, '>': 121, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '59331', 'license': 'MIT', 'n_tokens_mistral': 5337, 'n_tokens_neox': 4717, 'n_words': 1363, 'path': 'notebooks/3. Component Wear Anomalies.ipynb', 'repo_name': 'lompabo/aiiti-course-2021-04', 'size': '13890'}