🌾

Kaggle HMS competition Features+Head Starter - [LB 0.35] explained

2024/04/26に公開

0. 概要

https://www.kaggle.com/code/nartaa/features-head-starter-lb-0-35
公開ノートブック複数モデルのアンサンブルを使用した脳波分類タスクのソリューション

0.1 コンペ説明

患者の脳波が示す活動の種類を分類する。専門家の評価を正解として予測を行うコンペ

このコンペでは、生の脳波波形(eeg)とスペクトログラム(spectrogram)の両方が提供されており、スペクトログラムは10分、脳波波形は50秒です、またこれらのデータの中央の50秒は同じデータであり、この期間は同じデータを二つの表現方法で提示しています。

eeg: 一般に想像する波形データ(株価など)
spectrogram: 一般的にスペクトログラムとは、電極から得られた電気信号に対して、短い時間窓で区切ったうえで短時間フーリエ変換(STFT: Short-Time Fourier Transform)を適用して得られる三次元(時間、周波数、パワースペクトル密度)の情報を持ったデータのことである。つまり、datasetのspectrogramファイルはすべてeegデータのフーリエ変換後の情報である。下図の色で囲まれた枠内の電極からのデータで、計4つのスペクトログラム(LL、LP、RP、RR)の作成ができる。

より詳しく知りたい方は、コンペの概要と以下のデータ説明を読むと分かりやすいと思います。
https://www.kaggle.com/competitions/hms-harmful-brain-activity-classification/discussion/468010
https://www.kaggle.com/code/suzukisatsuki/hms-eda

1. アンサンブル(Ensumble)とは

複数のモデルの出力結果を(平均などして)統合する手法です。一般的に異なるモデルを組み合わせるとお互いのエラーを打ち消しあい、過学習を抑制することで性能が向上することが知られています。

2. Notebook解説

2.0 概要

以下の4つのモデルをアンサンブルする。

・Kaggle's spectrograms (CV 0.6365 - LB 0.43)
公式から提供されたスペクトログラムデータを使用したモデル
・Chris's EEG spectrograms(modified version) (CV 0.6336 - LB 0.41)
公式から提供されたeegをChrisさん考案の手法でスペクトログラム化したデータを使用したモデル
・Both Kaggle and EEG spectrograms (CV 0.5726 - LB 0.39)
両方のスペクトログラムを使用したモデル
・Chris's WaveNet (CV 0.6992 - LB 0.41)
WaveNetを使用したモデル

2.1 初期設定

・ライブラリのインポートとモデル設定

import os, random
import tensorflow as tf
import tensorflow
import tensorflow.keras.backend as K
import pandas as pd, numpy as np
import matplotlib.pyplot as plt
from tensorflow.keras.models import load_model

LOAD_BACKBONE_FROM = '/kaggle/input/efficientnetb-tf-keras/EfficientNetB2.h5'
LOAD_MODELS_FROM = '/kaggle/input/features-head-starter-models/'
VER = 37
DATA_TYPE = 'raw' # both|eeg|kaggle|raw
TEST_MODE = False
submission = True

# Setup for ensemble
ENSEMBLE = True
LBs = [0.39,0.41,0.43,0.41] # for weighted ensemble we use LBs of each model
VERK = 33.2 # Kaggle's spectrogram model version
VERB = 35 # Kaggle's and EEG's spectrogram model version
VERE = 34 # EEG's spectrogram model version
VERR = 37 # EEG's raw wavenet model version, trained on single GPU

np.random.seed(42)
random.seed(42)
tf.random.set_seed(42)

# USE SINGLE GPU, MULTIPLE GPUS 
gpus = tf.config.list_physical_devices('GPU')
# WE USE MIXED PRECISION
tf.config.optimizer.set_experimental_options({"auto_mixed_precision": True})
if len(gpus)>1:
    strategy = tf.distribute.MirroredStrategy()
    print(f'Using {len(gpus)} GPUs')
else:
    strategy = tf.distribute.OneDeviceStrategy(device="/gpu:0")
    print(f'Using {len(gpus)} GPU')
  1. LOAD_BACKBONE_FROM = '/kaggle/input/efficientnetb-tf-keras/EfficientNetB2.h5'
    EfficientNetB2モデルの保存先Path
  2. LOAD_MODELS_FROM = '/kaggle/input/features-head-starter-models/'
    モデルの重みの保存先Path
  3. VER = 37
    このNotebookのversion
  4. DATA_TYPE = 'raw' # both|eeg|kaggle|raw
    使用するデータの選択。
    • option
      kaggle: kaggleから提供されたスペクトログラムを使用してモデルを学習
      eeg: kaggleから提供されたeegから生成したスペクトログラムを使用してモデルを学習
      both: 上記の両方
      raw: kaggleから提供されたeegデータをそのままWaveNetで使用
  5. TEST_MODE = False
    Trueで500データをランダムに読み込んで実行される。高速に動作を確認したい場合に有用
  6. submission = True
    提出を行うかの可否
  7. ENSEMBLE = True
    アンサンブルの可否
  8. LBs = [0.39,0.41,0.43,0.41]
    各モデルのLBスコア。アンサンブルする際の重みとして使用
  9. VERK = 33.2 \ VERB = 35 \ VERE = 34 \ VERR = 37
    各モデルのversion
  10. np.random.seed(42) \ random.seed(42) \ tf.random.set_seed(42)
    乱数設定
  11. gpus = tf.config.list_physical_devices('GPU')
    使用可能なGPUの配列を取得
  12. tf.config.optimizer.set_experimental_options({"auto_mixed_precision": True})
    混合浮動小数点を使用(小数点以下の最適化)
  13. if len(gpus)>1: ...
    GPUの単一使用か、複数使用かの設定

2.2 学習用テーブルデータのロード

・eeg_idに重複のない学習用データの作成
推論用データでは、eeg_idの重複が存在しないため、それに合わせて学習用データの重複するeeg_idデータをまとめます。

TARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']
FEATS2 = ['Fp1','T3','C3','O1','Fp2','C4','T4','O2']
FEAT2IDX = {x:y for x,y in zip(FEATS2,range(len(FEATS2)))}

def eeg_from_parquet(parquet_path):

    eeg = pd.read_parquet(parquet_path, columns=FEATS2)
    rows = len(eeg)
    offset = (rows-10_000)//2
    eeg = eeg.iloc[offset:offset+10_000]
    data = np.zeros((10_000,len(FEATS2)))
    for j,col in enumerate(FEATS2):
        
        # FILL NAN
        x = eeg[col].values.astype('float32')
        m = np.nanmean(x)
        if np.isnan(x).mean()<1: x = np.nan_to_num(x,nan=m)
        else: x[:] = 0
        
        data[:,j] = x

    return data

def add_kl(data):
    import torch
    labels = data[TARGETS].values + 1e-5

    # compute kl-loss with uniform distribution by pytorch
    data['kl'] = torch.nn.functional.kl_div(
        torch.log(torch.tensor(labels)),
        torch.tensor([1 / 6] * 6),
        reduction='none'
    ).sum(dim=1).numpy()
    return data

def reset_seed(seed):
    np.random.seed(seed)
    random.seed(seed)
    tf.random.set_seed(seed)
    
if not submission:
    train = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')
    TARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']
    META = ['spectrogram_id','spectrogram_label_offset_seconds','patient_id','expert_consensus']
    train = train.groupby('eeg_id')[META+TARGETS
                           ].agg({**{m:'first' for m in META},**{t:'sum' for t in TARGETS}}).reset_index() 
    train[TARGETS] = train[TARGETS]/train[TARGETS].values.sum(axis=1,keepdims=True)
    train.columns = ['eeg_id','spec_id','offset','patient_id','target'] + TARGETS
    train = add_kl(train)
    print(train.head(1).to_string())
  1. TARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']
    予測すべきターゲット変数のリストを取得
  2. FEATS2 = ['Fp1','T3','C3','O1','Fp2','C4','T4','O2']
    使用する特徴量のリストを取得。この電極から得たデータを使用する
  3. FEAT2IDX = {x:y for x,y in zip(FEATS2,range(len(FEATS2)))}
    文字列の特徴量をインデックスに変換する辞書を作成
  4. def eeg_from_parquet(parquet_path):
    paquetファイルからeegデータを取得する関数
  5. def add_kl(data):
    入力データに"kl"列を追加し、入力データのターゲット列の列方向分布が一様分布(veteが全て1/6の分布)からどの程度離れているかの値を追加する関数。
  6. labels = data[TARGETS].values + 1e-5
    データのターゲット列に微小値を追加したものをlabelsとして定義
  7. data['kl'] = torch.nn.functional.kl_div(torch.log(torch.tensor(labels)), torch.tensor([1 / 6] * 6), reduction='none').sum(dim=1).numpy()
    入力データのターゲット列の各列の確率(スカラ値)が、一様分布(ここでは1/6(スカラ))とどの程度離れているか計算し、それらを行方向に合計することで、一様分布と各行の確率がどの程度離れているか計算する。全く同じ場合'kl'列の値は0になる
    • option
      reduction: 計算した分布を平均や総和でまとめるか否か
  8. return data
    eeg_idなどが記載された訓練用dfに'kl'列を追加して返す
  9. def reset_seed(seed):
    seed値をリセット
  10. if not submission:
    データ確認用スコープ。提出時には実行時間を早めるため通らない
  11. train = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')
    学習用データの取得
  12. TARGETS = ['seizure_vote', 'lpd_vote', 'gpd_vote', 'lrda_vote', 'grda_vote', 'other_vote']
    ターゲット列名のリスト作成
  13. META = ['spectrogram_id','spectrogram_label_offset_seconds','patient_id','expert_consensus']
    各データポイントに対するメタデータの列名リスト作成
  14. train = train.groupby('eeg_id')[META+TARGETS].agg({{m:'first' for m in META},{t:'sum' for t in TARGETS}}).reset_index()
    学習用データを一意のeeg_id毎にまとめる。メタデータは最初の値を、TARGETデータは合計値を取得する
  15. train[TARGETS] = train[TARGETS]/train[TARGETS].values.sum(axis=1,keepdims=True)
    TARGET列の各行をその行の合計値で割って値を確率に変更
  16. train.columns = ['eeg_id','spec_id','offset','patient_id','target'] + TARGETS
    列名をわかりやすいように変更
  17. train = add_kl(train)
    学習用データに'kl'列を追加

2.3 学習用スペクトログラム/eegのロード

・Chrisさん作成のデータセットから学習用データをロード

%%time
if not submission:
    # FOR TESTING SET READ_FILES TO TRUE
    if TEST_MODE:
        train = train.sample(500,random_state=42).reset_index(drop=True)
        spectrograms = {}
        for i,e in enumerate(train.spec_id.values):
            if i%100==0: print(i,', ',end='')
            x = pd.read_parquet(f'/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/{e}.parquet')
            spectrograms[e] = x.values
        all_eegs = {}
        for i,e in enumerate(train.eeg_id.values):
            if i%100==0: print(i,', ',end='')
            x = np.load(f'/kaggle/input/eeg-spectrograms/EEG_Spectrograms/{e}.npy')
            all_eegs[e] = x
        all_raw_eegs = {}
        for i,e in enumerate(train.eeg_id.values):
            if i%100==0: print(i,', ',end='')
            x = eeg_from_parquet(f'/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/{e}.parquet')              
            all_raw_eegs[e] = x
    else:
        spectrograms = None
        all_eegs = None
        all_raw_eegs = None
        if DATA_TYPE=='both' or DATA_TYPE=='kaggle':
            spectrograms = np.load('/kaggle/input/brain-spectrograms/specs.npy',allow_pickle=True).item()
        if DATA_TYPE=='both' or DATA_TYPE=='eeg':
            all_eegs = np.load('/kaggle/input/eeg-spectrograms/eeg_specs.npy',allow_pickle=True).item()
        if DATA_TYPE=='raw':
            all_raw_eegs = np.load('/kaggle/input/brain-eegs/eegs.npy',allow_pickle=True).item()

使用するスペクトログラム/eegをロードします。

  1. spectrograms = np.load('/kaggle/input/brain-spectrograms/specs.npy',allow_pickle=True).item()
    np.loadでnpyファイルをロード。リストや辞書など特定のnumpy配列以外のデータはpickle化されて書き込まれるため、allow_pickleで取得(pickleは実行可能コードを含むので、信頼性の観点からデフォルトはFalse)。.item()でnumpy配列から単一要素を抜き出す(.item()は[0]と似た使い方で、一つの要素を含むnumpy配列にのみ使用されるメソッド)。

2.4 データジェネレータ定義

・データジェネレータを定義
出力は512×512×3です。

import albumentations as albu
from scipy.signal import butter, lfilter

class DataGenerator():
    'Generates data for Keras'
    def __init__(self, data, specs=None, eeg_specs=None, raw_eegs=None, augment=False, mode='train', data_type=DATA_TYPE): 
        self.data = data
        self.augment = augment
        self.mode = mode
        self.data_type = data_type
        self.specs = specs
        self.eeg_specs = eeg_specs
        self.raw_eegs = raw_eegs
        self.on_epoch_end()
        
    def __len__(self):
        return self.data.shape[0]

    def __getitem__(self, index):
        X, y = self.data_generation(index)
        if self.augment: X = self.augmentation(X)
        return X, y
    
    def __call__(self):
        for i in range(self.__len__()):
            yield self.__getitem__(i)
            
            if i == self.__len__()-1:
                self.on_epoch_end()
                
    def on_epoch_end(self):
        if self.mode=='train': 
            self.data = self.data.sample(frac=1).reset_index(drop=True)
    
    def data_generation(self, index):
        if self.data_type == 'both':
            X,y = self.generate_all_specs(index)
        elif self.data_type == 'eeg' or self.data_type == 'kaggle':
            X,y = self.generate_specs(index)
        elif self.data_type == 'raw':
            X,y = self.generate_raw(index)

        return X,y
    
    def generate_all_specs(self, index):
        X = np.zeros((512,512,3),dtype='float32')
        y = np.zeros((6,),dtype='float32')
        
        row = self.data.iloc[index]
        if self.mode=='test': 
            offset = 0
        else:
            offset = int(row.offset/2)
            
        eeg = self.eeg_specs[row.eeg_id]
        spec = self.specs[row.spec_id]
        
        imgs = [spec[offset:offset+300,k*100:(k+1)*100].T for k in [0,2,1,3]] # to match kaggle with eeg
        img = np.stack(imgs,axis=-1)
        # LOG TRANSFORM SPECTROGRAM
        img = np.clip(img,np.exp(-4),np.exp(8))
        img = np.log(img)
            
        # STANDARDIZE PER IMAGE
        img = np.nan_to_num(img, nan=0.0)    
            
        mn = img.flatten().min()
        mx = img.flatten().max()
        ep = 1e-5
        img = 255 * (img - mn) / (mx - mn + ep)
        
        X[0_0+56:100+56,:256,0] = img[:,22:-22,0] # LL_k
        X[100+56:200+56,:256,0] = img[:,22:-22,2] # RL_k
        X[0_0+56:100+56,:256,1] = img[:,22:-22,1] # LP_k
        X[100+56:200+56,:256,1] = img[:,22:-22,3] # RP_k
        X[0_0+56:100+56,:256,2] = img[:,22:-22,2] # RL_k
        X[100+56:200+56,:256,2] = img[:,22:-22,1] # LP_k
        
        X[0_0+56:100+56,256:,0] = img[:,22:-22,0] # LL_k
        X[100+56:200+56,256:,0] = img[:,22:-22,2] # RL_k
        X[0_0+56:100+56,256:,1] = img[:,22:-22,1] # LP_k
        X[100+56:200+56,256:,1] = img[:,22:-22,3] # RP_K
        
        # EEG
        img = eeg
        mn = img.flatten().min()
        mx = img.flatten().max()
        ep = 1e-5
        img = 255 * (img - mn) / (mx - mn + ep)
        X[200+56:300+56,:256,0] = img[:,22:-22,0] # LL_e
        X[300+56:400+56,:256,0] = img[:,22:-22,2] # RL_e
        X[200+56:300+56,:256,1] = img[:,22:-22,1] # LP_e
        X[300+56:400+56,:256,1] = img[:,22:-22,3] # RP_e
        X[200+56:300+56,:256,2] = img[:,22:-22,2] # RL_e
        X[300+56:400+56,:256,2] = img[:,22:-22,1] # LP_e
        
        X[200+56:300+56,256:,0] = img[:,22:-22,0] # LL_e
        X[300+56:400+56,256:,0] = img[:,22:-22,2] # RL_e
        X[200+56:300+56,256:,1] = img[:,22:-22,1] # LP_e
        X[300+56:400+56,256:,1] = img[:,22:-22,3] # RP_e

        if self.mode!='test':
            y[:] = row[TARGETS]
        
        return X,y
    
    def generate_specs(self, index):
        X = np.zeros((512,512,3),dtype='float32')
        y = np.zeros((6,),dtype='float32')
        
        row = self.data.iloc[index]
        if self.mode=='test': 
            offset = 0
        else:
            offset = int(row.offset/2)
            
        if self.data_type == 'eeg':
            img = self.eeg_specs[row.eeg_id]
        elif self.data_type == 'kaggle':
            spec = self.specs[row.spec_id]
            imgs = [spec[offset:offset+300,k*100:(k+1)*100].T for k in [0,2,1,3]] # to match kaggle with eeg
            img = np.stack(imgs,axis=-1)
            # LOG TRANSFORM SPECTROGRAM
            img = np.clip(img,np.exp(-4),np.exp(8))
            img = np.log(img)
            
            # STANDARDIZE PER IMAGE
            img = np.nan_to_num(img, nan=0.0)    
            
        mn = img.flatten().min()
        mx = img.flatten().max()
        ep = 1e-5
        img = 255 * (img - mn) / (mx - mn + ep)
        
        X[0_0+56:100+56,:256,0] = img[:,22:-22,0]
        X[100+56:200+56,:256,0] = img[:,22:-22,2]
        X[0_0+56:100+56,:256,1] = img[:,22:-22,1]
        X[100+56:200+56,:256,1] = img[:,22:-22,3]
        X[0_0+56:100+56,:256,2] = img[:,22:-22,2]
        X[100+56:200+56,:256,2] = img[:,22:-22,1]
        
        X[0_0+56:100+56,256:,0] = img[:,22:-22,0]
        X[100+56:200+56,256:,0] = img[:,22:-22,1]
        X[0_0+56:100+56,256:,1] = img[:,22:-22,2]
        X[100+56:200+56,256:,1] = img[:,22:-22,3]
        
        X[200+56:300+56,:256,0] = img[:,22:-22,0]
        X[300+56:400+56,:256,0] = img[:,22:-22,1]
        X[200+56:300+56,:256,1] = img[:,22:-22,2]
        X[300+56:400+56,:256,1] = img[:,22:-22,3]
        X[200+56:300+56,:256,2] = img[:,22:-22,3]
        X[300+56:400+56,:256,2] = img[:,22:-22,2]
        
        X[200+56:300+56,256:,0] = img[:,22:-22,0]
        X[300+56:400+56,256:,0] = img[:,22:-22,2]
        X[200+56:300+56,256:,1] = img[:,22:-22,1]
        X[300+56:400+56,256:,1] = img[:,22:-22,3]
        
        if self.mode!='test':
            y[:] = row[TARGETS]
        
        return X,y
    
    def generate_raw(self,index):
        X = np.zeros((10_000,8),dtype='float32')
        y = np.zeros((6,),dtype='float32')
        
        row = self.data.iloc[index]
        eeg = self.raw_eegs[row.eeg_id]
            
        # FEATURE ENGINEER
        X[:,0] = eeg[:,FEAT2IDX['Fp1']] - eeg[:,FEAT2IDX['T3']]
        X[:,1] = eeg[:,FEAT2IDX['T3']] - eeg[:,FEAT2IDX['O1']]
            
        X[:,2] = eeg[:,FEAT2IDX['Fp1']] - eeg[:,FEAT2IDX['C3']]
        X[:,3] = eeg[:,FEAT2IDX['C3']] - eeg[:,FEAT2IDX['O1']]
            
        X[:,4] = eeg[:,FEAT2IDX['Fp2']] - eeg[:,FEAT2IDX['C4']]
        X[:,5] = eeg[:,FEAT2IDX['C4']] - eeg[:,FEAT2IDX['O2']]
            
        X[:,6] = eeg[:,FEAT2IDX['Fp2']] - eeg[:,FEAT2IDX['T4']]
        X[:,7] = eeg[:,FEAT2IDX['T4']] - eeg[:,FEAT2IDX['O2']]
            
        # STANDARDIZE
        X = np.clip(X,-1024,1024)
        X = np.nan_to_num(X, nan=0) / 32.0
            
        # BUTTER LOW-PASS FILTER
        X = self.butter_lowpass_filter(X)
        # Downsample
        X = X[::5,:]
        
        if self.mode!='test':
            y[:] = row[TARGETS]
                
        return X,y
        
    def butter_lowpass_filter(self, data, cutoff_freq=20, sampling_rate=200, order=4):
        nyquist = 0.5 * sampling_rate
        normal_cutoff = cutoff_freq / nyquist
        b, a = butter(order, normal_cutoff, btype='low', analog=False)
        filtered_data = lfilter(b, a, data, axis=0)
        return filtered_data
    
    def resize(self, img,size):
        composition = albu.Compose([
                albu.Resize(size[0],size[1])
            ])
        return composition(image=img)['image']
            
    def augmentation(self, img):
        composition = albu.Compose([
                albu.HorizontalFlip(p=0.4)
            ])
        return composition(image=img)['image']
  1. def __init__(self, data, specs=None, eeg_specs=None, raw_eegs=None, augment=False, mode='train', data_type=DATA_TYPE):
    クラスの初期化。
    • 引数
      data: 入力データ
      specs: idをkey、スペクトログラムデータをvalueとする辞書
      eeg_specs: idをkey、eegスペクトログラムデータをvalueとする辞書
      raw_eegs: idをkey、eegデータをvalueとする辞書
      augment: データ拡張の有無
      mode: ジェネレータの動作モードを指定
      data_type: どのデータを使用するか指定
  2. def __len__(self): \ return self.data.shape[0]
    データの長さ(データポイント数)を返す
  3. __getitem__(self, index): \ X, y = self.data_generation(index)
    指定indexの学習用データと正解ラベルをタプルで返す。バッチサイズはfit()により自動的に処理される
  4. def __call__(self):
    この関数によりDataGeneratorインスタンスはPythonのジェネレータとして機能する
  5. for i in range(self.__len__()): \ yield self.__getitem__(i)
    このインスタンスが呼ばれたときに、lenの数だけ繰り返し(X, y)を返す
  6. if i == self.__len__()-1: \ self.on_epoch_end()
    エポック終了時、on_epoch_end関数を呼ぶ
  7. def on_epoch_end(self):
    エポック終了時に呼ばれる関数
  8. if self.mode=='train': \ self.data = self.data.sample(frac=1).reset_index(drop=True)
    モデル学習時のみ、エポックの終了タイミングでデータのシャッフルを行う。data.sample()でデータをランダムにサンプリングして新しいデータセットを作成し、frac=1で100%のデータを選択するように指定。reset_index(drop=True)により、古いデータのインデックスは保持しないようにして新しいインデックスが割り振られる。
  9. def data_generation(self, index): ... return X,y
    使用するデータに合わせてジェネレーターを選択し、データとラベルを返す
  10. def generate_all_specs(self, index):
    kaggleから提供されたスペクトログラムデータと、eegから生成したスペクトログラムデータの両方を取得する関数
  11. X = np.zeros((512,512,3),dtype='float32')
    入力データを格納する変数の初期化。kaggleスペクトログラムデータとeegスペクトログラムデータの両方が含まれる。
  12. y = np.zeros((6,),dtype='float32')
    対応する正解ラベルを格納する変数の初期化
  13. row = self.data.iloc[index]
    指定されたインデックスの行を入力から取り出す
  14. if self.mode=='test': \ offset = 0
    テストデータにoffsetは存在しない。(複数のeeg_idが存在しないため、offsetという概念がない)
  15. else: \ offset = int(row.offset/2)
    学習データはoffsetの半分前からデータを取得
  16. eeg = self.eeg_specs[row.eeg_id] \ spec = self.specs[row.spec_id]
    それぞれのスペクトラムデータの辞書から指定のeeg_idのスペクトログラムを取得
  17. imgs = [spec[offset:offset+300,k*100:(k+1)*100].T for k in range(4)]
    specは縦軸が時間、横軸が周波数(Hz)、値がその周波数帯の信号の強さを表している。従って、ここでは時間軸300と周波数軸100要素を抜き出し、転置して、時間軸に対する周波数の強さを表すスペクトログラムの二次元配列を取得している。
  18. img = np.stack(imgs,axis=-1)
    スペクトログラムの二次元配列について、周波数要素を100ずつ0~400まで重ねて取得している
  19. img = np.clip(img,np.exp(-4),np.exp(8))
    スペクトログラム画像(二次元配列)の最小値をe^-4≒0.018、最大値をe^8≒2981に制限
  20. img = np.log(img)
    スペクトログラム画像を対数変換して、処理や視覚化を用意にする。またデータレンジを圧縮する
  21. img = np.nan_to_num(img, nan=0.0)
    スペクトログラム画像のNaNを0に置換
  22. mn = img.flatten().min() \ mx = img.flatten().max()
    画像を1次元に平滑化して、最大値と最小値を取得
  23. img = 255 * (img - mn) / (mx - mn + ep)
    スペクトログラム画像を0~255の範囲に正規化
  24. X[0_0+56:100+56,:256,0] = img[:,22:-22,0] ...
    得られた形状(100, 300)の4枚のスペクトログラム画像を、それぞれ中心の(100, 256)を切り抜いて、一つの画像を8つに分割した各部に左上から割り当てていく。最終的な入力は(512, 512, 3)にしたいので、kaggle_specとeeg_specを4枚ずつ割り当て。
    この順番はコンペ参加者との議論によって決定されたが、著者曰くあまり改善は見られなかったとのこと。
  25. img = eeg
    スペクトログラム画像をeegから生成したものに変更して同じ操作を行う。
  26. if self.mode!='test': \ y[:] = row[TARGETS]
    評価モードでない場合、正解ラベルを取得
  27. return X,y
    8つのスペクトログラム画像(numpy二次元配列)を一つの画像にまとめたデータをX,正解ラベル分布をyとして返す
  28. def generate_specs(self, index):
    kaggle_specもしくはeeg_specのみを使用する場合のジェネレータ関数
  29. X = np.zeros((512,512,3),dtype='float32') \ y = np.zeros((6,),dtype='float32')
    空の入力と出力を定義
  30. row = self.data.iloc[index]
    指定indexのtrainデータを取得
    以下、def generate_all_specs(self, index):と同じ
  31. def generate_raw(self,index):
    名前の1次元eegデータを使用する場合のジェネレータ関数。詳しくはChrisさんのWaveNet解法を確認して下さい。
    入力形状を定義してデータをロードして、電極間の差異のデータを取得して8つの系列データを作成する。各系列データは正規化、ローパスフィルター、ダウンサンプリングを経て波形データとなる。
  32. def butter_lowpass_filter(self, data, cutoff_freq=20, sampling_rate=200, order=4):
    ローパスフィルター。同様にWaveNet解法を確認して下さい。カットオフ周波数を計算してbutterローパスフィルターを使用。
  33. def resize(self, img,size):
    データ拡張に使用するリサイズ関数。size[0]からsize[1]にリサイズ。
  34. def augmentation(self, img)
    データ拡張に使用する反転関数。40%の確率で左右反転。

・データジェネレータの生成物を可視化

if not submission and DATA_TYPE!='raw':
    gen = DataGenerator(train, augment=False, specs=spectrograms, eeg_specs=all_eegs, data_type=DATA_TYPE)
    for x,y in gen:
        break
    plt.imshow(x[:,:,0])
    plt.title(f'Target = {y.round(1)}',size=12)
    plt.yticks([])
    plt.ylabel('Frequencies (Hz)',size=12)
    plt.xlabel('Time (sec)',size=12)
    plt.savefig('dataX_'+DATA_TYPE+'.png')
    plt.show()
    
if not submission and DATA_TYPE=='raw':
    gen = DataGenerator(train, raw_eegs=all_raw_eegs, data_type=DATA_TYPE)
    for x,y in gen:
        plt.figure(figsize=(20,4))
        offset = 0
        for j in range(x.shape[-1]):
            if j!=0: offset -= x[:,j].min()
            plt.plot(range(2_000),x[:,j]+offset,label=f'feature {j+1}')
            offset += x[:,j].max()
        plt.legend()
        plt.savefig('dataX_'+DATA_TYPE+'.png')
        plt.show()
        break

・output(X(画像),y(タイトル))
DATA_TYPE = 'both'

DATA_TYPE = 'eeg' (or 'kaggle')

DATA_TYPE = 'raw'

2.5 モデルの構築

・学習率の設定

if not submission:

    def lrfn(epoch):
        e3 = 1e-3 if DATA_TYPE=='raw' else 1e-4
        return [1e-3,1e-3,e3,1e-4,1e-5][epoch]

    LR = tf.keras.callbacks.LearningRateScheduler(lrfn, verbose = True)
    
    def lrfn2(epoch):
        return [1e-5,1e-5,1e-6][epoch]

    LR2 = tf.keras.callbacks.LearningRateScheduler(lrfn2, verbose = True)

学習率スケジューラーを定義

・モデル構築、関数定義

from tensorflow.keras.layers import Input, Dense, Multiply, Add, Conv1D, Concatenate

def build_model():  
    inp = tf.keras.layers.Input((512,512,3))
    base_model = load_model(f'{LOAD_BACKBONE_FROM}')    
    x = base_model(inp)
    x = tf.keras.layers.GlobalAveragePooling2D()(x)
    output = tf.keras.layers.Dense(6,activation='softmax', dtype='float32')(x)
    model = tf.keras.Model(inputs=inp, outputs=output)
    opt = tf.keras.optimizers.Adam(learning_rate = 1e-3)
    loss = tf.keras.losses.KLDivergence()
    model.compile(loss=loss, optimizer=opt)  
    return model

def score(y_true, y_pred):
    kl = tf.keras.metrics.KLDivergence()
    return kl(y_true, y_pred)

def wave_block(x, filters, kernel_size, n):
    dilation_rates = [2**i for i in range(n)]
    x = Conv1D(filters = filters,
               kernel_size = 1,
               padding = 'same')(x)
    res_x = x
    for dilation_rate in dilation_rates:
        tanh_out = Conv1D(filters = filters,
                          kernel_size = kernel_size,
                          padding = 'same', 
                          activation = 'tanh', 
                          dilation_rate = dilation_rate)(x)
        sigm_out = Conv1D(filters = filters,
                          kernel_size = kernel_size,
                          padding = 'same',
                          activation = 'sigmoid', 
                          dilation_rate = dilation_rate)(x)
        x = Multiply()([tanh_out, sigm_out])
        x = Conv1D(filters = filters,
                   kernel_size = 1,
                   padding = 'same')(x)
        res_x = Add()([res_x, x])
    return res_x

def build_wave_model():
        
    # INPUT 
    inp = tf.keras.Input(shape=(2_000,8))
    
    ############
    # FEATURE EXTRACTION SUB MODEL
    inp2 = tf.keras.Input(shape=(2_000,1))
    x = wave_block(inp2, 8, 4, 6)
    x = wave_block(x, 16, 4, 6)
    x = wave_block(x, 32, 4, 6)
    x = wave_block(x, 64, 4, 6)
    model2 = tf.keras.Model(inputs=inp2, outputs=x)
    ###########
    
    # LEFT TEMPORAL CHAIN
    x1 = model2(inp[:,:,0:1])
    x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)
    x2 = model2(inp[:,:,1:2])
    x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)
    z1 = tf.keras.layers.Average()([x1,x2])
    
    # LEFT PARASAGITTAL CHAIN
    x1 = model2(inp[:,:,2:3])
    x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)
    x2 = model2(inp[:,:,3:4])
    x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)
    z2 = tf.keras.layers.Average()([x1,x2])
    
    # RIGHT PARASAGITTAL CHAIN
    x1 = model2(inp[:,:,4:5])
    x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)
    x2 = model2(inp[:,:,5:6])
    x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)
    z3 = tf.keras.layers.Average()([x1,x2])
    
    # RIGHT TEMPORAL CHAIN
    x1 = model2(inp[:,:,6:7])
    x1 = tf.keras.layers.GlobalAveragePooling1D()(x1)
    x2 = model2(inp[:,:,7:8])
    x2 = tf.keras.layers.GlobalAveragePooling1D()(x2)
    z4 = tf.keras.layers.Average()([x1,x2])
    
    # COMBINE CHAINS
    y = tf.keras.layers.Concatenate()([z1,z2,z3,z4])
    y = tf.keras.layers.Dense(64, activation='relu')(y)
    y = tf.keras.layers.Dense(6,activation='softmax', dtype='float32')(y)
    
    # COMPILE MODEL
    model = tf.keras.Model(inputs=inp, outputs=y)
    opt = tf.keras.optimizers.Adam(learning_rate = 1e-3)
    loss = tf.keras.losses.KLDivergence()
    model.compile(loss=loss, optimizer = opt)
    
    return model

def plot_hist(hist):
    metrics = ['loss']
    for i,metric in enumerate(metrics):
        plt.figure(figsize=(10,4))
        plt.subplot(1,2,i+1)
        plt.plot(hist[metric])
        plt.plot(hist[f'val_{metric}'])
        plt.title(f'{metric}',size=12)
        plt.ylabel(f'{metric}',size=12)
        plt.xlabel('epoch',size=12)
        plt.legend(["train", "validation"], loc="upper left")
        plt.show()
  1. def build_model():
    ベースモデルを定義する関数
  2. inp = tf.keras.layers.Input((512,512,3))
    入力形状の定義
  3. base_model = load_model(f'{LOAD_BACKBONE_FROM}')
    ベースモデルの重みを取得
  4. x = base_model(inp)
    入力をベースモデルに通す
  5. x = tf.keras.layers.GlobalAveragePooling2D()(x)
    ベースモデル出力をグローバルアベレージプーリングに通す。高さ及び幅方向の全てのデータに対して平均をとり、(1,1,チャネル数)形状のデータを出力します。これにより、パラメータの削減による過学習リスクの削減と、検出物体の位置によらない画像分類が可能となる。
  6. output = tf.keras.layers.Dense(6,activation='softmax', dtype='float32')(x)
    出力数6、活性化関数がsoftmaxの全結合層で、出力を確率分布として扱えるようにする。
  7. model = tf.keras.Model(inputs=inp, outputs=output)
    入力と出力を指定し、モデルを定義
  8. opt = tf.keras.optimizers.Adam(learning_rate = 1e-3)
    最適化関数
  9. loss = tf.keras.losses.KLDivergence()
    損失関数
  10. model.compile(loss=loss, optimizer=opt)
    最適化関数と損失関数を指定してコンパイル
  11. return model
    モデルを返す
  12. def score(y_true, y_pred):
    二つの確率分布を受け取り、klDivスコアを返す関数
  13. def wave_block(x, filters, kernel_size, n):
    WaveNetを構成するブロックを定義。詳細はWaveNet解法を参照。
  • 引数
    x: 入力
    filters: 畳み込みフィルター(カーネル)数
    kernel_size: カーネルサイズ
    n: 拡張因果畳み込みの拡張数=畳み込み回数
  1. def build_wave_model():
    WaveNetを用いたモデルを構築。電極感のデータを使用して4つの系列を作成して畳み込みを行う。
    詳細は同上
  2. def plot_hist(hist):
    損失関数の履歴配列を受け取ってグラフ化する関数

2.6 転移学習

・学習及び交差検証

from sklearn.model_selection import KFold, GroupKFold
import tensorflow.keras.backend as K, gc

if not submission:
    # for CV scores setting random seed works for single GPU only
    reset_seed(42)
    all_oof = []
    all_true = []
    losses = []
    val_losses = []
    total_hist = {}

    gkf = GroupKFold(n_splits=5)
    for i, (train_index, valid_index) in enumerate(gkf.split(train, train.target, train.patient_id)):   
        
        print('#'*25)
        print(f'### Fold {i+1}')
        
        data, val = train.iloc[train_index],train.iloc[valid_index]
        train_gen = DataGenerator(data, augment=False, specs=spectrograms, eeg_specs=all_eegs, raw_eegs=all_raw_eegs)
        valid_gen = DataGenerator(val, mode='valid', specs=spectrograms, eeg_specs=all_eegs, raw_eegs=all_raw_eegs)
        data, val = data[data['kl']<5.5],val[val['kl']<5.5]
        train_gen2 = DataGenerator(data, augment=False, specs=spectrograms, eeg_specs=all_eegs, raw_eegs=all_raw_eegs)
        in_shape = (2000,8) if DATA_TYPE=='raw' else (512,512,3)
        EPOCHS = 5 if DATA_TYPE=='raw' else 4
        BATCH_SIZE_PER_REPLICA = 8 if DATA_TYPE=='raw' else 32
        BATCH_SIZE = BATCH_SIZE_PER_REPLICA * strategy.num_replicas_in_sync

        train_dataset = tf.data.Dataset.from_generator(generator=train_gen, 
                                                   output_signature=(tf.TensorSpec(shape=in_shape, dtype=tf.float32),
                                                                     tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(BATCH_SIZE).prefetch(tf.data.AUTOTUNE)
        val_dataset = tf.data.Dataset.from_generator(generator=valid_gen, 
                                                   output_signature=(tf.TensorSpec(shape=in_shape, dtype=tf.float32),
                                                                     tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(BATCH_SIZE).prefetch(tf.data.AUTOTUNE)
        train_dataset2 = tf.data.Dataset.from_generator(generator=train_gen2, 
                                                   output_signature=(tf.TensorSpec(shape=in_shape, dtype=tf.float32),
                                                                     tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(BATCH_SIZE).prefetch(tf.data.AUTOTUNE)
            
        print(f'### train size {len(train_index)}, valid size {len(valid_index)}')
        print('#'*25)
        
        K.clear_session()
        with strategy.scope():
            if DATA_TYPE=='raw':
                model = build_wave_model()
            else:
                model = build_model()
        
        hist = model.fit(train_dataset, validation_data = val_dataset, 
                         epochs=EPOCHS, callbacks=[LR])
        print(f'### seconds stage train size {len(data)}, valid size {len(val)}')
        print('#'*25)
        hist2 = model.fit(train_dataset2, validation_data = val_dataset, 
                         epochs=3, callbacks=[LR2])
        losses.append(hist.history['loss']+hist2.history['loss'])
        val_losses.append(hist.history['val_loss']+hist2.history['val_loss'])
        with strategy.scope():
            model.save_weights(f'model_{DATA_TYPE}_{VER}_{i}.weights.h5')
        oof = model.predict(val_dataset, verbose=1)
        all_oof.append(oof)
        all_true.append(train.iloc[valid_index][TARGETS].values)    
        del model, oof
        gc.collect()
        
    total_hist['loss'] = np.mean(losses,axis=0)
    total_hist['val_loss'] = np.mean(val_losses,axis=0)
    all_oof = np.concatenate(all_oof)
    all_true = np.concatenate(all_true)
    plot_hist(total_hist)
    print('#'*25)
    print(f'CV KL SCORE: {score(all_true,all_oof)}')
  1. if not submission:
    提出する場合は実行速度のため、交差検証を行わない
  2. gkf = GroupKFold(n_splits=5)
    GroupKFoldオブジェクトを定義。分割数は5
  3. for i, (train_index, valid_index) in enumerate(gkf.split(train, train.target, train.patient_id)):
    GroupKFoldオブジェクトをイテレートして学習と評価に使用するデータのインデックス配列を取得。繰り返し回数は分割数と同じ
  4. data, val = train.iloc[train_index],train.iloc[valid_index]
    trainデータから学習用及び評価用インデックスのデータを取得
  5. train_gen = DataGenerator(data, augment=False, specs=spectrograms, eeg_specs=all_eegs, raw_eegs=all_raw_eegs)
    学習用データジェネレータの定義
  6. valid_gen = DataGenerator(val, mode='valid', specs=spectrograms, eeg_specs=all_eegs, raw_eegs=all_raw_eegs)
    評価用データジェネレータの定義
  7. data, val = data[data['kl']<5.5],val[val['kl']<5.5]
    学習用、評価用データについて、kl距離が5.5以下のもののみ保持する
  8. train_gen2 = DataGenerator(data, augment=False, specs=spectrograms, eeg_specs=all_eegs, raw_eegs=all_raw_eegs)
    kl距離が5.5以下のデータのみを使用するジェネレータを定義
  9. in_shape = (2000,8) if DATA_TYPE=='raw' else (512,512,3)
    入力形状を取得。rawデータとspectramデータで区別
  10. EPOCHS = 5 if DATA_TYPE=='raw' else 4
    エポック数の指定
  11. BATCH_SIZE_PER_REPLICA = 8 if DATA_TYPE=='raw' else 32
    並列で動く各GPUでの学習に使用されるバッチサイズ
  12. BATCH_SIZE = BATCH_SIZE_PER_REPLICA * strategy.num_replicas_in_sync
    全てのGPUで処理されるバッチサイズの総量を計算
  13. train_dataset = tf.data.Dataset.from_generator(generator=train_gen, output_signature=(tf.TensorSpec(shape=in_shape, dtype=tf.float32), tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(BATCH_SIZE).prefetch(tf.data.AUTOTUNE)
    from_generator関数による学習用データセットの定義。
    ・引数
    generator: 定義したジェネレータ
    output_signature: 入力形状と出力形状を渡す
    .batch: バッチサイズを渡す
    .prefetch: GPUの計算リソース活用のために、学習中に次のデータを事前に読み込む際の挙動を指定。AUTOTUNEにするとバッファサイズを調整して、最適なパフォーマンスを自動的に決定する
  14. val_dataset, train_dataset2
    同上。評価用およびklフィルタ付きデータのデータセットを定義
  15. K.clear_session()
    kerasが生成したデータの初期化。メモリ解放
  16. with strategy.scope():
    GPUの設定を使用するスコープ。このスコープ内でモデルを定義。
  17. model = build_wave_model() \ model = build_model()
    モデルの定義。build_model()は内部で重みをロードしているため転移学習。
  18. hist = model.fit(train_dataset, validation_data = val_dataset, epochs=EPOCHS, callbacks=[LR])
    データセットを使用してモデルを学習させる。エポックごとにcallbacksが呼ばれて学習率が変更される
  19. hist2 = model.fit(train_dataset2, validation_data = val_dataset, epochs=3, callbacks=[LR2])
    klフィルタを通したデータで学習。エポック数は3
  20. losses.append(hist.history['loss']+hist2.history['loss'])
    損失関数の値を取得。2つのデータセットから取得した値を足して配列に追加
  21. val_losses.append(hist.history['val_loss']+hist2.history['val_loss'])
    評価用データに対する損失関数の値を取得
  22. model.save_weights(f'model_{DATA_TYPE}{VER}{i}.weights.h5')
    モデルの重みを保存
  23. oof = model.predict(val_dataset, verbose=1)
    評価用データセットを用いて推論、結果をoofに代入
  24. all_oof.append(oof)
    fold毎に、評価用データを用いた推論の結果をall_oof配列に追加
  25. all_true.append(train.iloc[valid_index][TARGETS].values)
    今回のfoldで使用した評価用データの正解分布をall_ture配列に追加
  26. del model, oof \ gc.collect()
    メモリ、ガベージコレクションの解放
  27. total_hist['loss'] = np.mean(losses,axis=0)
    fold全体の損失の平均を計算してtotal_histの'loss'列に代入
  28. total_hist['val_loss'] = np.mean(val_losses,axis=0)
    fold全体の評価用データ使用時の損失の平均を計算してtotal_histの'val_loss'列に代入
  29. all_oof = np.concatenate(all_oof)
    fold毎の評価用データを用いた推論の結果が格納されているall_oofを列方向に結合。二次元配列となる
  30. all_true = np.concatenate(all_true)
    同上。正解分布を列方向に結合
  31. plot_hist(total_hist)
    fold毎の損失を可視化

2.7 eegからスペクトログラム生成

・eegの前処理とスペクトログラムの作成

import pywt, librosa

USE_WAVELET = None 

NAMES = ['LL','LP','RP','RR']

FEATS = [['Fp1','F7','T3','T5','O1'],
         ['Fp1','F3','C3','P3','O1'],
         ['Fp2','F8','T4','T6','O2'],
         ['Fp2','F4','C4','P4','O2']]

# DENOISE FUNCTION
def maddest(d, axis=None):
    return np.mean(np.absolute(d - np.mean(d, axis)), axis)

def denoise(x, wavelet='haar', level=1):    
    coeff = pywt.wavedec(x, wavelet, mode="per")
    sigma = (1/0.6745) * maddest(coeff[-level])

    uthresh = sigma * np.sqrt(2*np.log(len(x)))
    coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])

    ret=pywt.waverec(coeff, wavelet, mode='per')
    
    return ret

import librosa

def spectrogram_from_eeg(parquet_path, display=False):
    
    # LOAD MIDDLE 50 SECONDS OF EEG SERIES
    eeg = pd.read_parquet(parquet_path)
    middle = (len(eeg)-10_000)//2
    eeg = eeg.iloc[middle:middle+10_000]
    
    # VARIABLE TO HOLD SPECTROGRAM
    img = np.zeros((100,300,4),dtype='float32')
    
    if display: plt.figure(figsize=(10,7))
    signals = []
    for k in range(4):
        COLS = FEATS[k]
        
        for kk in range(4):
            # FILL NANS
            x1 = eeg[COLS[kk]].values
            x2 = eeg[COLS[kk+1]].values
            m = np.nanmean(x1)
            if np.isnan(x1).mean()<1: x1 = np.nan_to_num(x1,nan=m)
            else: x1[:] = 0
            m = np.nanmean(x2)
            if np.isnan(x2).mean()<1: x2 = np.nan_to_num(x2,nan=m)
            else: x2[:] = 0
                
            # COMPUTE PAIR DIFFERENCES
            x = x1 - x2

            # DENOISE
            if USE_WAVELET:
                x = denoise(x, wavelet=USE_WAVELET)
            signals.append(x)

            # RAW SPECTROGRAM
            mel_spec = librosa.feature.melspectrogram(y=x, sr=200, hop_length=len(x)//300, 
                  n_fft=1024, n_mels=100, fmin=0, fmax=20, win_length=128)
            
            # LOG TRANSFORM
            width = (mel_spec.shape[1]//30)*30
            mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max).astype(np.float32)[:,:width]
            img[:,:,k] += mel_spec_db
                
        # AVERAGE THE 4 MONTAGE DIFFERENCES
        img[:,:,k] /= 4.0
        
        if display:
            plt.subplot(2,2,k+1)
            plt.imshow(img[:,:,k],aspect='auto',origin='lower')
            
    if display: 
        plt.show()
        plt.figure(figsize=(10,5))
        offset = 0
        for k in range(4):
            if k>0: offset -= signals[3-k].min()
            plt.plot(range(10_000),signals[k]+offset,label=NAMES[3-k])
            offset += signals[3-k].max()
        plt.legend()
        plt.show()
        
    return img
  1. USE_WAVELET = None
    WaveLet変換によるノイズ除去の使用有無。
  2. NAMES = ['LL','LP','RP','RR']
    脳波チェインの名称
  3. FEATS[[...]]
    電極場所の名称
  4. def maddest(d, axis=None)
    平均絶対偏差(Mean Absolute Deviation, mad)を計算する関数
  5. return np.mean(np.absolute(d - np.mean(d, axis)), axis)
    与えられたデータセットdに対して、madを計算する。madはデータの散らばり具合を示す統計量で、データセットからその平均値を引いた絶対値の平均である。
  6. def denoise(x, wavelet='haar', level=1):
    haarのマザーウェーブレットを使用してノイズ除去を行う関数。levelはウェーブレット係数で、大きいほど多くの周波数帯域に分割され、周波数特性を捉えられるが、計算コストと過学習のリスクが増加する。例えば16点のデータをlevel1で変換すると
  7. coeff = pywt.wavedec(x, wavelet, mode="per")
    ウェーブレット変換を行い、データを周波数成分に分解
  8. uthresh = sigma * np.sqrt(2*np.log(len(x)))
    ノイズを除去するためのしきい値を計算
  9. coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])
    しきい値以下の小さい係数を落とす(ノイズ除去)
  10. ret=pywt.waverec(coeff, wavelet, mode='per')
    元のデータに戻す
  11. return ret
    ノイズ除去後のデータを返す
  12. def spectrogram_from_eeg(parquet_path, display=False):
    eegからsupectrogramを作成する関数
  13. eeg = pd.read_parquet(parquet_path)
    parquetをdataframeで取得
  14. middle = (len(eeg)-10_000)//2
    10000行を除いた残りを2等分
  15. eeg = eeg.iloc[middle:middle+10_000]
    中心の10000行を抽出
  16. img = np.zeros((100,300,4),dtype='float32')
    値0のメル周波数の分割数100,時間の分割数300,チャネル4の4枚の画像を定義
  17. signals = []
    前処理後のeeg信号を格納するリスト
  18. for k in range(4):
    4回繰り返す
  19. COLS = FEATS[k]
    使用する電極を格納した配列を順番に選択
  20. for kk in range(4):
    さらに4回繰り返す
  21. x1 = eeg[COLS[kk]].values \ x2 = eeg[COLS[kk+1]].values
    電極配列からkk番目とkk+1番目を抽出。その電極のデータをeegから取得してx1とx2に格納
  22. m = np.nanmean(x1)
    kk番目の電極データ(x1)の欠損値を無視して平均を取得
  23. if np.isnan(x1).mean()<1: x1 = np.nan_to_num(x1,nan=m)
    全データが欠損値でなければ、欠損値を平均値で補完。
  24. else: x1[:] = 0
    全データ欠損なら初期化。x2も同様
  25. x = x1 - x2
    欠損値を処理した2つの電極データの差を取得。これを解析対象の信号とする
  26. if USE_WAVELET: \ x = denoise(x, wavelet=USE_WAVELET)
    Wavelet変換を用いてノイズ除去
  27. signals.append(x)
    前処理後のeeg信号を追加。可視化用
  28. mel_spec = librosa.feature.melspectrogram(y=x, sr=200, hop_length=len(x)//300, n_fft=1024, n_mels=100, fmin=0, fmax=20, win_length=128)
    xからスペクトログラムを取得。
    • 引数
      y: 解析対象
      sr: サンプルレート(Hz)。yが1秒に何回サンプリングされたかを渡す。今回は1秒に200サンプル
      hop_length: STFT(短時間フーリエ変換)を行う際にどの程度窓関数をずらすか。今回は入力信号の1/300にしているため、出力データの時間幅も300となる。STFTは解析範囲が重複するフーリエ変換。
      n_fft: FFT(高速フーリエ変換)の周波数分解能を決定する。1024なので、データを1024種類の周波数に分解する。ただしパワースペクトルは対照的なので実際には513個の周波数成分のパワーが取得できる
      n_mels: メルフィルタバンクの数。結果として周波数軸の数
      fmin: 解析周波数の最低値
      fmax: 解析周波数の最大値。20はeegなどの低周波数データに対して使用される
      win_legth: STFTの窓の幅。128データポイントごとに切り取ってフーリエ変換する
  29. width = (mel_spec.shape[1]//30)*30
    メルスペクトログラムの横幅を30で整数除算、30をかけることによって30の倍数にして扱いやすい値を求めておく
  30. mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max).astype(np.float32)[:,:width]
    メルスペクトログラムをデシベルスケールに変換して、幅の末尾を切り取って調整
  31. img[:,:,k] += mel_spec_db
    得られた二次元データをk番目のチャネルに追加。一つのkに対してkkが4回繰り返されるため、FEATSの5つの電極の、それぞれの差を取ったものを4つ足す形になる。
  32. img[:,:,k] /= 4.0
    5つの電極間の差を平均する。一連の電極の流れの中の電極間データの差を表す。
  33. return img
    4cheinの電極配列における、電極と電極の間の差を表した画像を4枚返す。サイズは(100,300,4)(メル周波数、時間(最大10000のステップ10000/300)、各電極チェイン)

2.8 推論

・推論用データの確認

if submission:
    test = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/test.csv')
    print('Test shape',test.shape)
    test.head()
  1. if submission: \ test = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/test.csv')
    推論用データの確認

・推論用スペクトログラムデータのロード

# READ ALL SPECTROGRAMS
if submission:
    PATH2 = '/kaggle/input/hms-harmful-brain-activity-classification/test_spectrograms/'
    files2 = os.listdir(PATH2)
    print(f'There are {len(files2)} test spectrogram parquets')
    
    spectrograms2 = {}
    for i,f in enumerate(files2):
        if i%100==0: print(i,', ',end='')
        tmp = pd.read_parquet(f'{PATH2}{f}')
        name = int(f.split('.')[0])
        spectrograms2[name] = tmp.iloc[:,1:].values
    
    # RENAME FOR DATA GENERATOR
    test = test.rename({'spectrogram_id':'spec_id'},axis=1)
  1. files2 = os.listdir(PATH2)
    推論用スペクトログラムデータのファイル名を配列で取得
  2. spectrograms2 = {}
    スペクトログラムデータ補完用の辞書を定義
  3. tmp = pd.read_parquet(f'{PATH2}{f}')
    スペクトログラムデータをdataframeで取得
  4. name = int(f.split('.')[0])
    スペクトログラムデータのファイル名を取得
  5. spectrograms2[name] = tmp.iloc[:,1:].values
    {ファイル名:スペクトログラムのDataframe}の辞書を作成。ilocは時間方向の余分な情報の可能性のある最初の行を削除。(一般的にメタデータなどの可能性有)
  6. test = test.rename({'spectrogram_id':'spec_id'},axis=1)
    推論用データ(test)の'spectrogram_id'列を'spec_id'に名称変更

・推論用eegデータをスペクトログラムにしてロード

# READ ALL EEG SPECTROGRAMS
if submission:
    PATH2 = '/kaggle/input/hms-harmful-brain-activity-classification/test_eegs/'
    DISPLAY = 0
    EEG_IDS2 = test.eeg_id.unique()
    all_eegs2 = {}

    print('Converting Test EEG to Spectrograms...'); print()
    for i,eeg_id in enumerate(EEG_IDS2):
        
        # CREATE SPECTROGRAM FROM EEG PARQUET
        img = spectrogram_from_eeg(f'{PATH2}{eeg_id}.parquet', i<DISPLAY)
        all_eegs2[eeg_id] = img
  1. EEG_IDS2 = test.eeg_id.unique()
    EEG_IDS2として、推論用データから一意のeeg_idを抜き出す。
  2. all_eegs2 = {}
    eeg格納用の辞書を定義
  3. img = spectrogram_from_eeg(f'{PATH2}{eeg_id}.parquet', i<DISPLAY)
    parquet形式のeegをスペクトログラムに変換して抜き出す
  4. all_eegs2[eeg_id] = img
    {eeg_id: スペクトログラム}の辞書を作成

・非アンサンブル推論

# Submission ON TEST without ensemble
if submission and not ENSEMBLE:
    preds = []
    
    if DATA_TYPE=='raw':
        test_gen = DataGenerator(test, mode='test', raw_eegs=all_raw_eegs2)
        in_shape = (2000,8)
    else:
        test_gen = DataGenerator(test, mode='test', specs = spectrograms2, eeg_specs = all_eegs2)
        in_shape = (512,512,3)
    
    test_dataset = tf.data.Dataset.from_generator(generator=test_gen, 
                                               output_signature=(tf.TensorSpec(shape=in_shape, dtype=tf.float32),
                                                                 tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)
    if DATA_TYPE=='raw':
        model = build_wave_model()
    else:
        model = build_model()

    for i in range(5):
        print(f'Fold {i+1}')
        model.load_weights(f'{LOAD_MODELS_FROM}model_{DATA_TYPE}_{VER}_{i}.weights.h5')
        pred = model.predict(test_dataset, verbose=1)
        preds.append(pred)
        
    pred = np.mean(preds,axis=0)
    print('Test preds shape',pred.shape)
  1. preds = []
    推論結果格納用リスト
  2. if DATA_TYPE=='raw':
    生のeegデータをWaveNetで利用する場合
  3. test_gen = DataGenerator(test, mode='test', raw_eegs=all_raw_eegs2)
    WaveNet用データジェネレーターの定義
  4. in_shape = (2000,8)
    入力データ形状を定義
  5. else:
    スペクトログラムデータをEfficientNetで利用する場合
  6. test_gen = DataGenerator(test, mode='test', specs = spectrograms2, eeg_specs = all_eegs2)
    EffNet用データジェネレーターの定義
  7. in_shape = (512,512,3)
    入力データ形状を定義
  8. test_dataset = tf.data.Dataset.from_generator(generator=test_gen, output_signature=(tf.TensorSpec(shape=in_shape, dtype=tf.float32), tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)
    ・引数
    generator: 定義したジェネレータ
    output_signature: 入力形状と出力形状を渡す
    .batch: バッチサイズを渡す
    .prefetch: GPUの計算リソース活用のために、学習中に次のデータを事前に読み込む際の挙動を指定。AUTOTUNEにするとバッファサイズを調整して、最適なパフォーマンスを自動的に決定する
  9. if DATA_TYPE=='raw': \ model = build_wave_model() \ else: \ model = build_model()
    使用データに合わせてモデルを定義
  10. for i in range(5):
    foldの回数繰り返し
  11. model.load_weights(f'{LOAD_MODELS_FROM}model_{DATA_TYPE}{VER}{i}.weights.h5')
    モデルの重みをロード
  12. pred = model.predict(test_dataset, verbose=1)
    推論用データセットを使用して推論
  13. preds.append(pred)
    予測結果をリストに格納
  14. pred = np.mean(preds,axis=0)
    列方向に平均。5foldの結果を平均する

・アンサンブルを使用した推論

# Submission ON TEST with ensemble
if submission and ENSEMBLE:
    preds = []
    test_gen_kaggle = DataGenerator(test, mode='test', data_type='kaggle', specs = spectrograms2, eeg_specs = all_eegs2)
    test_dataset_kaggle = tf.data.Dataset.from_generator(generator=test_gen_kaggle, 
                                               output_signature=(tf.TensorSpec(shape=(512,512,3), dtype=tf.float32),
                                                                 tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)
    test_gen_both = DataGenerator(test, mode='test', data_type='both', specs = spectrograms2, eeg_specs = all_eegs2)
    test_dataset_both = tf.data.Dataset.from_generator(generator=test_gen_both, 
                                               output_signature=(tf.TensorSpec(shape=(512,512,3), dtype=tf.float32),
                                                                 tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)

    test_gen_eeg = DataGenerator(test, mode='test', data_type='eeg', specs = spectrograms2, eeg_specs = all_eegs2)
    test_dataset_eeg = tf.data.Dataset.from_generator(generator=test_gen_eeg, 
                                               output_signature=(tf.TensorSpec(shape=(512,512,3), dtype=tf.float32),
                                                                 tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)
    test_gen_raw = DataGenerator(test, mode='test', data_type='raw', raw_eegs=all_raw_eegs2)
    test_dataset_raw = tf.data.Dataset.from_generator(generator=test_gen_raw, 
                                               output_signature=(tf.TensorSpec(shape=(2000,8), dtype=tf.float32),
                                                                 tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)
 
    # LB SCORE FOR EACH MODEL
    lbs = 1 - np.array(LBs)
    weights = lbs/lbs.sum()
    model = build_model()
    model_wave = build_wave_model()

    for i in range(5):
        print(f'Fold {i+1}')
        
        model.load_weights(f'{LOAD_MODELS_FROM}model_kaggle_{VERK}_{i}.weights.h5')
        pred_kaggle = model.predict(test_dataset_kaggle, verbose=1)
        
        model.load_weights(f'{LOAD_MODELS_FROM}model_both_{VERB}_{i}.weights.h5')
        pred_both = model.predict(test_dataset_both, verbose=1)
        
        model.load_weights(f'{LOAD_MODELS_FROM}model_eeg_{VERE}_{i}.weights.h5')
        pred_eeg = model.predict(test_dataset_eeg, verbose=1)
        
        model_wave.load_weights(f'{LOAD_MODELS_FROM}model_raw_{VERR}_{i}.weights.h5')
        pred_raw = model_wave.predict(test_dataset_raw, verbose=1)
        
        pred = np.array([pred_both,pred_eeg,pred_kaggle,pred_raw])
        pred = np.average(pred,axis=0,weights=weights)
        preds.append(pred)
        
    pred = np.mean(preds,axis=0)
    print('Test preds shape',pred.shape)
  1. preds = []
    予測結果格納用リスト
  2. test_gen_kaggle = DataGenerator(test, mode='test', data_type='kaggle', specs = spectrograms2, eeg_specs = all_eegs2)
    kaggleから提供されたスペクトログラムのデータジェネレータを定義
  3. test_dataset_kaggle = tf.data.Dataset.from_generator(generator=test_gen_kaggle, output_signature=(tf.TensorSpec(shape=(512,512,3), dtype=tf.float32), tf.TensorSpec(shape=(6,), dtype=tf.float32))).batch(64).prefetch(tf.data.AUTOTUNE)
    kaggleスペクトログラムのデータセット定義
  4. test_gen_both \ test_gen_eeg \ test_gen_raw
    上記同様にデータジェネレータとデータセットを定義
    both: kaggleスペクトログラムとeegから生成したスペクトログラムの両方を使用
    eeg: eegから生成したスペクトログラムを使用
    raw: eegデータを使用
  5. lbs = 1 - np.array(LBs)
    それぞれのモデルのLBスコア配列を逆転(大きいほど良いモデルとする)
  6. weights = lbs/lbs.sum()
    逆転したLBスコア配列を元にモデルをアンサンブルする際の重み(最終推論への影響度)を決定
  7. model = build_model() \ model_wave = build_wave_model()
    使用モデルの定義
  8. for i in range(5):
    fold数繰り返し
  9. model.load_weights \ model.predict
    3種類のスペクトログラムデータセット(kaggle, eeg, both)とEfficientNetの対応する重みを使用して推論を行う
  10. model_wave.load_weights \ model_wave.predict
    生のeegデータセットとWaveNetの重みを使用して推論を行う
  11. pred = np.array([pred_both,pred_eeg,pred_kaggle,pred_raw])
    各モデルの推論結果を二次元配列で格納
  12. pred = np.average(pred,axis=0,weights=weights)
    モデルのLBスコアによる重みを利用して平均を取る
  13. preds.append(pred)
    4モデルの推論結果を平均したものを追加。5flod分
  14. pred = np.mean(preds,axis=0)
    5fold分の推論結果を列方向に平均

2.9 提出

if submission:
    sub = pd.DataFrame({'eeg_id':test.eeg_id.values})
    sub[TARGETS] = pred
    sub.to_csv('submission.csv',index=False)
    print('Submissionn shape',sub.shape)
    print()
    print(sub.head().to_string())
  1. sub = pd.DataFrame({'eeg_id':test.eeg_id.values})
    列名'eeg_id'として推論用eeg_idを格納したデータフレームsubを定義
  2. sub[TARGETS] = pred
    推論結果をsubの[TRAGET]列として追加
  3. sub.to_csv('submission.csv',index=False)
    csvファイルに変更

・最終確認

if submission:
    print(sub.iloc[:,-6:].sum(axis=1).to_string())
  1. print(sub.iloc[:,-6:].sum(axis=1).to_string())
    推論結果の合計が1になることを確認

試したいこと
・スペクトログラム画像の三チャネル目使ってないのなんで?有効活用できないかな
・2.6の7のkl距離でフィルタかけてるやつ、予測したいのはあくまで専門家の評価やから、「評価が割れてるデータ」も学習する価値があるように思う。のでフィルター消してみる→ちゃうわ、これフィルタかけてるやつ(gen2)とかけてないやつ(gen)両方使ってる。
・他ノートブックのspectrogram_from_eeg関数を使ってみる

Discussion