<jupyter_start><jupyter_text>DataInputs (X) consist of segments of waveform (i.e. a time series of amplitude measurements) that represents the vibrations of physical material, and is sampled at a rate of 4MHz.  Train data is one contiguous segment of ~630MM samples.  Test data are segments of length 150,000.Here we plot a sample of the original waveform viewed at successive levels of magnification (on the time axis).  The macro features (first graph) appear as aperiodic pulses of varying magnitude.  Within the pulses (fourth graph), there appears to be a fairly regular sinusoid oscillation at roughly 1 cycle every 16 samples, corresponding to a frequency of 4MHz / 16 = 250kHz.Outside of the pulses (last graph), the base pattern appears to be mostly noise.<jupyter_code>start = 130000000
t = np.arange(150000)

offset1, range1 = 60000, 15000
offset2, range2 = 63000, 1500
offset3, range3 = 63400, 150
noise_offset, noise_range = 63100, 150

fig = plt.figure(figsize=(18, 16))

ax = fig.add_subplot(5, 1, 1)
plt.plot(t, data[start:start+150000, 0])
plt.plot(t[offset1:offset1+range1], data[start+offset1:start+offset1+range1, 0], c='r')
ax.text(2, 120, 'Signal Segment, 150,000 Samples', fontsize=14)

ax = fig.add_subplot(5, 1, 2)
plt.plot(t[offset1:offset1+range1], data[start+offset1:start+offset1+range1, 0])
plt.plot(t[offset2:offset2+range2], data[start+offset2:start+offset2+range2, 0], c='r')
ax.text(offset1+2, 120, 'Zoom 10x', fontsize=14)

ax = fig.add_subplot(5, 1, 3)
plt.plot(t[offset2:offset2+range2], data[start+offset2:start+offset2+range2, 0])
plt.plot(t[offset3:offset3+range3], data[start+offset3:start+offset3+range3, 0], c='r')
plt.plot(t[noise_offset:noise_offset+noise_range],
         data[start+noise_offset:start+noise_offset+noise_range, 0], c='m')
ax.text(offset2+2, 120, 'Zoom 100x', fontsize=14)
ax.text(63130, 25, 'noise', fontsize=14)
ax.text(63410, 80, 'pulse', fontsize=14)

ax = fig.add_subplot(5, 1, 4)
plt.plot(t[offset3:offset3+range3], data[start+offset3:start+offset3+range3, 0], c='r')
ax.text(offset3+2, 95, 'Pulse Detail', fontsize=14)

ax = fig.add_subplot(5, 1, 5)
plt.plot(t[noise_offset:noise_offset+range3], data[start+noise_offset:start+noise_offset+range3, 0], c='m')
ax.text(noise_offset+2, 11, 'Noise Detail', fontsize=14)

for ax in fig.axes[:2]:
    ax.margins(x=0.003)
    ax.axis('off')

for ax in fig.axes[2:]:
    ax.margins(x=0.003)
    for s in ax.spines:
        ax.spines[s].set_visible(False);<jupyter_output><empty_output><jupyter_text>LabelsThe label/target is "time to failure" (TTF), which is the amount of time until the next quake.  The training data is unbalanced: samples with TTF > 11 are relatively scarce.<jupyter_code>fig = plt.figure(figsize=(16,4))
plt.plot(y_train_raw[::1000000]);<jupyter_output><empty_output><jupyter_text>DownsamplingUse scipy.signal.decimate to downsample the raw waveform.  This reduces the compute time to process spectrograms.  decimate also applies an anti-aliasing (i.e. low-pass) filter that reduces the noise present at the highest frequencies.<jupyter_code>def downsample(data, downsample_rate):
    trunc = int((data.shape[0] // downsample_rate) * downsample_rate)
    data = data[:trunc, :]
    x = signal.decimate(data[:, 0].ravel(), downsample_rate)
    y_ = data[:, 1].reshape((-1, downsample_rate)).T
    y = y_.mean(axis=0)
    return x, y

RAW_SEQ_LEN = 150000
DOWNSAMPLE_RATE = 4

assert (RAW_SEQ_LEN / DOWNSAMPLE_RATE) % 1 == 0
downsampled_seq_len = int(RAW_SEQ_LEN/DOWNSAMPLE_RATE)

X_train_raw, y_train_raw = downsample(data, DOWNSAMPLE_RATE)
del data

print(f'downsample rate: {DOWNSAMPLE_RATE}')
print(f'downsampled seq len: {downsampled_seq_len}')
print(f'X shape: {X_train_raw.shape} y shape: {y_train_raw.shape}')<jupyter_output>downsample rate: 4
downsampled seq len: 37500
X shape: (157286370,) y shape: (157286370,)<jupyter_text>SpectrogramsSpectrograms sample a signal at regular intervals and convert the waveform into power measurements on a frequency scale.  In other words, it is a time series of Fourier Transforms.  The scipy.signal.spectrogram function generates spectrograms from signal data.Parameters for the spectrogram function control the size of the sampling window, the type of convolution function ("window function") used for sampling, and other aspects of the FFT algorithm.  Several combinations of parameters were tested.  The logarithm of the power spectrum was used to generate features.<jupyter_code>def spectrogram(sig_in, dsamp):
    nperseg = 128 # default 256
    noverlap = nperseg // 2 # default: nperseg // 8
    fs = 4000000 // dsamp # sampling frequency.  raw signal sample rate is 4MHz
    window = 'hann'
    scaling = 'density' # {'density', 'spectrum'}
    detrend = 'linear' # {'linear', 'constant', False}
    eps = 1e-11  # epsilon
    f, t, Sxx = signal.spectrogram(sig_in, nperseg=nperseg, noverlap=noverlap,
                                   fs=fs, window=window,
                                   scaling=scaling, detrend=detrend)
    return f, t, np.log(Sxx + eps)<jupyter_output><empty_output><jupyter_text>Plot signal samples and their corresponding spectrograms:<jupyter_code>def plot_spec(sig_in):
    f, t, Sxx = spectrogram(sig_in, DOWNSAMPLE_RATE)
    print(f'spectrogram shape: {Sxx.shape}, size: {Sxx.size}')
    fig = plt.figure(figsize=(18, 8))
    ax = fig.add_subplot(2, 1, 1)
    ax.margins(x=0.003)
    plt.plot(sig_in)
    plt.title('downsampled signal:', fontsize=16, loc='left')
    plt.axis('off')
    ax = fig.add_subplot(2, 1, 2)
    cmap = plt.get_cmap('magma')
    spec = plt.pcolormesh(t, f, Sxx, cmap=cmap)
    plt.title('log spectrogram:',
              fontsize=16, loc='left')
    plt.axis('off');

start = rng.randint(0, X_train_raw.shape[0] - downsampled_seq_len)
sig_in = X_train_raw[start:start+downsampled_seq_len]
plot_spec(sig_in)
def get_Sxx_batch(source, n_samples, seq_len, dsamp_rate, mode=None):
    if source == 'train':
        X = X_train_; y = y_train_
    elif source == 'val':
        X = X_val_; y = y_val_
    else:
        X = source; y = [0]
    if source in ('train', 'val'):
        start_list = [rng.randint(0, X.shape[0]-seq_len)
                      for i in range(n_samples)]
    else:
        start_list = [0]
    X_batch = []
    for start in start_list:
        sig_in = X[start:start+seq_len]
        f, _, Sxx = spectrogram(sig_in, dsamp_rate)
        X_batch.append(Sxx)
    sxx_h, sxx_w = (X_batch[0].shape)
    X_batch = np.array(X_batch).reshape(n_samples, sxx_h, sxx_w)
    if source in ('train', 'val'):
        y_batch = [y[start+seq_len] for start in start_list]
    else:
        y_batch = [0]
    return X_batch, y_batch, f<jupyter_output><empty_output><jupyter_text>FeaturesFeatures generated from the spectrograms comprise two tiers of aggregations:Tier 1: Aggregations across all frequency bands: - Range (max - min) was a more powerful predictor than standard deviation. - Mean scored low in feature importance and was not included.  The tier-2 features were more powerful than tier-1 mean. - 10 and 90 percentiles, and the range/difference between the percentiles improved CV scores.Tier 2: Aggregations on individual frequency bands: - Mean and standard deviation - Quantiles/partitions within frequency band contributed most to accuracy.  Cut points were established by experimentation.  The most effective was top 4, then next 4, then next 8, 16, 32, etc.  In addition, the mean of the lowest 64 values were included. - Frequency-band quantiles in the 200kHz-300kHz range generated the most predictive power, corresponding to the resonant frequency range seen in the plots above.<jupyter_code>def get_feats_batch(args_in):
    source, batch_size = args_in
    X_out, y_out = [], []
    X_batch, y_batch, freqs = get_Sxx_batch(
        source, batch_size, downsampled_seq_len, DOWNSAMPLE_RATE)
    Sxx_w = X_batch.shape[2]
    
    freqs = ['%03d' % (i/1000) for i in freqs]
    
    # aggregations
    X_2d = X_batch.reshape((batch_size, -1))
    X_out.append(pd.DataFrame(X_2d.mean(axis=1), columns=['mean_']))
    X_out.append(pd.DataFrame(X_2d.std(axis=1), columns=['std_']))
    max_ = pd.DataFrame(X_2d.max(axis=1), columns=['max_'])
    min_ = pd.DataFrame(X_2d.min(axis=1), columns=['min_'])
    range_ = pd.DataFrame((max_.values-min_.values), columns=['range_'])
    X_out.extend([min_, max_, range_])

    # quantiles
    q90 = pd.DataFrame(np.quantile(X_2d, 0.9, axis=1), columns=['q90_'])
    q10 = pd.DataFrame(np.quantile(X_2d, 0.1, axis=1), columns=['q10_'])
    range80 = pd.DataFrame((q90.values-q10.values), columns=['range80_'])
    X_out.extend([q90, q10, range80])
    
    # aggregations by frequency band.
    X_out.append(pd.DataFrame(X_batch.mean(axis=2), columns=freqs).add_prefix('mean_'))
    X_out.append(pd.DataFrame(X_batch.std(axis=2), columns=freqs).add_prefix('std_'))
    
    # partition the power levels in each frequency band by quantiles,
    # then generate the mean of each quantile
    top_k = (0, 4, 8, 16, 32, 64, 128) # cut points for the quantiles
    cp_inv = [Sxx_w - i for i in top_k][::-1]
    part = np.partition(X_batch, cp_inv[:-1], axis=2)
    
    t1 = (pd.DataFrame(part[:, :, cp_inv[-2]:cp_inv[-1]].mean(axis=2),
                      columns=freqs)
         .add_prefix('top_' + str(top_k[1]) + '_mean_'))

    t2 = (pd.DataFrame(part[:, :, cp_inv[-3]:cp_inv[-2]].mean(axis=2),
                      columns=freqs)
         .add_prefix('top_' + str(top_k[2]) + '_mean_'))
    
    t3 = (pd.DataFrame(part[:, :, cp_inv[-4]:cp_inv[-3]].mean(axis=2),
                      columns=freqs)
         .add_prefix('top_' + str(top_k[3]) + '_mean_'))

    t4 = (pd.DataFrame(part[:, :, cp_inv[-5]:cp_inv[-4]].mean(axis=2),
                      columns=freqs)
         .add_prefix('top_' + str(top_k[4]) + '_mean_'))

    t5 = (pd.DataFrame(part[:, :, cp_inv[-6]:cp_inv[-5]].mean(axis=2),
                      columns=freqs)
         .add_prefix('top_' + str(top_k[5]) + '_mean_'))

    t6 = (pd.DataFrame(part[:, :, cp_inv[-7]:cp_inv[-6]].mean(axis=2),
                      columns=freqs)
         .add_prefix('top_' + str(top_k[6]) + '_mean_'))

    # a final partition is the mean of the lowest power levels in the spectrogram    
    low64_mean = pd.DataFrame((np.partition(X_batch, 64, axis=2)[:, :, :64]
                              .mean(axis=2)),
                             columns=freqs)
    low64_mean = low64_mean.add_prefix('low64_mean_')
    
    X_out.extend([low64_mean, t6, t5, t4, t3, t2, t1])
    
    X_out = pd.concat(X_out, axis=1)
    y_out = pd.Series(y_batch)
    return (X_out, y_out)<jupyter_output><empty_output><jupyter_text>Divide the training data into 5 folds, and generate random samples from each of the folds, for a total of 500,000 training features:<jupyter_code>def joblib_map(func, args, n_jobs=None, prefer='threads', verbose=10):
    from joblib import Parallel, delayed
    output = Parallel(n_jobs=n_jobs, prefer=prefer, verbose=verbose)(
    (delayed(func)(arg) for arg in args))
    return output

def get_feats_parallel(source, n_samples, batch_size=500, n_jobs=8):
    print('n_tasks =', n_samples // batch_size)
    targs = [(source, batch_size) for i in range(n_samples // batch_size)]
    feats = joblib_map(get_feats_batch, targs, n_jobs=n_jobs)
    X_out = pd.concat([batch[0] for batch in feats], axis=0)
    y_out = pd.concat([batch[1] for batch in feats])
    return X_out, y_out

N_FOLDS = 5
train_samples = 500000
interval = X_train_raw.shape[0] / N_FOLDS
splits = [int(interval * i) for i in range(N_FOLDS+1)]

X_folds, y_folds = [], []
for fold in range(N_FOLDS):
    X_train_ = X_train_raw[splits[fold]:splits[fold+1]]
    y_train_ = y_train_raw[splits[fold]:splits[fold+1]]
    print(f'building fold {fold+1} of {N_FOLDS}...')
    X_train, y_train = get_feats_parallel(
        'train', train_samples//N_FOLDS, n_jobs=8)
    X_train['fold'] = int(fold + 1)
    X_folds.append(X_train)
    y_folds.append(y_train)
X_train = pd.concat(X_folds, axis=0, ignore_index=True)
y_train = pd.concat(y_folds, ignore_index=True)
print('X_train.shape:', X_train.shape)<jupyter_output>building fold 1 of 5...
n_tasks = 200<jupyter_text>The resulting dataset comprises 593 features (last column is a fold designator):<jupyter_code>print(f'X_train size: {X_train.memory_usage().sum()/1000000000:.3f}GB')
print(f'shape {X_train.shape}')
X_train.head()<jupyter_output>X_train size: 2.376GB
shape (500000, 594)<jupyter_text>Re-weighting/Augmenting Training DataInitial CV results exhibit a fair amount of underfitting at the far ends of the labels, i.e. TTF values above 11 and below 1, most likely reflecting the imbalance of the training data.  To address the imbalance, we augment the training data by duplicating training samples outside the range of 1.0-11.0.  Augmentation/upsampling resulted in modest improvements in CV results.<jupyter_code>high_cutoff = 11
low_cutoff = 1
augment_rate = 1

def augment_train(X_train, y_train):
    high_idx = y_train[y_train>high_cutoff].index
    low_idx = y_train[y_train<low_cutoff].index
    print('high', len(high_idx), 'low', len(low_idx), 'total', len(high_idx) + len(low_idx))
    X_high = X_train.loc[high_idx, :]
    X_low = X_train.loc[low_idx, :]
    y_high = y_train.loc[high_idx]
    y_low = y_train.loc[low_idx]
    X_augment = [X_high] * augment_rate + [X_low] * augment_rate
    y_augment = [y_high] * augment_rate + [y_low] * augment_rate
    X_augment = pd.concat(X_augment, axis=0, ignore_index=True)
    y_augment = pd.concat(y_augment, ignore_index=True)
    return X_augment, y_augment

X_augment, y_augment = augment_train(X_train, y_train)
print('X_augment', X_augment.shape, 'y_augment', y_augment.shape)
X_train = pd.concat([X_train, X_augment], axis=0, ignore_index=True)
y_train = pd.concat([y_train, y_augment], ignore_index=True)
del X_augment; del y_augment
print(f'X_train size: {X_train.memory_usage().sum()/1000000000:.3f}GB, ' +
     f'X_train.shape: {X_train.shape}')<jupyter_output>high 47142 low 48887 total 96029
X_augment (96029, 594) y_augment (96029,)
X_train size: 2.832GB, X_train.shape: (596029, 594)<jupyter_text>Generate features from the test dataset:<jupyter_code>test_path = os.path.join(data_path, 'test')
files = os.listdir(test_path)
print('total files', len(files)) # 2624
preds = {'seg_id': [], 'time_to_failure': []}
feats_list = []
for fname in tqdm_notebook(files):
    path = os.path.join(test_path, fname)
    obs = np.array(pd.read_csv(path).values)
    data = np.c_[obs, np.zeros(len(obs))]
    x, y = downsample(data, DOWNSAMPLE_RATE)
    (feats, _) = get_feats_batch((x, 1))
    seg_id = fname.split('.')[0]
    feats['seg_id'] = seg_id
    feats_list.append(feats)
    
X_test = pd.concat(feats_list, axis=0).reset_index(drop=True)
import pickle
with open(os.path.join(DATA_PATH, OUTPUT_FILENAME), 'wb') as file:
    pickle.dump([X_train, y_train, X_test], file)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1316264456180226, 'fraction_numerical': 0.03924059094859282, 'mean_word_length': 4.238095238095238, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 39, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '193284', 'license': 'MIT', 'n_tokens_mistral': 5884, 'n_tokens_neox': 5431, 'n_words': 1310, 'path': 'features.ipynb', 'repo_name': 'epsdg/quake-prediction', 'size': '14959'}