{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":59093,"databundleVersionId":7469972,"sourceType":"competition"},{"sourceId":7392733,"sourceType":"datasetVersion","datasetId":4297749},{"sourceId":7414022,"sourceType":"datasetVersion","datasetId":4312784}],"dockerImageVersionId":30635,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.12"},"papermill":{"default_parameters":{},"duration":102.485133,"end_time":"2024-01-16T12:14:21.10455","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-01-16T12:12:38.619417","version":"2.4.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd, numpy as np, os\nimport random\nimport matplotlib.pyplot as plt, gc\n\ntrain = pd.read_csv('/kaggle/input/hms-harmful-brain-activity-classification/train.csv')\nprint('Train shape', train.shape )\ndisplay( train.head() )\n\nCREATE_SPECTROGRAMS = True","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":0.693827,"end_time":"2024-01-16T12:12:42.606147","exception":false,"start_time":"2024-01-16T12:12:41.91232","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-25T02:46:12.373302Z","iopub.execute_input":"2024-01-25T02:46:12.373979Z","iopub.status.idle":"2024-01-25T02:46:13.140453Z","shell.execute_reply.started":"2024-01-25T02:46:12.373942Z","shell.execute_reply":"2024-01-25T02:46:13.139321Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### v1\nLL = ( (Fp1 - F7) + (F7 - T3) + (T3 - T5) + (T5 - O1) )/4.\n### v2\nfrom https://www.kaggle.com/code/cdeotte/how-to-make-spectrogram-from-eeg v4\n\nLL Spec = ( spec(Fp1 - F7) + spec(F7 - T3) + spec(T3 - T5) + spec(T5 - O1) )/4.\n\n","metadata":{}},{"cell_type":"code","source":"NAMES = ['LL','LP','RP','RR']\n\nFEATS = [['Fp1','F7','T3','T5','O1'],\n         ['Fp1','F3','C3','P3','O1'],\n         ['Fp2','F8','T4','T6','O2'],\n         ['Fp2','F4','C4','P4','O2']]\n\ndirectory_path = 'EEG_Spectrograms/'\nif not os.path.exists(directory_path):\n    os.makedirs(directory_path)","metadata":{"papermill":{"duration":0.017131,"end_time":"2024-01-16T12:12:42.644893","exception":false,"start_time":"2024-01-16T12:12:42.627762","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-25T02:46:13.142800Z","iopub.execute_input":"2024-01-25T02:46:13.143582Z","iopub.status.idle":"2024-01-25T02:46:13.152318Z","shell.execute_reply.started":"2024-01-25T02:46:13.143537Z","shell.execute_reply":"2024-01-25T02:46:13.150950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pywt\nprint(\"The wavelet functions we can use:\")\nprint(pywt.wavelist())\n\nUSE_WAVELET = None #or \"db8\" or anything below","metadata":{"papermill":{"duration":0.540601,"end_time":"2024-01-16T12:12:43.206759","exception":false,"start_time":"2024-01-16T12:12:42.666158","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-25T02:46:13.153876Z","iopub.execute_input":"2024-01-25T02:46:13.154374Z","iopub.status.idle":"2024-01-25T02:46:13.697146Z","shell.execute_reply.started":"2024-01-25T02:46:13.154333Z","shell.execute_reply":"2024-01-25T02:46:13.695985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# DENOISE FUNCTION\ndef maddest(d, axis=None):\n    return np.mean(np.absolute(d - np.mean(d, axis)), axis)\n\ndef denoise(x, wavelet='haar', level=1):    \n    coeff = pywt.wavedec(x, wavelet, mode=\"per\")\n    sigma = (1/0.6745) * maddest(coeff[-level])\n\n    uthresh = sigma * np.sqrt(2*np.log(len(x)))\n    coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])\n\n    ret=pywt.waverec(coeff, wavelet, mode='per')\n    \n    return ret","metadata":{"papermill":{"duration":0.018392,"end_time":"2024-01-16T12:12:43.23335","exception":false,"start_time":"2024-01-16T12:12:43.214958","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-01-25T02:46:13.700188Z","iopub.execute_input":"2024-01-25T02:46:13.700975Z","iopub.status.idle":"2024-01-25T02:46:13.711600Z","shell.execute_reply.started":"2024-01-25T02:46:13.700926Z","shell.execute_reply":"2024-01-25T02:46:13.710364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# From https://github.com/tomrunia/PyTorchWavelets/blob/master/wavelets_pytorch/wavelets.py\n\n## from G2net1 \nimport torch\nfrom scipy import signal\nfrom scipy import optimize\nimport torch.nn as nn\nfrom timm.layers.conv2d_same import conv2d_same\n\nclass Morlet(object):\n    def __init__(self, w0=6):\n        \"\"\"w0 is the nondimensional frequency constant. If this is\n        set too low then the wavelet does not sample very well: a\n        value over 5 should be ok; Terrence and Compo set it to 6.\n        \"\"\"\n        self.w0 = w0\n        if w0 == 6:\n            # value of C_d from TC98\n            self.C_d = 0.776\n\n    def __call__(self, *args, **kwargs):\n        return self.time(*args, **kwargs)\n\n    def time(self, t, s=1.0, complete=True):\n        \"\"\"\n        Complex Morlet wavelet, centred at zero.\n        Parameters\n        ----------\n        t : float\n            Time. If s is not specified, this can be used as the\n            non-dimensional time t/s.\n        s : float\n            Scaling factor. Default is 1.\n        complete : bool\n            Whether to use the complete or the standard version.\n        Returns\n        -------\n        out : complex\n            Value of the Morlet wavelet at the given time\n        See Also\n        --------\n        scipy.signal.gausspulse\n        Notes\n        -----\n        The standard version::\n            pi**-0.25 * exp(1j*w*x) * exp(-0.5*(x**2))\n        This commonly used wavelet is often referred to simply as the\n        Morlet wavelet.  Note that this simplified version can cause\n        admissibility problems at low values of `w`.\n        The complete version::\n            pi**-0.25 * (exp(1j*w*x) - exp(-0.5*(w**2))) * exp(-0.5*(x**2))\n        The complete version of the Morlet wavelet, with a correction\n        term to improve admissibility. For `w` greater than 5, the\n        correction term is negligible.\n        Note that the energy of the return wavelet is not normalised\n        according to `s`.\n        The fundamental frequency of this wavelet in Hz is given\n        by ``f = 2*s*w*r / M`` where r is the sampling rate.\n        \"\"\"\n        w = self.w0\n\n        x = t / s\n\n        output = np.exp(1j * w * x)\n\n        if complete:\n            output -= np.exp(-0.5 * (w ** 2))\n\n        output *= np.exp(-0.5 * (x ** 2)) * np.pi ** (-0.25)\n\n        return output\n\n    # Fourier wavelengths\n    def fourier_period(self, s):\n        \"\"\"Equivalent Fourier period of Morlet\"\"\"\n        return 4 * np.pi * s / (self.w0 + (2 + self.w0 ** 2) ** 0.5)\n\n    def scale_from_period(self, period):\n        \"\"\"\n        Compute the scale from the fourier period.\n        Returns the scale\n        \"\"\"\n        # Solve 4 * np.pi * scale / (w0 + (2 + w0 ** 2) ** .5)\n        #  for s to obtain this formula\n        coeff = np.sqrt(self.w0 * self.w0 + 2)\n        return (period * (coeff + self.w0)) / (4.0 * np.pi)\n\n    # Frequency representation\n    def frequency(self, w, s=1.0):\n        \"\"\"Frequency representation of Morlet.\n        Parameters\n        ----------\n        w : float\n            Angular frequency. If `s` is not specified, i.e. set to 1,\n            this can be used as the non-dimensional angular\n            frequency w * s.\n        s : float\n            Scaling factor. Default is 1.\n        Returns\n        -------\n        out : complex\n            Value of the Morlet wavelet at the given frequency\n        \"\"\"\n        x = w * s\n        # Heaviside mock\n        Hw = np.array(w)\n        Hw[w <= 0] = 0\n        Hw[w > 0] = 1\n        return np.pi ** -0.25 * Hw * np.exp((-((x - self.w0) ** 2)) / 2)\n\n    def coi(self, s):\n        \"\"\"The e folding time for the autocorrelation of wavelet\n        power at each scale, i.e. the timescale over which an edge\n        effect decays by a factor of 1/e^2.\n        This can be worked out analytically by solving\n            |Y_0(T)|^2 / |Y_0(0)|^2 = 1 / e^2\n        \"\"\"\n        return 2 ** 0.5 * s\n\n\nclass CWT(nn.Module):\n    def __init__(\n        self,\n        dj=0.0625,\n        dt=1 / 200,\n        wavelet=Morlet(),\n        fmin: int = 20,\n        fmax: int = 500,\n        output_format=\"Magnitude\",\n        trainable=False,\n        hop_length: int = 1,\n    ):\n        super().__init__()\n        self.wavelet = wavelet\n\n        self.dt = dt\n        self.dj = dj\n        self.fmin = fmin\n        self.fmax = fmax\n        self.output_format = output_format\n        self.trainable = trainable  # TODO make kernel a trainable parameter\n        self.stride = (1, hop_length)\n        # self.padding = 0  # \"same\"\n\n        self._scale_minimum = self.compute_minimum_scale()\n\n        self.signal_length = None\n        self._channels = None\n\n        self._scales = None\n        self._kernel = None\n        self._kernel_real = None\n        self._kernel_imag = None\n\n    def compute_optimal_scales(self):\n        \"\"\"\n        Determines the optimal scale distribution (see. Torrence & Combo, Eq. 9-10).\n        :return: np.ndarray, collection of scales\n        \"\"\"\n        if self.signal_length is None:\n            raise ValueError(\n                \"Please specify signal_length before computing optimal scales.\"\n            )\n        J = int(\n            (1 / self.dj) * np.log2(self.signal_length * self.dt / self._scale_minimum)\n        )\n        scales = self._scale_minimum * 2 ** (self.dj * np.arange(0, J + 1))\n\n        # Remove high and low frequencies\n        frequencies = np.array([1 / self.wavelet.fourier_period(s) for s in scales])\n        if self.fmin:\n            frequencies = frequencies[frequencies >= self.fmin]\n            scales = scales[0 : len(frequencies)]\n        if self.fmax:\n            frequencies = frequencies[frequencies <= self.fmax]\n            scales = scales[len(scales) - len(frequencies) : len(scales)]\n\n        return scales\n\n    def compute_minimum_scale(self):\n        \"\"\"\n        Choose s0 so that the equivalent Fourier period is 2 * dt.\n        See Torrence & Combo Sections 3f and 3h.\n        :return: float, minimum scale level\n        \"\"\"\n        dt = self.dt\n\n        def func_to_solve(s):\n            return self.wavelet.fourier_period(s) - 2 * dt\n\n        return optimize.fsolve(func_to_solve, 1)[0]\n\n    def _build_filters(self):\n        self._filters = []\n        for scale_idx, scale in enumerate(self._scales):\n            # Number of points needed to capture wavelet\n            M = 10 * scale / self.dt\n            # Times to use, centred at zero\n            t = torch.arange((-M + 1) / 2.0, (M + 1) / 2.0) * self.dt\n            if len(t) % 2 == 0:\n                t = t[0:-1]  # requires odd filter size\n            # Sample wavelet and normalise\n            norm = (self.dt / scale) ** 0.5\n            filter_ = norm * self.wavelet(t, scale)\n            self._filters.append(torch.conj(torch.flip(filter_, [-1])))\n\n        self._pad_filters()\n\n    def _pad_filters(self):\n        filter_len = self._filters[-1].shape[0]\n        padded_filters = []\n\n        for f in self._filters:\n            pad = (filter_len - f.shape[0]) // 2\n            padded_filters.append(nn.functional.pad(f, (pad, pad)))\n\n        self._filters = padded_filters\n\n    def _build_wavelet_bank(self):\n        \"\"\"This function builds a 2D wavelet filter using wavelets at different scales\n\n        Returns:\n            tensor: Tensor of shape (num_widths, 1, channels, filter_len)\n        \"\"\"\n        self._build_filters()\n        wavelet_bank = torch.stack(self._filters)\n        wavelet_bank = wavelet_bank.view(\n            wavelet_bank.shape[0], 1, 1, wavelet_bank.shape[1]\n        )\n        # See comment by tez6c32\n        # https://www.kaggle.com/anjum48/continuous-wavelet-transform-cwt-in-pytorch/comments#1499878\n        # wavelet_bank = torch.cat([wavelet_bank] * self.channels, 2)\n        return wavelet_bank\n\n    def forward(self, x):\n        \"\"\"Compute CWT arrays from a batch of multi-channel inputs\n\n        Args:\n            x (torch.tensor): Tensor of shape (batch_size, channels, time)\n\n        Returns:\n            torch.tensor: Tensor of shape (batch_size, channels, widths, time)\n        \"\"\"\n        if self.signal_length is None:\n            self.signal_length = x.shape[-1]\n            self.channels = x.shape[-2]\n            self._scales = self.compute_optimal_scales()\n            self._kernel = self._build_wavelet_bank()\n\n            if self._kernel.is_complex():\n                self._kernel_real = self._kernel.real\n                self._kernel_imag = self._kernel.imag\n\n        x = x.unsqueeze(1)\n        if self._kernel.is_complex():\n            if (\n                x.dtype != self._kernel_real.dtype\n                or x.device != self._kernel_real.device\n            ):\n                self._kernel_real = self._kernel_real.to(device=x.device, dtype=x.dtype)\n                self._kernel_imag = self._kernel_imag.to(device=x.device, dtype=x.dtype)\n\n            # Strides > 1 not yet supported for \"same\" padding\n            # output_real = nn.functional.conv2d(\n            #     x, self._kernel_real, padding=self.padding, stride=self.stride\n            # )\n            # output_imag = nn.functional.conv2d(\n            #     x, self._kernel_imag, padding=self.padding, stride=self.stride\n            # )\n            output_real = conv2d_same(x, self._kernel_real, stride=self.stride)\n            output_imag = conv2d_same(x, self._kernel_imag, stride=self.stride)\n            output_real = torch.transpose(output_real, 1, 2)\n            output_imag = torch.transpose(output_imag, 1, 2)\n\n            if self.output_format == \"Magnitude\":\n                return torch.sqrt(output_real ** 2 + output_imag ** 2)\n            else:\n                return torch.stack([output_real, output_imag], -1)\n\n        else:\n            if x.device != self._kernel.device:\n                self._kernel = self._kernel.to(device=x.device, dtype=x.dtype)\n\n            # output = nn.functional.conv2d(\n            #     x, self._kernel, padding=self.padding, stride=self.stride\n            # )\n            output = conv2d_same(x, self._kernel, stride=self.stride)\n            return torch.transpose(output, 1, 2)","metadata":{"execution":{"iopub.status.busy":"2024-01-25T02:46:13.713625Z","iopub.execute_input":"2024-01-25T02:46:13.714359Z","iopub.status.idle":"2024-01-25T02:46:21.346440Z","shell.execute_reply.started":"2024-01-25T02:46:13.714314Z","shell.execute_reply":"2024-01-25T02:46:21.344655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(\"/kaggle/input/hms-harmful-brain-activity-classification/train.csv\")\n\n\n","metadata":{"execution":{"iopub.status.busy":"2024-01-25T02:46:21.349853Z","iopub.execute_input":"2024-01-25T02:46:21.351190Z","iopub.status.idle":"2024-01-25T02:46:21.587785Z","shell.execute_reply.started":"2024-01-25T02:46:21.351136Z","shell.execute_reply":"2024-01-25T02:46:21.586752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport librosa\nfor i in range(6):\n    df_ = df[df[\"expert_consensus\"]==df[\"expert_consensus\"].unique()[i]]\n    print(\"LABEL==\",df_[\"expert_consensus\"].unique()[0])\n    for j in random.sample(range(len(df_)),5):\n        row = df_.iloc[j]\n\n        sp = pd.read_parquet(f\"/kaggle/input/hms-harmful-brain-activity-classification/train_spectrograms/{row.spectrogram_id}.parquet\")\n        spec_offset = int( row.spectrogram_label_offset_seconds )\n        sp = sp.loc[(sp.time>=spec_offset)\n                             &(sp.time<spec_offset+600)]\n        img = sp.fillna(0).values[:, 1:].T.astype(\"float32\") \n        img = np.clip(img,np.exp(-4), np.exp(8))\n        img = np.log(img)\n\n        # normalize per image\n        eps = 1e-6\n        img_mean = img.mean(axis=(0, 1))\n        img = img - img_mean\n        img_std = img.std(axis=(0, 1))\n        img = img / (img_std + eps)\n        plt.figure()\n\n        plt.subplot(1, 4, 1)\n        plt.imshow(img, aspect=\"auto\")\n        plt.title(\"kaggle spectrogram\")\n        plt.axis('off')  \n\n        eeg = pd.read_parquet(f\"/kaggle/input/hms-harmful-brain-activity-classification/train_eegs/{row.eeg_id}.parquet\")\n        #middle = (len(eeg)-10_000)//2\n        #eeg = eeg.iloc[middle:middle+10_000]\n        eeg_offset = int( row.eeg_label_offset_seconds )\n        eeg = eeg.iloc[eeg_offset*200:(eeg_offset+50)*200]\n        signals = []\n        img = np.zeros((128,256,4),dtype='float32')\n        cwt_img = np.zeros((166,527,4),dtype='float32')\n        cwt_img1 = np.zeros((166,527,4),dtype='float32')\n\n        pycwt = CWT(fmin=0, fmax=25, hop_length=10000//512)\n\n        signals = []\n        for k in range(4):\n            COLS = FEATS[k]\n\n            for kk in range(4):\n\n                # COMPUTE PAIR DIFFERENCES\n                x = eeg[COLS[kk]].values - eeg[COLS[kk+1]].values\n\n                # FILL NANS\n                m = np.nanmean(x)\n                if np.isnan(x).mean()<1: x = np.nan_to_num(x,nan=m)\n                else: x[:] = 0\n\n                # DENOISE\n                if USE_WAVELET:\n                    x = denoise(x, wavelet=USE_WAVELET)\n                #signals.append(x)\n                S = torch.tensor(x)[None,:]\n\n                # RAW SPECTROGRAM\n                mel_spec = librosa.feature.melspectrogram(y=x, sr=200, hop_length=len(x)//256, \n                      n_fft=1024, n_mels=128, fmin=0, fmax=20, win_length=128)\n\n                # LOG TRANSFORM\n                width = (mel_spec.shape[1]//32)*32\n                mel_spec_db = librosa.power_to_db(mel_spec, ref=np.max).astype(np.float32)[:,:width]\n\n                # STANDARDIZE TO -1 TO 1\n                mel_spec_db = (mel_spec_db+40)/40 \n                img[:,:,k] += mel_spec_db\n                \n                out = pycwt(S).numpy()\n                cwt_img[:,:,k] +=out[:,:,0]\n                out1 = librosa.power_to_db(out, ref=np.max).astype(np.float32)\n                out1 = (out1+40)/40 \n                cwt_img1[:,:,k] +=out1[:,:,0]\n\n            # AVERAGE THE 4 MONTAGE DIFFERENCES\n            img[:,:,k] /= 4.0\n            cwt_img[:,:,k] /= 4.0\n            cwt_img1[:,:,k] /= 4.0\n\n        \n        \n        cwt_spec = np.concatenate([cwt_img[:,:,0],cwt_img[:,:,1],cwt_img[:,:,2],cwt_img[:,:,3]])\n        plt.subplot(1,4,2)\n        plt.imshow(cwt_spec, aspect=\"auto\")\n        plt.title(\"eeg → CWT\")\n        plt.axis('off')  \n        cwt_spec = np.concatenate([cwt_img1[:,:,0],cwt_img1[:,:,1],cwt_img1[:,:,2],cwt_img1[:,:,3]])\n        plt.subplot(1,4,3)\n        plt.imshow(cwt_spec, aspect=\"auto\")\n        plt.title(\"eeg → CWT norm\")\n        plt.axis('off')  \n\n    \n\n        spec = np.concatenate([img[:,:,0],img[:,:,1],img[:,:,2],img[:,:,3]])\n        plt.subplot(1, 4, 4)\n        plt.imshow(spec, aspect=\"auto\")\n        plt.title(\"eeg → STFT\")\n        plt.axis('off')  \n        plt.tight_layout()\n\n        plt.show()\n    ","metadata":{"execution":{"iopub.status.busy":"2024-01-25T02:46:21.589579Z","iopub.execute_input":"2024-01-25T02:46:21.590548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"CWTs seem to be able to extract different information from STFTs, but my local CV has not improved...","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}