DEV Community: yuto

【Parallel processing in Python】Joblib explained

yuto — Sat, 29 Jun 2024 07:21:13 +0000

1. What is Joblib

Joblib is a Python library that provides tools for efficiently saving and loading Python objects, particularly useful for machine learning workflows.

・Install

pip install joblib

2. Key Features

2.1 Parallel Processing

Joblib provides easy-to-use parallel processing capabilities through its Parallel and delayed functions. This is useful for tasks that can be parallelized, such as parameter grid searches or data preprocessing.

from joblib import Parallel, delayed

def process_data(data):
    # Simulate a time-consuming data processing step
    import time
    time.sleep(1)
    return data ** 2

data = [1, 2, 3, 4, 5]

# Parallel(n_jobs=workjob_num)(delayed(func_be_applied)(aug) for elem in array
results = Parallel(n_jobs=2)(delayed(process_data)(d) for d in data)
print(results)

We can use it smiply with list comprehensions as like above. If you specify n_jobs=-1, all available CPU cores will be used for the parallel computation. This can significantly speed up processing time for tasks that are CPU-bound and can be effectively parallelized.
However, it may affect other application using CPU or memory, so should be careful this.

:::details Speed test
・Test

from joblib import Parallel, delayed
import time

def process_data(data):
    # Simulate a time-consuming data processing step
    time.sleep(1)
    return data ** 2

data = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

# Normal calculation
start_time = time.time()
results_normal = [process_data(d) for d in data]
end_time = time.time()
normal_duration = end_time - start_time
print("Normal Calculation Results:", results_normal)
print("Normal Calculation Duration:", normal_duration, "seconds")

# Parallel calculation  with n_jobs=2
start_time = time.time()
results_parallel = Parallel(n_jobs=2)(delayed(process_data)(d) for d in data)
end_time = time.time()
parallel_duration = end_time - start_time
print("Parallel Calculation Results:", results_parallel)
print("Parallel Calculation Duration:", parallel_duration, "seconds")

# Parallel calculation with n_jobs=-1
start_time = time.time()
results_parallel = Parallel(n_jobs=-1)(delayed(process_data)(d) for d in data)
end_time = time.time()
parallel_duration = end_time - start_time
print("Parallel Calculation Results:", results_parallel)
print("Parallel Calculation Duration:", parallel_duration, "seconds")

・Result

# Normal Calculation Results: [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
# Normal Calculation Duration: 10.011737823486328 seconds
# Parallel Calculation Results: [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
# Parallel Calculation Duration: 5.565693616867065 seconds
# Parallel Calculation Results: [1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
# Parallel Calculation Duration: 3.627182722091675 seconds

:::

As you can see the results of test, parallel processing provides 2x or more faster calculation.

2.2 Serialization/Compression

Joblib using binary format when saves and loads python objects to disk. It provide efficient and faster operation.
Also, it supports various compression methods like zlib, gzip, bz2, and xz, allowing you to reduce the storage size of saved objects.

・Serialization

import joblib

data = [i for i in range(1000000)]

compression = False
if compression:
    joblib.dump(data, 'data.pkl', compress=('gzip', 3))
else:
    joblib.dump(data, 'data.pkl')

data = joblib.load('data.pkl')
print(len(data))
# 1000000

the '3' specifies the compression level (typically from 1 to 9, where higher numbers indicate more compression but slower speeds).

2.3 Memory Mapping

For large NumPy arrays, Joblib can use memory mapping to save memory by keeping a reference to the data on disk instead of loading it all into memory.

When you memory-map a file, parts of the file are loaded into RAM as needed, which can result in slower access times compared to having the entire dataset in RAM. However, it allows you to handle datasets larger than your available RAM.

If you wanna use the data, you have to load the data from storage disk.

from joblib import Memory
import math

cachedir = "./memory_cache"
memory = Memory(cachedir, verbose=0)

@memory.cache
def calc(x):
    print("RUNNING......")
    return math.sqrt(x)

print(calc(2))
print(calc(2))
print(calc(5))

# RUNNING......
# 1.41421356237
# 1.41421356237
# RUNNING......
# 2.23606797749979

As shown, the same calculation's result is returned without calculation(not through the func).
This is useful when doing same calculation, like Fibonacci sequence.

3. Summary

Joblib is so useful liblary in python. Espacially, parallel processing is crucial impact for like data preprocessing that must be faster.

Reference

[1] Joblib
[2] Joblibの様々な便利機能を把握する

【Kaggle】HMS 2nd Solution explaind

yuto — Sat, 29 Jun 2024 07:19:57 +0000

This time, I'll explain the 2nd solution of HMS competition on kaggle.
the original solution is here:
https://www.kaggle.com/competitions/hms-harmful-brain-activity-classification/discussion/492254

0. Competiton explanation

This competition's purpose is classification to the 182 type bird calls.
The validation method is AUC-ROC, so we had to predict all of birds classes.

Competition host provide the random length audio calling data over 240000 samples as training data, and data not labeled for utilize to augmentation or otehrs.

As the problem must be addressed were:
unbalanced data, probably noisy and mixed test data, data while bird not calling, etc..

Additionaly point to be noted is the limitation of submittion. This competition only allow CPU and under 120 minutes for inference. We have to both infer and process the data in time.

If you wanna know details, reccoemnd you check the competition overview and dataset.
https://www.kaggle.com/c/birdclef-2024

0.1 Technique used

Use MNE-tool for signal(eeg) processing
Use scipy.signal for filtering(scipy has various fiter)
hgnet: is this?

0.2 Point to see

2D-CNN is not well-equipped to capture positional information with in the channels. Because there is not pad in channel direction.
Using 3D-CNN model with 16channnel input fitered by 0.5~20Hz
Same model but different library's(mne and scipy) filter for diversity.
Double feature head(eeg + spectrum) (chap 1.4)
2 stage learning (chap 3)
Random offset sampling and average as label (chap 3)

After here is the solution.

1. Model

In early time, feeding a bsx4xHxM input to the 2D-CNN produce worse result than the one-image method. Then author start to think why?

Due to the position-sensitive nature of labels(LPD, GPD, etc.), We should take care of the channel dimention. 2D-CNN is not well-equipped to capture positional information with in the channels. Because there is not pad in channel direction. And that's why they need to concat the spectrum into one image for a vision model other than a bsx16xHxW image(double banana montage: spectrum created from 16 types of signals depending on the location of the brain).

So author decided to use 3d-CNN for spectrum, and 2D-CNN model for raw-eeg signal.

total solution like this:

mne 0.5-20Hz means use MNE-tool do filter.
scypy.signal means use scipy.signal to do filter.
For the reshape operator and the stft params please refer to the codes below.
And the number in () means final weights of the ensemble.[0.2, 0.2, 0.2, 0.1, 0.1, 0.2]

1.1 ×3d-l (spectrum model)

After double banana montage, +- 1024 clip and 0.5-20Hz filter, use stft to extract the spectrum, then feed to a 3d-CNN(x3d-l).

explanation of words(My interpretation):
・Double banana montage: This is a specific montage setup used in EEG signal processing where signals from pairs of electrodes are subtracted from each other to emphasize differences and reduce common noise.
・±1024 clip: Limits the signal amplitude to the range between -1024 and +1024.
・0.5-20Hz filter: Applies a bandpass filter to retain frequencies between 0.5 Hz and 20 Hz, which is typical for isolating certain EEG signal components.
・STFT (Short-Time Fourier Transform): Transforms the time-domain signal into the frequency domain over short time segments to extract the spectrum.
・3D-CNN (x3d-l): A 3D Convolutional Neural Network model that processes the spectral data.

The input data is a 16-channls spectrum image.
X3d-l cv 0.21+, public lb 0.25, private lb 0.29.

stft code like below, use it as a nn.Module:

class Transform50s(nn.Module):
    def __init__(self, ):
        super().__init__()
        self.wave_transform = torchaudio.transforms.Spectrogram(n_fft=512, win_length=128, hop_length=50, power=1)

    def forward(self, x):
        image = self.wave_transform(x)
        image = torch.clip(image, min=0, max=10000) / 1000
        n, c, h, w = image.size()
        image = image[:, :, :int(20 / 100 * h + 10), :]
        return image


class Transform10s(nn.Module):
    def __init__(self, ):
        super().__init__()
        self.wave_transform = torchaudio.transforms.Spectrogram(n_fft=512, win_length=128, hop_length=10, power=1)

    def forward(self, x):
        image = self.wave_transform(x)
        image = torch.clip(image, min=0, max=10000) / 1000
        n, c, h, w = image.size()
        image = image[:, :, :int(20 / 100 * h + 10), :]
        return image


class Model(nn.Module):
    def __init__(self):
        super().__init__()

        model_name = "x3d_l"
        self.net = torch.hub.load('facebookresearch/pytorchvideo',
                                  model_name, pretrained=True)

        self.net.blocks[5].pool.pool = nn.AdaptiveAvgPool3d(1)
        # self.net.blocks[5]=nn.Identity()
        # self.net.avgpool = nn.Identity()
        self.net.blocks[5].dropout = nn.Identity()
        self.net.blocks[5].proj = nn.Identity()
        self.net.blocks[5].activation = nn.Identity()
        self.net.blocks[5].output_pool = nn.Identity()

    def forward(self, x):
        x = self.net(x)
        return x


class Net(nn.Module):
    def __init__(self, num_classes=1):
        super().__init__()

        self.preprocess50s = Transform50s()
        self.preprocess10s = Transform10s()

        self.model = Model()

        self.pool = nn.AdaptiveAvgPool3d(1)
        self.fc = nn.Linear(2048, 6, bias=True)
        self.drop = nn.Dropout(0.5)

    def forward(self, eeg):
        # do preprocess
        bs = eeg.size(0)
        eeg_50s = eeg
        eeg_10s = eeg[:, :, 4000:6000]
        x_50 = self.preprocess50s(eeg_50s)
        x_10 = self.preprocess10s(eeg_10s)
        x = torch.cat([x_10, x_50], dim=1)
        x = torch.unsqueeze(x, dim=1)
        x = torch.cat([x, x, x], dim=1)
        x = self.model(x)
        # x = self.pool(x)
        x = x.view(bs, -1)
        x = self.drop(x)
        x = self.fc(x)
        return x

1.2 one image

It's a 2d vision model, efficientnetb5, public lb 0.262490 prvate lb 0.304877.
For 2d model, author concat 16 channels spectrum image like many peple does.

image = torch.reshape(image, shape=[n, 2, -1, w])
x1 = image[:, 0:1, ...]
x2 = image[:, 1:2, ...]
image = torch.cat([x1, x2], dim=-1)
image = torch.cat([image, image, image], dim=1)

Author check the submission score, combine these two spectrum model, get 0.28 private lb.

1.3 eeg model

View eeg (bsx16x10000) as an image. So expand dim=1(bsx1x16x10000), but the time dimention is too large. Then author make a reshape.

The model define as below:

rclass Net(nn.Module):
    def __init__(self,):
        super(Net, self).__init__()
        self.model = timm.create_model('efficientnet_b5', pretrained=True, in_chans=3)
        self.pool = nn.AdaptiveAvgPool2d(1)
        self.fc = nn.Linear(2048, out_features=6, bias=True)
        self.dropout = nn.Dropout(p=0.5)

    def extract_features(self, x):
        feature1 = self.model.forward_features(x)
        return feature1

    def forward(self, x):
        # extarct batch size
        bs = x.size(0)
        # reshape (bs, 16, 10000) -> (bs, 16, 1000, 10)
        reshaped_tensor = x.view(bs, 16, 1000, 10)
        # pemute to (bs, 16, 10, 1000)
        reshaped_and_permuted_tensor = reshaped_tensor.permute(0, 1, 3, 2)
        # reshape (bs, 16, 10, 1000) -> (bs, 160, 1000)
        reshaped_and_permuted_tensor = reshaped_and_permuted_tensor.reshape(bs, 16 * 10, 1000)
        # add new dimention at index 1, shape will be (bs, 1, 160, 1000)
        x = torch.unsqueeze(reshaped_and_permuted_tensor, dim=1)

        # concat the image to simulate RBG channel and increasing capacity, shape is (bs, 3, 160, 1000)
        x = torch.cat([x, x, x], dim=1)
        bs = x.size(0)

        x = self.extract_features(x)
        x = self.pool(x)
        x = x.view(bs, -1)
        x = self.dropout(x)
        x = self.fc(x)
        return x

forward_features() function in timm can extarct feature as like CNN, and it is used here.

The eeg 'reshaped' to bsx3x160x1000. With efficientnetb5, actives public lb 0.230978, private lb 0.282873.

There are 3models, but slightly different, mne-filter-efficientnetb5, scipy-filter-efficientnetb5 and mne-filter-hgnetb5

And the scires are slightly different,also.

1.4 Doublehead (eeg+spectrum)

Author use x3d-l to extarct the spectrum feature (with transform50s only), efficientnetb5 to extract the raw eeg feature.

Like the solution figure, concat the last feature, public lb 0.24, private lb 0.29? not sure.

2. Preprocess

2.1 Double banana montage, eeg as 16x10000
2.2 Filter with 0.5-20Hz
2.3 Clip by +-1024

3. Train

3.1 Stage 1, 15 epochs, with loss weght voters_num/20, Adamw lr=0.001, cos schedule
3.2 Stage 2, 5epoch, loss weight=1, voters_num>=6, Adamw lr=0.0001, cos schedule
3.3 use eeg_label_offset_seconds. Author random choose an offset for each eegid, and the target is average according to eegid for each train iter.
3.4 Augmentation: mirror eeg, flip between left brain data and right brain data
3.5 10 folds, then move 1000 samples from val to train, left 709 samples in val set. And use vote_num>=6 to do validation.

* voters_num equal number of vote, the competion data has number of voter for determine what the symptoms are.

4. Ensemble

By combining these models, author thinks it could get their current score. However their score is 6 model emsemble, but not improvement that much (0.28->0.27 private lb).

Final ensemble including:

2 spectrum model(x3d, efficientnetb5)
3 raw eeg model(efficientnetb5 with mne.filter, efficientnetb5 -butter filter, 1 ghnetb5 mne.filter)
1 eeg-spectrum mix model, butter filter

With weigts = [0.1,0.1,0.2,0.2,0.2,0.2].

ps. with-mne.filter means use mne lib to do filter, butter filter means use scipy.signal, just to add some diversity.

5 Some thought

Author thinks raw eeg is more important in this task. Feeding raw EEG data into a 2D visual model is somewhat similar to how humans observe EEG signals. Observe in time and channel dimensions！

6. Summary

Author using

raw eeg as 160x1000 data
spectrum as input of 2d-cnn and 3d-cnn
extract feature from both of eeg(efficientnetb5) and spectrum(x3d-l) for mix model
2stage train(only high quality data in stage 2), augmentation
different filter and impotant others. This is a great work and solution!

I think top-tier kaggler do all of basically things, and trying new methods with solid evidence of some experimental results. This is the important thing to win the competiton.

But, also enjoying competition is one of the most important thing, let we don't forget that!

Reference

[1] coolz, 2nd place solution, kaggle

What is CWT/Scarogram?

yuto — Wed, 17 Apr 2024 13:18:42 +0000

1. CWT/Scarogram

Scarogram is a 2-dimentional image created by applying CWT(Continous Wavelet Transform) to 1-dimentional signal.

・Scarogram

2. Wavelet Transform

2.1 Wavelet

What is wavelet? It is a function that take time as an argument. Generally, be represented as $Ψ(t)\Psi(t)$ .

Wavelet has a variety of types. Morlet is a famous one.
・Wavelet

To be a proper wavelet, function has to satisfy 2 main constraints.

Zero mean: $∫−∞+∞Ψ(t)dt=0\int^{+ \infty}_{- \infty} \Psi (t)dt = 0$
Finite energy: $∫−∞+∞∣Ψ(t)∣2dt<∞\int^{+ \infty}_{- \infty} |\Psi(t)|^2 dt < \infty$

Example, sin function not satisfy condition 2, because it has infinite energy.
You can think of those constraints as indicating a converging fanction with zero mean.

2.2 Adhust Wavelet

I use wavelet for this transfrom.

Morlet

From here on, I will use Morlet as an example.

・Morlet

When we adjust $t$ , the morlet changes shape accordingly.
・ $(t - b)$
b = 2

・ $(ta)(\dfrac{t}{a})$
a = 0.5

Then, we can make various morlet by use $Ψ(t−ba)\Psi(\dfrac{t-b}{a})$ ,
b = -2, a = 0.25

2.3 Wavelet Transform

Wavelet Transform create image that one axis represent frequency, and one another represent time from signal.

・Wavelet Transform

\rightarrow T(t,f)

The right arrow represent Wavelet Transform, what is going on there? Let's check in below.

・Calculation 1

\cdot \Psi_{a,b}(t)

y

: Be analyzed signal

This calculation creates a ew function that shows how similar each function is. At every point, if the values of each function are close, the new

It represent similarity of two waves y(t) and Ψ.

・Calculation 2

\int^{+\infty}{-\infty} y(t) \cdot \Psi{a,b}(t) dt

This formula calculates the total area of overlap between the Morlet and the original formula, assuming that if the overlapped parts have the same sign, it is a plus, and if they have different signs, it is a minus.
In other words, it shows the degree of overlap between morlet and the original signal.

2.4 Imaginary part

Wavelet has a imaginary part. And as same with real part, T_i(a,b) is calculated by imaginary part.

2.5 Result of Transform

The T(a, b) is the result of WT.
T is a function representing the degree of matching between the Morlet wavelet and the original function.

If you change parameter b, T shows similarity between two functions in the time direction. It makes a one-dimensional wave show where there is a wave that has the same frequency as T (based on parameter a).

When you repeat the same process by changing parameter a, you can obtain various waves that show where each wave has the same frequency corresponding to parameter a.

After obtaining T(a,b) and T_i(a,b) in this way, obtain the square of their absolute value as power.

\sqrt{T(a,b)^2 + T_i(a,b)^2}^2 = T(a,b)^2 + T_i(a,b)^2

When you get this far, all you have to do is color map the frequency components of each power and arrange them in time. And you'll get Scarogram.
It's same operation when you create spectrogram.

It shows at what time in the original waveform a waveform with the same frequency component (varied by a) as the wavelet exists.

summary

This time, I explained about CWT/Scarogram.
thank you for reading.

Reference

(1) Wavelets: a mathematical microscope, Artem Kirsanov

【Analytics】Mel Spectrogram explanation

yuto — Mon, 15 Apr 2024 03:09:37 +0000

Assuming you understand normal spectrograms.

1. Mel Spectrogram

Mel spectrogram is adjusted spectrogram to be easy for humans to understand.
It was made by applying some mel-band filters.

Simply put, it is an enhancement of the low frequency components of the spectrogram.

2. Formula

The process to create Mel Spectrogram contains transform to Mel scale and Hz scale.
Show both formula.

・To Mel scale
$m = 2595 \dot log(1 + \dfrac{f}{700})$
・To Hz scale
$f = 700(10^{m/2575}) - 1$

f: signal(Hz)
m: signal(Mel)

3. Operation

Transform the spectrogram to Mel sacle by above formula.
choose the number of mel-bands.The appropriate number changes depending on the task.(I feel like 128 is generally used.)
Create # bands equally spaced points.
Back apectrogram to Hz scale by above formula.
Create triangle filter based on points(with Mel sacle). This is a Mel-filter.

Finally, we've could obtain Mel Spectrogram.

As you can show, the greater the number of Melbands, the more detailed analysis becomes possible. (Because It makes decrease loss of amplitude with respect to frequency)

Example
・10 Mel-bands

・40 Mel-bands

・80 Mel-bands

::: details code

import librosa
import librosa.display
import matplotlib.pyplot as plt
import numpy as np

# Sample rate (Hz)
sr = 22050
n_fft = 2048  # FFT window size

# Different numbers of Mel bands
mels_list = [40]

plt.figure(figsize=(10, 6))
for n_mels in mels_list:
    mel_filters = librosa.filters.mel(sr=sr, n_fft=n_fft, n_mels=n_mels)
    for i in range(n_mels):
        plt.plot(mel_filters[i], label=f'{n_mels} Mel bands' if i == 0 else "")

plt.title('Effect of Different Numbers of Mel Bands on Mel Filters')
plt.xlabel('Frequency Bin')
plt.ylabel('Filter Amplitude')
plt.legend()
plt.savefig('mel-spectrogram-with-40menbands.png')
plt.show()

:::

:::message
This is a official imprementation of librosa. By decreasing the amplitude of the Mel filter as the frequency increases, the low frequencies become more distinctive.
:::

4. More details

The human sense of sound can capture more in the low frequency range than in the high.
In this time, let's show the mel scale vs Hz scale.
:::details code

import numpy as np
import matplotlib.pyplot as plt

def hz_to_mel(hz):
    """Convert a value in Hertz to Mels."""
    return 2595 * np.log10(1 + hz / 700)

# Generate a range of frequencies from 20 Hz to 20,000 Hz
frequencies = np.linspace(20, 20000, 400)
mel_values = hz_to_mel(frequencies)

# Plot Frequency vs. Mel
plt.figure(figsize=(10, 5))
plt.plot(frequencies, mel_values, label='Hz to Mel', color='blue')
plt.title('Frequency (Hz) to Mel Scale Conversion')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Mel')
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend()
plt.show()

:::

Against the moves in low Hz frequency, Mel frequency moves sensitively. But in high Hz frequency, it don't react that much.

This is the reason of use Mel scale when make filters.

Summary

This time, I explaned about Mel spectrogram.
Thank you for reading!

Reference

(1) Mel Spectrograms Explained Easily