SoundDr_cough.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook to create a overall system, include Unsupervised(Isolation Forest, XGBOD) at last this notebook.\n",
    "\n",
    "+ Choose dataset ```dataset_type='SoundDr'```:\n",
    "+ Choose pretrain model to extract Feature ```PRETRAIN=\"FRILL\"```\n",
    "+ Choose model to classify ```Classifier=\"SVM\"```\n",
    "+ Set seed ```seed=\"2022\"```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "SoX could not be found!\n",
      "\n",
      "    If you do not have SoX, proceed here:\n",
      "     - - - http://sox.sourceforge.net/ - - -\n",
      "\n",
      "    If you do (or think that you should) have SoX, double-check your\n",
      "    path variables.\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "import warnings\n",
    "from glob import glob\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import os, time, math, random, cv2\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import zipfile, pickle, h5py, joblib, json\n",
    "import multiprocessing\n",
    "\n",
    "import librosa\n",
    "import opensmile\n",
    "import xgboost as xgb\n",
    "# import tensorflow_hub as hub\n",
    "\n",
    "from math import pi\n",
    "from tqdm import tqdm\n",
    "from pathlib import Path\n",
    "from functools import partial\n",
    "from scipy.fftpack import fft, hilbert\n",
    "from scipy.spatial.distance import cdist\n",
    "from scipy.spatial import cKDTree\n",
    "from sklearn.svm import LinearSVC\n",
    "from sklearn.ensemble import IsolationForest\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import StratifiedKFold, cross_val_predict, train_test_split\n",
    "from sklearn.metrics import f1_score, confusion_matrix, roc_auc_score, auc, precision_recall_curve, roc_curve, average_precision_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Configuration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Classifier: Model need to classify. Currently, We support SVM and XGBoost Classify. We show Unsupervised(Isolation Forest, XGBOD) at last this notebook.\n",
    "+ PRETRAIN: Feature type need to extract for classify model. We support these: TRILL, FRILL, OpenSmile, OpenSmileBoAW, DeepSpectrum\n",
    "+ dataset_type: dataset type which we need run on such as SoundDr, CoughVid, Coswara"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sex_choice</th>\n",
       "      <th>age_choice</th>\n",
       "      <th>current_city</th>\n",
       "      <th>symptoms_status_choice</th>\n",
       "      <th>medical_condition_choice</th>\n",
       "      <th>insomnia_status_choice</th>\n",
       "      <th>smoke_status_choice</th>\n",
       "      <th>cov19_status_choice</th>\n",
       "      <th>hospital_choice</th>\n",
       "      <th>cough_noise</th>\n",
       "      <th>device_model</th>\n",
       "      <th>file_cough</th>\n",
       "      <th>label</th>\n",
       "      <th>cough_duration</th>\n",
       "      <th>nose_duration</th>\n",
       "      <th>mouth_duration</th>\n",
       "      <th>file_path</th>\n",
       "      <th>label_symptom</th>\n",
       "      <th>label_abnormal</th>\n",
       "      <th>label_covid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1305</th>\n",
       "      <td>Female</td>\n",
       "      <td>20</td>\n",
       "      <td>thanh hóa</td>\n",
       "      <td>['No']</td>\n",
       "      <td>['No']</td>\n",
       "      <td>No</td>\n",
       "      <td>never</td>\n",
       "      <td>never</td>\n",
       "      <td>No</td>\n",
       "      <td>True</td>\n",
       "      <td>OPPO CPH1933</td>\n",
       "      <td>cough/good_cough_2021-08-15T13:43:33.132Z</td>\n",
       "      <td>0</td>\n",
       "      <td>25.856000</td>\n",
       "      <td>17.664000</td>\n",
       "      <td>22.357333</td>\n",
       "      <td>cough/good_cough_2021-08-15T13:43:33.132Z.wav</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1306</th>\n",
       "      <td>Male</td>\n",
       "      <td>32</td>\n",
       "      <td>Ho Chi Minh</td>\n",
       "      <td>['fever', 'headache']</td>\n",
       "      <td>['No']</td>\n",
       "      <td>No</td>\n",
       "      <td>never</td>\n",
       "      <td>last14</td>\n",
       "      <td>No</td>\n",
       "      <td>True</td>\n",
       "      <td>Laptop/Desktop</td>\n",
       "      <td>cough/bad_cough_2021-09-16T07:18:48.594Z</td>\n",
       "      <td>1</td>\n",
       "      <td>29.350042</td>\n",
       "      <td>29.814438</td>\n",
       "      <td>29.814438</td>\n",
       "      <td>cough/bad_cough_2021-09-16T07:18:48.594Z.wav</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1307</th>\n",
       "      <td>Male</td>\n",
       "      <td>23</td>\n",
       "      <td>Ho Chi Minh</td>\n",
       "      <td>['fever', 'chills', 'sorethroat', 'drycough', ...</td>\n",
       "      <td>['No']</td>\n",
       "      <td>1</td>\n",
       "      <td>1to10</td>\n",
       "      <td>last14</td>\n",
       "      <td>No</td>\n",
       "      <td>True</td>\n",
       "      <td>Laptop/Desktop</td>\n",
       "      <td>cough/bad_cough_2021-09-06T11:04:49.842Z</td>\n",
       "      <td>1</td>\n",
       "      <td>18.432000</td>\n",
       "      <td>17.152000</td>\n",
       "      <td>15.530667</td>\n",
       "      <td>cough/bad_cough_2021-09-06T11:04:49.842Z.wav</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1308</th>\n",
       "      <td>Male</td>\n",
       "      <td>18</td>\n",
       "      <td>Ho Chi Minh</td>\n",
       "      <td>['wetcough', 'sorethroat']</td>\n",
       "      <td>['No']</td>\n",
       "      <td>No</td>\n",
       "      <td>never</td>\n",
       "      <td>never</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "      <td>Laptop/Desktop</td>\n",
       "      <td>cough/bad_cough_2021-08-24T08:08:28.798Z</td>\n",
       "      <td>1</td>\n",
       "      <td>17.322667</td>\n",
       "      <td>17.749333</td>\n",
       "      <td>18.517333</td>\n",
       "      <td>cough/bad_cough_2021-08-24T08:08:28.798Z.wav</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1309</th>\n",
       "      <td>Female</td>\n",
       "      <td>28</td>\n",
       "      <td>Ho Chi Minh</td>\n",
       "      <td>['No']</td>\n",
       "      <td>['No']</td>\n",
       "      <td>No</td>\n",
       "      <td>never</td>\n",
       "      <td>never</td>\n",
       "      <td>No</td>\n",
       "      <td>True</td>\n",
       "      <td>iPhone 8</td>\n",
       "      <td>cough/good_cough_2021-09-24T06:33:54.423Z</td>\n",
       "      <td>0</td>\n",
       "      <td>21.504000</td>\n",
       "      <td>26.624000</td>\n",
       "      <td>25.600000</td>\n",
       "      <td>cough/good_cough_2021-09-24T06:33:54.423Z.wav</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     sex_choice  age_choice current_city  \\\n",
       "1305     Female          20    thanh hóa   \n",
       "1306       Male          32  Ho Chi Minh   \n",
       "1307       Male          23  Ho Chi Minh   \n",
       "1308       Male          18  Ho Chi Minh   \n",
       "1309     Female          28  Ho Chi Minh   \n",
       "\n",
       "                                 symptoms_status_choice  \\\n",
       "1305                                             ['No']   \n",
       "1306                              ['fever', 'headache']   \n",
       "1307  ['fever', 'chills', 'sorethroat', 'drycough', ...   \n",
       "1308                         ['wetcough', 'sorethroat']   \n",
       "1309                                             ['No']   \n",
       "\n",
       "     medical_condition_choice insomnia_status_choice smoke_status_choice  \\\n",
       "1305                   ['No']                     No               never   \n",
       "1306                   ['No']                     No               never   \n",
       "1307                   ['No']                      1               1to10   \n",
       "1308                   ['No']                     No               never   \n",
       "1309                   ['No']                     No               never   \n",
       "\n",
       "     cov19_status_choice hospital_choice  cough_noise    device_model  \\\n",
       "1305               never              No         True    OPPO CPH1933   \n",
       "1306              last14              No         True  Laptop/Desktop   \n",
       "1307              last14              No         True  Laptop/Desktop   \n",
       "1308               never             NaN         True  Laptop/Desktop   \n",
       "1309               never              No         True        iPhone 8   \n",
       "\n",
       "                                     file_cough  label  cough_duration  \\\n",
       "1305  cough/good_cough_2021-08-15T13:43:33.132Z      0       25.856000   \n",
       "1306   cough/bad_cough_2021-09-16T07:18:48.594Z      1       29.350042   \n",
       "1307   cough/bad_cough_2021-09-06T11:04:49.842Z      1       18.432000   \n",
       "1308   cough/bad_cough_2021-08-24T08:08:28.798Z      1       17.322667   \n",
       "1309  cough/good_cough_2021-09-24T06:33:54.423Z      0       21.504000   \n",
       "\n",
       "      nose_duration  mouth_duration  \\\n",
       "1305      17.664000       22.357333   \n",
       "1306      29.814438       29.814438   \n",
       "1307      17.152000       15.530667   \n",
       "1308      17.749333       18.517333   \n",
       "1309      26.624000       25.600000   \n",
       "\n",
       "                                          file_path  label_symptom  \\\n",
       "1305  cough/good_cough_2021-08-15T13:43:33.132Z.wav              0   \n",
       "1306   cough/bad_cough_2021-09-16T07:18:48.594Z.wav              1   \n",
       "1307   cough/bad_cough_2021-09-06T11:04:49.842Z.wav              1   \n",
       "1308   cough/bad_cough_2021-08-24T08:08:28.798Z.wav              1   \n",
       "1309  cough/good_cough_2021-09-24T06:33:54.423Z.wav              0   \n",
       "\n",
       "      label_abnormal  label_covid  \n",
       "1305               0            0  \n",
       "1306               1            1  \n",
       "1307               1            1  \n",
       "1308               1            0  \n",
       "1309               0            0  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Classifier = 'SVM'\n",
    "PRETRAIN = 'FRILL'  # TRILL, FRILL, OpenSmile, OpenSmileBoAW, DeepSpectrum\n",
    "dataset_type = 'SoundDr'\n",
    "\n",
    "#sample rate\n",
    "SR = 44100\n",
    "#100 ms\n",
    "FRAME_LEN = int(SR / 10)\n",
    "#50% overlap, meaning 5ms hop length\n",
    "HOP = int(FRAME_LEN / 2)\n",
    "#the MFCC dimension\n",
    "MFCC_dim = 13\n",
    "codebook_size = 1000\n",
    "fold_num = 5\n",
    "seed = 2022\n",
    "\n",
    "if dataset_type == 'SoundDr':\n",
    "    PERIOD = 15\n",
    "    SR = 48000\n",
    "    DIR_DATA = \"./sounddr_data/\"\n",
    "\n",
    "    df = pd.read_csv(DIR_DATA + 'data.csv')\n",
    "    df['file_path'] = df['file_cough'] + '.wav'\n",
    "\n",
    "    df['label_symptom'] = (df['symptoms_status_choice'].map(str) != \"['No']\").astype(int)\n",
    "    df['label_abnormal'] = ((df['symptoms_status_choice'].map(str) != \"['No']\") | (df['cov19_status_choice'] != 'never')).astype(int)\n",
    "    df['label_covid'] = (df['cov19_status_choice'] != 'never').astype(int)\n",
    "elif dataset_type == 'CoughVid':\n",
    "    PERIOD = 10\n",
    "    SR = 22050\n",
    "\n",
    "    DIR_DATA = './coughvid_data/'\n",
    "\n",
    "    VidData   = pd.read_csv(os.path.join(DIR_DATA, 'public_dataset/metadata_compiled.csv'), header=0)\n",
    "    VidData   = VidData.loc[VidData['cough_detected'] >= 0.9][['uuid','fever_muscle_pain','respiratory_condition','status', 'quality_1', 'age', 'gender']]\n",
    "    VidData.dropna(subset=['uuid','fever_muscle_pain','respiratory_condition','status'], inplace=True)\n",
    "    VidData = VidData[(VidData['quality_1'] != 'no_cough') & (VidData['quality_1'] != 'poor')]\n",
    "    VidData = VidData[(VidData['status'] != 'symptomatic') & (VidData['status'].notna())]\n",
    "    VidData['label_covid'] = (VidData['status'] == 'COVID-19').astype(int)\n",
    "\n",
    "    extradata = VidData.loc[VidData['status']=='COVID-19']\n",
    "    notradata = VidData.loc[VidData['status']!='COVID-19']\n",
    "\n",
    "    df = pd.concat([extradata, notradata], ignore_index= True)\n",
    "    df['file_path'] = df['uuid'].apply(lambda x: 'public_dataset/' + x + '.webm')\n",
    "    def g(x):\n",
    "        for i in x:\n",
    "            if i is True:\n",
    "                return 1\n",
    "        return 0\n",
    "    df['label_abnormal'] = df[['fever_muscle_pain', 'respiratory_condition', 'label_covid']].apply(lambda x: g(x), axis=1)\n",
    "else:\n",
    "    PERIOD = 5\n",
    "    SR = 44100\n",
    "    DIR_DATA = \"./coswara_data/\"\n",
    "\n",
    "    join_by = pd.read_csv(os.path.join(DIR_DATA, 'combined_data.csv'))\n",
    "    df_list = []\n",
    "    for each in os.listdir(DIR_DATA):\n",
    "        for path in tqdm(glob(DIR_DATA + each + '/*/cough-shallow.wav')):\n",
    "            temp = pd.DataFrame(columns=['id', 'DIR'])\n",
    "            temp['id'] = [path.split('/')[-2]]\n",
    "            temp['DIR'] = [path]\n",
    "            temp = pd.merge(left=temp,right=join_by,on='id',how='inner')\n",
    "\n",
    "            temp['label_cough'] = (temp['cough'] == True).astype(int)\n",
    "\n",
    "            temp['file_path'] = each + '/' + temp['id'] + '/cough-shallow.wav'\n",
    "            temp['label_covid'] = temp['covid_status'].apply(lambda x: 1 if x == 'positive_mild' or x =='positive_moderate' or x == 'COVID-19' else 0)\n",
    "            df_list.append(temp)\n",
    "    df = pd.concat(df_list)\n",
    "    def g(x):\n",
    "        for i in x:\n",
    "            if i is True:\n",
    "                return 1\n",
    "        return 0\n",
    "    df['label_abnormal'] = df[['st', 'bd', 'cld', 'pneumonia', 'others_resp', 'asthma', 'label_covid']].apply(lambda x: g(x), axis=1)\n",
    "\n",
    "target_col = 'label_abnormal'\n",
    "OUTPUT_DIR = DIR_DATA + 'output/'\n",
    "os.makedirs(OUTPUT_DIR, exist_ok = True)\n",
    "\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def get_duration(filename, mono=True, res_type=\"kaiser_fast\"):\n",
    "    duration = 0\n",
    "    sr = SR\n",
    "    try:\n",
    "        y, sr = librosa.load(filename, sr=None, mono=mono, res_type=res_type)\n",
    "        duration = librosa.get_duration(y=y, sr=sr)\n",
    "    except:\n",
    "        print('Error file:' + filename)\n",
    "    return duration, sr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    825\n",
       "1    485\n",
       "Name: label_abnormal, dtype: int64"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[target_col].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def crop_or_pad(y, length):\n",
    "    \"\"\"\n",
    "    Crop audio random by length\n",
    "    \"\"\"\n",
    "    if len(y) < length:\n",
    "        y = np.concatenate([y, np.zeros(length-len(y))])\n",
    "    elif len(y) > length:\n",
    "        cut = random.randint(0, len(y) - length)\n",
    "        y = y[cut:cut+length]\n",
    "    return y\n",
    "\n",
    "def merge_feature(list_features):\n",
    "    \"\"\"\n",
    "      Merge numpy array features\n",
    "      Args:\n",
    "        - list_features: list of numpy array features\n",
    "                         :type: a list of numpy arrays \n",
    "      Returns:\n",
    "        - features: the concatenate numpy array along axis=1\n",
    "                    :type: a numpy array                 \n",
    "    \"\"\"      \n",
    "    features = np.concatenate(list_features, axis=1)\n",
    "    features = np.nan_to_num(features)\n",
    "    features = np.clip(features, -np.finfo(np.float32).max, np.finfo(np.float32).max)\n",
    "    return features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def compute_metrics(cfs_matrix):\n",
    "    \"\"\"\n",
    "      Calculate common metrics based on the confusion matrix\n",
    "      Args:\n",
    "        - cfs_matrix: a sklearn confusion matrix \n",
    "                      :type: a ndarray of shape (n_classes, n_classes)\n",
    "      Returns:\n",
    "        - precision: the precision of the prediction\n",
    "                     :type: float  \n",
    "        - recall: the recall of the prediction\n",
    "                  :type: float  \n",
    "        - f1: the f1-score of the prediction\n",
    "              :type: float                       \n",
    "    \"\"\"     \n",
    "    precision = cfs_matrix[1,1] / (cfs_matrix[1,1] + cfs_matrix[0,1])\n",
    "    recall = cfs_matrix[1,1] / (cfs_matrix[1,1] + cfs_matrix[1,0])\n",
    "    f1 = 2 * (precision * recall) / (precision + recall)\n",
    "    return precision, recall, f1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extract Audio Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions to extract feature follow by feature type in **Configuration** Cell(PRETRAIN) and save to cache file\n",
    "If exist file, we pass over it for save time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(OUTPUT_DIR + PRETRAIN + '.pickle'):\n",
    "    if PRETRAIN == 'TRILL':\n",
    "        import tensorflow.compat.v2 as tf\n",
    "        tf.enable_v2_behavior()\n",
    "        import tensorflow_hub as hub\n",
    "        import time\n",
    "        trill_model = hub.load('https://tfhub.dev/google/nonsemantic-speech-benchmark/trill/3')\n",
    "\n",
    "        def make_nonsemantic_trill_feat(filename):\n",
    "            try:\n",
    "                waveform, _ = librosa.load(os.path.join(DIR_DATA, filename), sr=16000, mono=True, res_type=\"kaiser_fast\")\n",
    "                if 2048 > waveform.shape[-1]:\n",
    "                    print('File length < 2048')\n",
    "                    return None, filename\n",
    "                embeddings = trill_model(samples=waveform, sample_rate=16000)['embedding']\n",
    "                mean_emb = embeddings.numpy().mean(axis=0)\n",
    "                std_emb = embeddings.numpy().std(axis=0)\n",
    "            except Exception as e:\n",
    "                print('Error: ' + str(e))\n",
    "                return None, filename\n",
    "            return np.concatenate((mean_emb, std_emb)), filename\n",
    "    elif PRETRAIN == 'FRILL':\n",
    "        import tensorflow.compat.v2 as tf\n",
    "        tf.enable_v2_behavior()\n",
    "        import tensorflow_hub as hub\n",
    "\n",
    "        frill_nofrontend_model = hub.load('https://tfhub.dev/google/nonsemantic-speech-benchmark/frill-nofrontend/1')\n",
    "\n",
    "        def stabilized_log(data, additive_offset, floor):\n",
    "          \"\"\"TF version of mfcc_mel.StabilizedLog.\"\"\"\n",
    "          return tf.math.log(tf.math.maximum(data, floor) + additive_offset)\n",
    "\n",
    "\n",
    "        def log_mel_spectrogram(data,\n",
    "                                audio_sample_rate,\n",
    "                                num_mel_bins=64,\n",
    "                                log_additive_offset=0.001,\n",
    "                                log_floor=1e-12,\n",
    "                                window_length_secs=0.025,\n",
    "                                hop_length_secs=0.010,\n",
    "                                fft_length=None):\n",
    "            \"\"\"TF version of mfcc_mel.LogMelSpectrogram.\"\"\"\n",
    "            window_length_samples = int(round(audio_sample_rate * window_length_secs))\n",
    "            hop_length_samples = int(round(audio_sample_rate * hop_length_secs))\n",
    "            if not fft_length:\n",
    "                fft_length = 2 ** int(np.ceil(np.log(window_length_samples) / np.log(2.0)))\n",
    "\n",
    "            spectrogram = tf.abs(\n",
    "                tf.signal.stft(\n",
    "                    tf.cast(data, tf.dtypes.float64),\n",
    "                    frame_length=window_length_samples,\n",
    "                    frame_step=hop_length_samples,\n",
    "                    fft_length=fft_length,\n",
    "                    window_fn=tf.signal.hann_window,\n",
    "                )\n",
    "            )\n",
    "\n",
    "            to_mel = tf.signal.linear_to_mel_weight_matrix(\n",
    "                num_mel_bins=num_mel_bins,\n",
    "                num_spectrogram_bins=fft_length // 2 + 1,\n",
    "                sample_rate=audio_sample_rate,\n",
    "                lower_edge_hertz=125.0,\n",
    "                upper_edge_hertz=7500.0,\n",
    "                dtype=tf.dtypes.float64\n",
    "            )\n",
    "\n",
    "            mel = spectrogram @ to_mel\n",
    "            log_mel = stabilized_log(mel, log_additive_offset, log_floor)\n",
    "            return log_mel\n",
    "\n",
    "        def compute_frontend_features(samples, sr, frame_hop, n_required=16000, num_mel_bins=64, frame_width=96):\n",
    "            if samples.dtype == np.int16:\n",
    "                samples = tf.cast(samples, np.float32) / np.iinfo(np.int16).max\n",
    "            if samples.dtype == np.float64:\n",
    "                samples = tf.cast(samples, np.float32)\n",
    "            assert samples.dtype == np.float32, samples.dtype\n",
    "            n = tf.size(samples)\n",
    "            samples = tf.cond(\n",
    "                n < n_required,\n",
    "                lambda: tf.pad(samples, [(0, n_required - n)]),\n",
    "                lambda: samples\n",
    "            )\n",
    "            mel = log_mel_spectrogram(samples, sr, num_mel_bins=num_mel_bins)\n",
    "            mel = tf.signal.frame(mel, frame_length=frame_width, frame_step=frame_hop, axis=0)\n",
    "            return mel\n",
    "\n",
    "        def make_nonsemantic_frill_nofrontend_feat(filename):\n",
    "            try:\n",
    "                waveform, _ = librosa.load(os.path.join(DIR_DATA, filename), sr=16000, mono=True, res_type=\"kaiser_fast\")\n",
    "                if 2048 > waveform.shape[-1]:\n",
    "                    print('File length < 2048')\n",
    "                    return None, filename\n",
    "                frontend_feats = tf.expand_dims(compute_frontend_features(waveform, 16000, frame_hop=17), axis=-1).numpy().astype(np.float32)\n",
    "                assert frontend_feats.shape[1:] == (96, 64, 1)\n",
    "\n",
    "                embeddings = frill_nofrontend_model(frontend_feats)['embedding']\n",
    "                mean_emb = embeddings.numpy().mean(axis=0)\n",
    "                std_emb = embeddings.numpy().std(axis=0)\n",
    "            except Exception as e:\n",
    "                print('Error: ' + str(e))\n",
    "                return None, filename\n",
    "            return np.concatenate((mean_emb, std_emb)), filename\n",
    "    elif PRETRAIN == 'OpenSmile':\n",
    "        SMILE = opensmile.Smile(feature_set=opensmile.FeatureSet.ComParE_2016, feature_level=opensmile.FeatureLevel.Functionals)\n",
    "\n",
    "        def make_opensmile_feat(filename):\n",
    "            try:\n",
    "                if dataset_type == 'CoughVid':\n",
    "                    # convert audio to wav file\n",
    "                    # for i in public_dataset/*.webm; do ffmpeg -i \"$i\" -c:a pcm_f32le \"wav_dataset/${i%.*}.wav\"; done\n",
    "                    path = os.path.join(DIR_DATA, 'wav_dataset', filename[:-5]+'.wav')\n",
    "                else:\n",
    "                    path = os.path.join(DIR_DATA, filename)\n",
    "                waveform, _ = librosa.load(path, sr=SR, mono=True, res_type=\"kaiser_fast\")\n",
    "                if 2048 > waveform.shape[-1]:\n",
    "                    print('File length < 2048')\n",
    "                feature = SMILE.process_signal(waveform, SR).values.reshape(-1)  # process_file, process_signal\n",
    "            except Exception as e:\n",
    "                print('Error: ' + str(e))\n",
    "                return None, filename\n",
    "            return feature, filename\n",
    "    elif PRETRAIN == 'OpenSmileBoAW':\n",
    "        lld = opensmile.Smile(feature_set=opensmile.FeatureSet.ComParE_2016, feature_level='lld')\n",
    "        lld_deltas = opensmile.Smile(feature_set=opensmile.FeatureSet.ComParE_2016, feature_level='lld_de')\n",
    "\n",
    "        def make_opensmileboaw_feat(filename):\n",
    "            try:\n",
    "                if dataset_type == 'CoughVid':\n",
    "                    # convert audio to wav file\n",
    "                    # for i in public_dataset/*.webm; do ffmpeg -i \"$i\" -c:a pcm_f32le \"wav_dataset/${i%.*}.wav\"; done\n",
    "                    path = os.path.join(DIR_DATA, 'wav_dataset', filename[:-5]+'.wav')\n",
    "                else:\n",
    "                    path = os.path.join(DIR_DATA, filename)\n",
    "                waveform, _ = librosa.load(path, sr=SR, mono=True, res_type=\"kaiser_fast\")\n",
    "                if 2048 > waveform.shape[-1]:\n",
    "                    print('File length < 2048')\n",
    "                # shape of feature 2991*65 + 2993*65 = 194610\n",
    "                waveform = crop_or_pad(waveform, SR*PERIOD)\n",
    "                feature = np.concatenate([\n",
    "                    lld.process_signal(waveform, SR).values.reshape(-1),\n",
    "                    lld_deltas.process_signal(waveform, SR).values.reshape(-1)\n",
    "                ], axis=0)\n",
    "            except Exception as e:\n",
    "                print('Error: ' + str(e))\n",
    "                return None, filename\n",
    "            return feature, filename\n",
    "    elif PRETRAIN == 'DeepSpectrum':\n",
    "        import tensorflow as tf\n",
    "        \n",
    "        base_model = tf.keras.applications.densenet.DenseNet121(weights='imagenet')\n",
    "        layer = base_model.layers[-2].name\n",
    "        outputs = (base_model.get_layer(layer) if not hasattr(base_model.get_layer(layer), \"output\") else base_model.get_layer(layer).output)\n",
    "        deep_model = tf.keras.models.Model(inputs=base_model.input, outputs=outputs)\n",
    "\n",
    "        def make_deepspect_feat(filename):\n",
    "            try:\n",
    "                waveform, _ = librosa.load(os.path.join(DIR_DATA, filename), sr=SR, mono=True, res_type=\"kaiser_fast\")\n",
    "                if 2048 > waveform.shape[-1]:\n",
    "                    print('File length < 2048')\n",
    "                    return None, filename\n",
    "                waveform = crop_or_pad(waveform, SR*PERIOD)\n",
    "                melspec = librosa.feature.melspectrogram(waveform, sr=SR, n_mels=128, n_fft=int(0.032*SR), hop_length=int(0.016*SR))\n",
    "                melspec = cv2.resize(melspec, (224, 224))\n",
    "                image = np.stack([melspec, melspec, melspec], axis=-1)\n",
    "                image = tf.keras.applications.densenet.preprocess_input(image)\n",
    "                feature = deep_model.predict(np.array([image]))[0]\n",
    "                return feature, filename\n",
    "            except Exception as e:\n",
    "                print('Error: ' + str(e))\n",
    "                return None, filename\n",
    "            return feature, filename"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Extract Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def get_features_of_list_audio(df):\n",
    "    X_features = []\n",
    "    df['error'] = 0\n",
    "    for idx, r in tqdm(df.iterrows(), total=len(df)):\n",
    "        if PRETRAIN == 'TRILL':\n",
    "            feature, filename = make_nonsemantic_trill_feat(r['file_path'])\n",
    "        elif PRETRAIN == 'OpenSmileBoAW':\n",
    "            feature, filename = make_opensmileboaw_feat(r['file_path'])\n",
    "        elif PRETRAIN == 'OpenSmile':\n",
    "            feature, filename = make_opensmile_feat(r['file_path'])\n",
    "        elif PRETRAIN == 'DeepSpectrum':\n",
    "            feature, filename = make_deepspect_feat(r['file_path'])\n",
    "        else:\n",
    "            feature, filename = make_nonsemantic_frill_nofrontend_feat(r['file_path'])\n",
    "        \n",
    "        if feature is None:\n",
    "            df['error'][df['file_path'] == filename] = 1\n",
    "        else:\n",
    "            X_features.append(feature)\n",
    "    return np.array(X_features), df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data feature shape: ((1310, 4096), 1310)\n"
     ]
    }
   ],
   "source": [
    "if not os.path.exists(OUTPUT_DIR + PRETRAIN + '.pickle'):\n",
    "    X_features, df = get_features_of_list_audio(df)\n",
    "    df.to_csv(os.path.join(OUTPUT_DIR, 'data.csv'), index=False)\n",
    "    pickle.dump({\n",
    "        'X_trill_features': X_features\n",
    "    }, open(OUTPUT_DIR + PRETRAIN + '.pickle', \"wb\" ))\n",
    "else:\n",
    "    df = pd.read_csv(os.path.join(OUTPUT_DIR, 'data.csv'))\n",
    "    f = pickle.load(open(OUTPUT_DIR + PRETRAIN + '.pickle', \"rb\" ))\n",
    "    X_features = f['X_trill_features']\n",
    "df = df[df['error'] == 0].reset_index(drop=True)\n",
    "\n",
    "# if PRETRAIN == 'OpenSmile':\n",
    "scaler = StandardScaler()\n",
    "X = scaler.fit_transform(merge_feature([X_features]))\n",
    "# else:\n",
    "# X = merge_feature([X_features])\n",
    "print(f\"Data feature shape: {X.shape, len(df)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# calculate TF\n",
    "def compute_TF(word_dict, bow):\n",
    "    tf_dict = {}\n",
    "    bow_count = len(bow)\n",
    "    for word, count in word_dict.iteritems():\n",
    "        tf_dict[word] = count/float(bow_count)\n",
    "\n",
    "    return tf_dict\n",
    "\n",
    "def compute_IDF(doc_list):\n",
    "    import math\n",
    "    idf_dict = {}\n",
    "    N = len(doc_list)\n",
    "\n",
    "    #count number of documents that contain this word\n",
    "    idf_dict = dict.fromkeys(doc_list[0].keys(), 0)\n",
    "    for doc in doc_list:\n",
    "        for word, count in doc.items():\n",
    "            if count > 0:\n",
    "                idf_dict[word] += 1\n",
    "\n",
    "    for word, count in idf_dict.items():\n",
    "        idf_dict[word] = math.log(N/float(count))\n",
    "\n",
    "    return idf_dict\n",
    "\n",
    "def compute_TFIDF(tf_bow, idfs):\n",
    "    tfidf = {}\n",
    "    for word, val in tf_bow.items():\n",
    "        tfidf[word] = val*idfs[word]\n",
    "    return tfidf\n",
    "\n",
    "# getting bag of words as sample\n",
    "def bow(feature, kmeans):\n",
    "    unique, counts = np.unique(kmeans.predict(feature), return_counts=True)\n",
    "    no_of_centers = kmeans.cluster_centers_.shape[0]\n",
    "    sample = [0]*no_of_centers\n",
    "    for j in range(no_of_centers):\n",
    "        if j in unique:\n",
    "            sample[j] = counts[np.where(unique == j)][0]\n",
    "    return sample\n",
    "\n",
    "def boaw(feature, codebook, num_assign=10):\n",
    "    # Faster for large dataset\n",
    "    x_idx = cKDTree(codebook).query(feature, k=num_assign)[1].reshape(-1)\n",
    "\n",
    "    word_dict = dict.fromkeys([j for j in range(len(codebook))], 0)\n",
    "    for word in x_idx:\n",
    "        word_dict[word] += 1\n",
    "\n",
    "    boaw = np.array([*word_dict.values()])\n",
    "    return boaw, word_dict\n",
    "\n",
    "def generate_boaw(X, codebook):\n",
    "    boaws = []\n",
    "    pool = multiprocessing.Pool(multiprocessing.cpu_count())\n",
    "    with tqdm(total=len(X)) as t:\n",
    "        for x_boaw, word_dict in pool.imap(partial(boaw, codebook=codebook), X):\n",
    "            tf_dict = {}\n",
    "            for word, count in word_dict.items():\n",
    "                tf_dict[word] = count\n",
    "            boaws.append(x_boaw*np.log(np.array([*tf_dict.values()])))\n",
    "            t.update(1)\n",
    "    pool.close()\n",
    "    pool.join()\n",
    "    return boaws\n",
    "\n",
    "def evaluate(ensem_preds, targets):\n",
    "    \"\"\"\n",
    "      Evaluate the prediction by providing metrics & also the best threshold (to get the highest f1-score)\n",
    "      Ex: AUC, Accurary, Precision, Recall, F1-Score.\n",
    "      Then print these metrics\n",
    "      Args:\n",
    "        - ensem_preds: predictions for ids \n",
    "                       :type: a numpy array\n",
    "        - targets: the actual results of ids \n",
    "                   :type: a numpy array                 \n",
    "      Returns:\n",
    "        - None                  \n",
    "    \"\"\"     \n",
    "    best_th = 0\n",
    "    best_score = 0\n",
    "\n",
    "    for th in np.arange(0.0, 0.6, 0.01):\n",
    "        pred = (ensem_preds > th).astype(int)\n",
    "        score = f1_score(targets, pred)\n",
    "        if score > best_score:\n",
    "            best_th = th\n",
    "            best_score = score\n",
    "\n",
    "    print(f\"\\nAUC score: {roc_auc_score(targets, ensem_preds):12.4f}\")\n",
    "    print(f\"Best threshold {best_th:12.4f}\")\n",
    "\n",
    "    preds = (ensem_preds > best_th).astype(int)\n",
    "    # print(classification_report(targets, preds, digits=3))\n",
    "\n",
    "    cm1 = confusion_matrix(targets, preds)\n",
    "    print('\\nConfusion Matrix : \\n', cm1)\n",
    "    precision, recall, f1 = compute_metrics(cm1)\n",
    "    \n",
    "    print('\\n=============')\n",
    "    print (f'Precision    : {precision:12.4f}')\n",
    "    \n",
    "    print(f'Recall : {recall:12.4f}')\n",
    "    \n",
    "    print(f'F1 Score : {f1:12.4f}')\n",
    "    \n",
    "    total1=sum(sum(cm1))\n",
    "\n",
    "    print('\\n=============')\n",
    "    accuracy1=(cm1[0,0]+cm1[1,1])/total1\n",
    "    print (f'Accuracy    : {accuracy1:12.4f}')\n",
    "\n",
    "def get_model(pos_scale, c=1):\n",
    "    if Classifier == 'SVM':\n",
    "        model = LinearSVC(C=c, class_weight='balanced', random_state=seed)\n",
    "    else:\n",
    "        model = xgb.XGBClassifier(max_depth=6, learning_rate=0.07,\n",
    "            scale_pos_weight=pos_scale,\n",
    "            n_estimators=200,\n",
    "            subsample=1,\n",
    "            colsample_bytree=1,\n",
    "            eta=1, objective='binary:logistic',\n",
    "            eval_metric='auc'\n",
    "        )\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Train COVID-19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    " if not os.path.exists(OUTPUT_DIR + 'df_label_covid_5fold.csv'):\n",
    "    folds = df.copy()\n",
    "    Fold = StratifiedKFold(n_splits=fold_num, shuffle=True, random_state=seed)\n",
    "    for n, (train_index, val_index) in enumerate(Fold.split(folds, folds['label_covid'])):\n",
    "        folds.loc[val_index, 'fold'] = int(n)\n",
    "    folds['fold'] = folds['fold'].astype(int)\n",
    "    folds.to_csv(OUTPUT_DIR + 'df_label_covid_5fold.csv', index=False)\n",
    "else:\n",
    "    folds = pd.read_csv(OUTPUT_DIR + 'df_label_covid_5fold.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.786127167630058\n",
      "0.7471726190476189\n",
      "0.8586393331831494\n",
      "0.8021701584440941\n",
      "0.8539460839528423\n",
      "0.8228955470451302\n",
      "(!) cv5 AUC  0.816964748334567 0.04056941158112057\n",
      "\n",
      "AUC score:       0.8168\n",
      "Best threshold       0.0000\n",
      "\n",
      "Confusion Matrix : \n",
      " [[803 161]\n",
      " [ 69 277]]\n",
      "\n",
      "=============\n",
      "Precision    :       0.6324\n",
      "Recall :       0.8006\n",
      "F1 Score :       0.7066\n",
      "\n",
      "=============\n",
      "Accuracy    :       0.8244\n"
     ]
    }
   ],
   "source": [
    "y = folds['label_covid']\n",
    "pos_scale = (y == 0).sum() / (y == 1).sum()\n",
    "print(pos_scale)\n",
    "targets = []\n",
    "preds = []\n",
    "aucs = []\n",
    "\n",
    "for fold in range(5):\n",
    "    train_idx = folds['fold'] != fold\n",
    "    valid_idx = folds['fold'] == fold\n",
    "    X_train = X[train_idx]\n",
    "    y_train = y[train_idx]\n",
    "    X_val = X[valid_idx]\n",
    "    y_val = y[valid_idx]\n",
    "\n",
    "    targets.append(y_val)\n",
    "    \n",
    "    if PRETRAIN == 'OpenSmileBoAW':\n",
    "        from sklearn.cluster import MiniBatchKMeans\n",
    "\n",
    "        p = int(X_train.shape[1]/130)\n",
    "        X_train1 = X_train[:, :(p-1)*65].reshape(-1, p-1, 65)\n",
    "        X_train2 = X_train[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "        X_val1 = X_val[:, :(p-1)*65].reshape(-1, (p-1), 65)\n",
    "        X_val2 = X_val[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "\n",
    "        t = time.time()\n",
    "        vocab1 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train1.reshape(-1, 65))\n",
    "        codebook1 = vocab1.cluster_centers_\n",
    "        print(\"Kmean =\", time.time()-t)\n",
    "        X_train1 = generate_boaw(X_train1, codebook1)\n",
    "        X_val1 = generate_boaw(X_val1, codebook1)\n",
    "        \n",
    "        vocab2 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train2.reshape(-1, 65))\n",
    "        codebook2 = vocab2.cluster_centers_\n",
    "        X_train2 = generate_boaw(X_train2, codebook2)\n",
    "        X_val2 = generate_boaw(X_val2, codebook2)\n",
    "\n",
    "        X_train = np.nan_to_num(np.concatenate([X_train1, X_train2], axis=1))  # [num_sample, codebook_size*2]\n",
    "        X_val = np.nan_to_num(np.concatenate([X_val1, X_val2], axis=1))\n",
    "\n",
    "    model = get_model(pos_scale=pos_scale)\n",
    "    model.fit(X_train, y_train)\n",
    "\n",
    "    pred = model.predict(X_val)\n",
    "    preds.append(pred)\n",
    "    auc = roc_auc_score(y_val, pred)\n",
    "    aucs.append(auc)\n",
    "    del model\n",
    "\n",
    "targets = np.concatenate(targets)\n",
    "preds = np.concatenate(preds)\n",
    "\n",
    "print(\"(!) cv5 AUC \", np.mean(aucs), np.std(aucs))\n",
    "evaluate(preds, targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Train abnormal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    " if not os.path.exists(OUTPUT_DIR + 'df_label_abnormal_5fold.csv'):\n",
    "    folds = df.copy()\n",
    "    Fold = StratifiedKFold(n_splits=fold_num, shuffle=True, random_state=seed)\n",
    "    for n, (train_index, val_index) in enumerate(Fold.split(folds, folds['label_abnormal'])):\n",
    "        folds.loc[val_index, 'fold'] = int(n)\n",
    "    folds['fold'] = folds['fold'].astype(int)\n",
    "    folds.to_csv(OUTPUT_DIR + 'df_label_abnormal_5fold.csv', index=False)\n",
    "else:\n",
    "    folds = pd.read_csv(OUTPUT_DIR + 'df_label_abnormal_5fold.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.7010309278350515\n",
      "0.7823180256169946\n",
      "0.7259293970634177\n",
      "0.7965635738831615\n",
      "0.7219931271477663\n",
      "0.723805060918463\n",
      "(!) cv5 AUC  0.7501218369259606 0.032442258421842146\n",
      "\n",
      "AUC score:       0.7501\n",
      "Best threshold       0.0000\n",
      "\n",
      "Confusion Matrix : \n",
      " [[741  84]\n",
      " [193 292]]\n",
      "\n",
      "=============\n",
      "Precision    :       0.7766\n",
      "Recall :       0.6021\n",
      "F1 Score :       0.6783\n",
      "\n",
      "=============\n",
      "Accuracy    :       0.7885\n"
     ]
    }
   ],
   "source": [
    "y = folds['label_abnormal']\n",
    "pos_scale = (y == 0).sum() / (y == 1).sum()\n",
    "print(pos_scale)\n",
    "targets = []\n",
    "preds = []\n",
    "aucs = []\n",
    "\n",
    "def get_model(pos_scale):\n",
    "    model = xgb.XGBClassifier(\n",
    "        max_depth=7,\n",
    "        scale_pos_weight=pos_scale,\n",
    "        learning_rate=0.3,\n",
    "        n_estimators=200,\n",
    "        subsample=1,\n",
    "        colsample_bytree=1,\n",
    "        nthread=-1,\n",
    "        seed=seed,\n",
    "        eval_metric='logloss'\n",
    "    )\n",
    "    return model\n",
    "\n",
    "for fold in range(5):\n",
    "    train_idx = folds['fold'] != fold\n",
    "    valid_idx = folds['fold'] == fold\n",
    "    X_train = X[train_idx]\n",
    "    y_train = y[train_idx]\n",
    "    X_val = X[valid_idx]\n",
    "    y_val = y[valid_idx]\n",
    "\n",
    "    targets.append(y_val)\n",
    "    \n",
    "    if PRETRAIN == 'OpenSmileBoAW':\n",
    "        from sklearn.cluster import MiniBatchKMeans\n",
    "\n",
    "        p = int(X_train.shape[1]/130)\n",
    "        X_train1 = X_train[:, :(p-1)*65].reshape(-1, p-1, 65)\n",
    "        X_train2 = X_train[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "        X_val1 = X_val[:, :(p-1)*65].reshape(-1, (p-1), 65)\n",
    "        X_val2 = X_val[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "\n",
    "        t = time.time()\n",
    "        vocab1 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train1.reshape(-1, 65))\n",
    "        codebook1 = vocab1.cluster_centers_\n",
    "        print(\"Kmean =\", time.time()-t)\n",
    "        X_train1 = generate_boaw(X_train1, codebook1)\n",
    "        X_val1 = generate_boaw(X_val1, codebook1)\n",
    "        \n",
    "        vocab2 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train2.reshape(-1, 65))\n",
    "        codebook2 = vocab2.cluster_centers_\n",
    "        X_train2 = generate_boaw(X_train2, codebook2)\n",
    "        X_val2 = generate_boaw(X_val2, codebook2)\n",
    "\n",
    "        X_train = np.nan_to_num(np.concatenate([X_train1, X_train2], axis=1))  # [num_sample, codebook_size*2]\n",
    "        X_val = np.nan_to_num(np.concatenate([X_val1, X_val2], axis=1))\n",
    "\n",
    "    model = get_model(pos_scale=pos_scale)\n",
    "\n",
    "    pred = model.predict(X_val)\n",
    "    preds.append(pred)\n",
    "    auc = roc_auc_score(y_val, pred)\n",
    "    aucs.append(auc)\n",
    "    print(auc)\n",
    "    del model\n",
    "\n",
    "targets = np.concatenate(targets)\n",
    "preds = np.concatenate(preds)\n",
    "\n",
    "print(\"(!) cv5 AUC \", np.mean(aucs), np.std(aucs))\n",
    "evaluate(preds, targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test other model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IsolationForest\n",
    "Isolation Forest use Sklearn library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.7010309278350515\n"
     ]
    }
   ],
   "source": [
    "y = folds['label_abnormal']\n",
    "pos_scale = (y == 0).sum() / (y == 1).sum()\n",
    "print(pos_scale)\n",
    "\n",
    "def get_model(pos_scale):\n",
    "    model = IsolationForest(n_estimators=500, max_samples='auto', contamination=0.1, n_jobs=-1, random_state=seed) \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4850359262730396\n",
      "0.544704779756326\n",
      "0.5323961262105591\n",
      "0.5684473601999376\n",
      "0.5526398000624805\n",
      "(!) cv5 AUC  0.5366447985004685 0.02833278388113743\n",
      "\n",
      "AUC score:       0.5352\n",
      "Best threshold       0.0000\n",
      "\n",
      "Confusion Matrix : \n",
      " [[740  85]\n",
      " [437  48]]\n",
      "\n",
      "=============\n",
      "Precision    :       0.3609\n",
      "Recall :       0.0990\n",
      "F1 Score :       0.1553\n",
      "\n",
      "=============\n",
      "Accuracy    :       0.6015\n"
     ]
    }
   ],
   "source": [
    "targets = []\n",
    "preds = []\n",
    "aucs = []\n",
    "\n",
    "for fold in range(5):\n",
    "    train_idx = folds['fold'] != fold\n",
    "    valid_idx = folds['fold'] == fold\n",
    "    X_train = X[train_idx]\n",
    "    y_train = y[train_idx]\n",
    "    X_val = X[valid_idx]\n",
    "    y_val = y[valid_idx]\n",
    "\n",
    "    targets.append(y_val)\n",
    "    \n",
    "    if PRETRAIN == 'OpenSmileBoAW':\n",
    "        from sklearn.cluster import MiniBatchKMeans\n",
    "\n",
    "        p = int(X_train.shape[1]/130)\n",
    "        X_train1 = X_train[:, :(p-1)*65].reshape(-1, p-1, 65)\n",
    "        X_train2 = X_train[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "        X_val1 = X_val[:, :(p-1)*65].reshape(-1, (p-1), 65)\n",
    "        X_val2 = X_val[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "\n",
    "        t = time.time()\n",
    "        vocab1 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train1.reshape(-1, 65))\n",
    "        codebook1 = vocab1.cluster_centers_\n",
    "        print(\"Kmean =\", time.time()-t)\n",
    "        X_train1 = generate_boaw(X_train1, codebook1)\n",
    "        X_val1 = generate_boaw(X_val1, codebook1)\n",
    "        \n",
    "        vocab2 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train2.reshape(-1, 65))\n",
    "        codebook2 = vocab2.cluster_centers_\n",
    "        X_train2 = generate_boaw(X_train2, codebook2)\n",
    "        X_val2 = generate_boaw(X_val2, codebook2)\n",
    "\n",
    "        X_train = np.nan_to_num(np.concatenate([X_train1, X_train2], axis=1))  # [num_sample, codebook_size*2]\n",
    "        X_val = np.nan_to_num(np.concatenate([X_val1, X_val2], axis=1))\n",
    "\n",
    "    model = get_model(pos_scale)\n",
    "    model.fit(X_train)\n",
    "\n",
    "    scores = (-1.0) * model.decision_function(X_val)\n",
    "    pred = scores.flatten()\n",
    "    preds.append(pred)\n",
    "    auc = roc_auc_score(y_val, pred)\n",
    "    aucs.append(auc)\n",
    "    print(auc)\n",
    "    del model\n",
    "\n",
    "targets = np.concatenate(targets)\n",
    "preds = np.concatenate(preds)\n",
    "\n",
    "print(\"(!) cv5 AUC \", np.mean(aucs), np.std(aucs))\n",
    "evaluate(preds, targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XGBOD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.7010309278350515\n"
     ]
    }
   ],
   "source": [
    "y = folds['label_abnormal']\n",
    "pos_scale = (y == 0).sum() / (y == 1).sum()\n",
    "print(pos_scale)\n",
    "\n",
    "from pyod.models.xgbod import XGBOD\n",
    "\n",
    "def get_model(pos_scale):\n",
    "    # Model candidate\n",
    "    model = XGBOD(max_depth=7,\n",
    "        scale_pos_weight=pos_scale,\n",
    "        learning_rate=0.3,\n",
    "        n_estimators=200,\n",
    "        subsample=1,\n",
    "        colsample_bytree=1,\n",
    "        nthread=-1,\n",
    "        seed=seed,\n",
    "        eval_metric='logloss')\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[01:00:05] WARNING: ../src/learner.cc:516: \n",
      "Parameters: { silent } might not be used.\n",
      "\n",
      "  This may not be accurate due to some parameters are only used in language bindings but\n",
      "  passed down to XGBoost core.  Or some parameters are not used but slip through this\n",
      "  verification. Please open an issue if you find above cases.\n",
      "\n",
      "\n",
      "0.8647922524211183\n",
      "[01:05:11] WARNING: ../src/learner.cc:516: \n",
      "Parameters: { silent } might not be used.\n",
      "\n",
      "  This may not be accurate due to some parameters are only used in language bindings but\n",
      "  passed down to XGBoost core.  Or some parameters are not used but slip through this\n",
      "  verification. Please open an issue if you find above cases.\n",
      "\n",
      "\n",
      "0.8091221493283349\n",
      "[01:10:14] WARNING: ../src/learner.cc:516: \n",
      "Parameters: { silent } might not be used.\n",
      "\n",
      "  This may not be accurate due to some parameters are only used in language bindings but\n",
      "  passed down to XGBoost core.  Or some parameters are not used but slip through this\n",
      "  verification. Please open an issue if you find above cases.\n",
      "\n",
      "\n",
      "0.8431740081224617\n",
      "[01:15:19] WARNING: ../src/learner.cc:516: \n",
      "Parameters: { silent } might not be used.\n",
      "\n",
      "  This may not be accurate due to some parameters are only used in language bindings but\n",
      "  passed down to XGBoost core.  Or some parameters are not used but slip through this\n",
      "  verification. Please open an issue if you find above cases.\n",
      "\n",
      "\n",
      "0.8338019368947204\n",
      "[01:20:21] WARNING: ../src/learner.cc:516: \n",
      "Parameters: { silent } might not be used.\n",
      "\n",
      "  This may not be accurate due to some parameters are only used in language bindings but\n",
      "  passed down to XGBoost core.  Or some parameters are not used but slip through this\n",
      "  verification. Please open an issue if you find above cases.\n",
      "\n",
      "\n",
      "0.8003748828491096\n",
      "(!) cv5 AUC  0.830253045923149 0.023288289316191306\n",
      "\n",
      "AUC score:       0.8297\n",
      "Best threshold       0.3400\n",
      "\n",
      "Confusion Matrix : \n",
      " [[712 113]\n",
      " [151 334]]\n",
      "\n",
      "=============\n",
      "Precision    :       0.7472\n",
      "Recall :       0.6887\n",
      "F1 Score :       0.7167\n",
      "\n",
      "=============\n",
      "Accuracy    :       0.7985\n"
     ]
    }
   ],
   "source": [
    "targets = []\n",
    "preds = []\n",
    "aucs = []\n",
    "\n",
    "for fold in range(5):\n",
    "    train_idx = folds['fold'] != fold\n",
    "    valid_idx = folds['fold'] == fold\n",
    "    X_train = X[train_idx]\n",
    "    y_train = y[train_idx]\n",
    "    X_val = X[valid_idx]\n",
    "    y_val = y[valid_idx]\n",
    "\n",
    "    targets.append(y_val)\n",
    "    \n",
    "    if PRETRAIN == 'OpenSmileBoAW':\n",
    "        from sklearn.cluster import MiniBatchKMeans\n",
    "\n",
    "        p = int(X_train.shape[1]/130)\n",
    "        X_train1 = X_train[:, :(p-1)*65].reshape(-1, p-1, 65)\n",
    "        X_train2 = X_train[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "        X_val1 = X_val[:, :(p-1)*65].reshape(-1, (p-1), 65)\n",
    "        X_val2 = X_val[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "\n",
    "        t = time.time()\n",
    "        vocab1 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train1.reshape(-1, 65))\n",
    "        codebook1 = vocab1.cluster_centers_\n",
    "        print(\"Kmean =\", time.time()-t)\n",
    "        X_train1 = generate_boaw(X_train1, codebook1)\n",
    "        X_val1 = generate_boaw(X_val1, codebook1)\n",
    "        \n",
    "        vocab2 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train2.reshape(-1, 65))\n",
    "        codebook2 = vocab2.cluster_centers_\n",
    "        X_train2 = generate_boaw(X_train2, codebook2)\n",
    "        X_val2 = generate_boaw(X_val2, codebook2)\n",
    "\n",
    "        X_train = np.nan_to_num(np.concatenate([X_train1, X_train2], axis=1))  # [num_sample, codebook_size*2]\n",
    "        X_val = np.nan_to_num(np.concatenate([X_val1, X_val2], axis=1))\n",
    "\n",
    "    model = get_model(pos_scale)\n",
    "    model.fit(X_train, y_train)\n",
    "\n",
    "    pred = model.decision_function(X_val)  # predict raw outlier scores on test\n",
    "    auc = roc_auc_score(y_val, pred)\n",
    "    preds.append(pred)\n",
    "    aucs.append(auc)\n",
    "    print(auc)\n",
    "    del model\n",
    "\n",
    "targets = np.concatenate(targets)\n",
    "preds = np.concatenate(preds)\n",
    "\n",
    "print(\"(!) cv5 AUC \", np.mean(aucs), np.std(aucs))\n",
    "evaluate(preds, targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IForest\n",
    "Isolation Forest use pyod library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.7010309278350515\n"
     ]
    }
   ],
   "source": [
    "y = folds['label_abnormal']\n",
    "pos_scale = (y == 0).sum() / (y == 1).sum()\n",
    "print(pos_scale)\n",
    "\n",
    "from pyod.models.iforest import IForest\n",
    "\n",
    "def get_model(pos_scale):\n",
    "    model = IForest(n_estimators=500, max_samples='auto', contamination=0.1, n_jobs=-1, random_state=seed) \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4850359262730396\n",
      "0.544704779756326\n",
      "0.5323961262105591\n",
      "0.5684473601999376\n",
      "0.5526398000624805\n",
      "(!) cv5 AUC  0.5366447985004685 0.02833278388113743\n",
      "\n",
      "AUC score:       0.5352\n",
      "Best threshold       0.0000\n",
      "\n",
      "Confusion Matrix : \n",
      " [[740  85]\n",
      " [437  48]]\n",
      "\n",
      "=============\n",
      "Precision    :       0.3609\n",
      "Recall :       0.0990\n",
      "F1 Score :       0.1553\n",
      "\n",
      "=============\n",
      "Accuracy    :       0.6015\n"
     ]
    }
   ],
   "source": [
    "targets = []\n",
    "preds = []\n",
    "aucs = []\n",
    "\n",
    "for fold in range(5):\n",
    "    train_idx = folds['fold'] != fold\n",
    "    valid_idx = folds['fold'] == fold\n",
    "    X_train = X[train_idx]\n",
    "    y_train = y[train_idx]\n",
    "    X_val = X[valid_idx]\n",
    "    y_val = y[valid_idx]\n",
    "\n",
    "    targets.append(y_val)\n",
    "    \n",
    "    if PRETRAIN == 'OpenSmileBoAW':\n",
    "        from sklearn.cluster import MiniBatchKMeans\n",
    "\n",
    "        p = int(X_train.shape[1]/130)\n",
    "        X_train1 = X_train[:, :(p-1)*65].reshape(-1, p-1, 65)\n",
    "        X_train2 = X_train[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "        X_val1 = X_val[:, :(p-1)*65].reshape(-1, (p-1), 65)\n",
    "        X_val2 = X_val[:, (p-1)*65:].reshape(-1, (p+1), 65)\n",
    "\n",
    "        t = time.time()\n",
    "        vocab1 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train1.reshape(-1, 65))\n",
    "        codebook1 = vocab1.cluster_centers_\n",
    "        print(\"Kmean =\", time.time()-t)\n",
    "        X_train1 = generate_boaw(X_train1, codebook1)\n",
    "        X_val1 = generate_boaw(X_val1, codebook1)\n",
    "        \n",
    "        vocab2 = MiniBatchKMeans(n_clusters=codebook_size).fit(X_train2.reshape(-1, 65))\n",
    "        codebook2 = vocab2.cluster_centers_\n",
    "        X_train2 = generate_boaw(X_train2, codebook2)\n",
    "        X_val2 = generate_boaw(X_val2, codebook2)\n",
    "\n",
    "        X_train = np.nan_to_num(np.concatenate([X_train1, X_train2], axis=1))  # [num_sample, codebook_size*2]\n",
    "        X_val = np.nan_to_num(np.concatenate([X_val1, X_val2], axis=1))\n",
    "\n",
    "    model = get_model(pos_scale)\n",
    "    model.fit(X_train)\n",
    "\n",
    "    pred = model.decision_function(X_val)  # predict raw outlier scores on test\n",
    "    auc = roc_auc_score(y_val, pred)\n",
    "    preds.append(pred)\n",
    "    aucs.append(auc)\n",
    "    print(auc)\n",
    "    del model\n",
    "\n",
    "targets = np.concatenate(targets)\n",
    "preds = np.concatenate(preds)\n",
    "\n",
    "print(\"(!) cv5 AUC \", np.mean(aucs), np.std(aucs))\n",
    "evaluate(preds, targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp \"SoundDr_cough.ipynb\" \"$OUTPUT_DIR/SoundDr_cough_frill_svm_2022.ipynb\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}