lostanlen_waspaa2019_fig3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "\n",
    "from librosa.display import specshow\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import tqdm\n",
    "\n",
    "data_dir = \"/beegfs/vl1019/waspaa2019_data\"\n",
    "ccb_dir = os.path.join(data_dir, \"ccb18\")\n",
    "csv_name = 'CCB_2009Feb18to22Apr17.csv'\n",
    "csv_path = os.path.join(ccb_dir, csv_name)\n",
    "\n",
    "hop_length = 2**7\n",
    "n_fft = 2 * hop_length\n",
    "sr = 2000\n",
    "half_clip_n_cols = 16\n",
    "low_freq_bin = 5\n",
    "\n",
    "pcen_version_id = 2\n",
    "pcen_version_str = str(pcen_version_id)\n",
    "hdf5_dir = os.path.join(\n",
    "    ccb_dir, \"ccb18_h5_v-\" + pcen_version_str)\n",
    "distances = []\n",
    "dates = []\n",
    "channels = []\n",
    "\n",
    "stft_max_list = []\n",
    "stft_avg_list = []\n",
    "log_pcen_max_list = []\n",
    "log_pcen_avg_list = []\n",
    "\n",
    "\n",
    "df = pd.read_csv(csv_path)\n",
    "df = df.dropna()\n",
    "df = df.sort_values(by=\"DistanceKm\")\n",
    "\n",
    "\n",
    "import h5py\n",
    "\n",
    "from librosa.display import specshow\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import tqdm\n",
    "\n",
    "data_dir = \"/beegfs/vl1019/waspaa2019_data\"\n",
    "ccb_dir = os.path.join(data_dir, \"ccb18\")\n",
    "csv_name = 'CCB_2009Feb18to22Apr17.csv'\n",
    "csv_path = os.path.join(ccb_dir, csv_name)\n",
    "\n",
    "hop_length = 2**7\n",
    "n_fft = 2 * hop_length\n",
    "sr = 2000\n",
    "half_clip_n_cols = 16\n",
    "low_freq_bin = 5\n",
    "\n",
    "pcen_version_id = 2\n",
    "pcen_version_str = str(pcen_version_id)\n",
    "hdf5_dir = os.path.join(\n",
    "    ccb_dir, \"ccb18_h5_v-\" + pcen_version_str)\n",
    "distances = []\n",
    "dates = []\n",
    "channels = []\n",
    "\n",
    "stft_max_list = []\n",
    "stft_avg_list = []\n",
    "log_pcen_max_list = []\n",
    "log_pcen_avg_list = []\n",
    "\n",
    "\n",
    "df = pd.read_csv(csv_path)\n",
    "df = df.dropna()\n",
    "df = df.sort_values(by=\"DistanceKm\")\n",
    "\n",
    "for i, row_id in tqdm.tqdm(enumerate(range(len(df)))):\n",
    "    row = df.iloc[row_id]\n",
    "    \n",
    "    # Read date.\n",
    "    date_int = int(row[\"Date\"])\n",
    "    date_str = \"20090\" + str(date_int)\n",
    "    prefix_str = \"CCB18_\" + date_str\n",
    "\n",
    "    # Read channel.\n",
    "    channel_id = int(row[\"Channel\"])\n",
    "    channel_str = str(channel_id)\n",
    "\n",
    "    # Read onset time\n",
    "    onset_time_int = int(row[\"Begin.Time..s.\"])\n",
    "    hour_str = str(onset_time_int // 3600).zfill(2)\n",
    "    minute_str = str(((onset_time_int % 3600) // (60*15)) * 15).zfill(2)\n",
    "    second_str = \"00\"\n",
    "    time_str = \"\".join([hour_str, minute_str, second_str])\n",
    "    \n",
    "    # Locate HDF5 \n",
    "    hdf5_name = \"_\".join([\n",
    "        \"CCB18\",\n",
    "        date_str,\n",
    "        time_str,\n",
    "        pcen_version_str\n",
    "    ])\n",
    "    hdf5_path = os.path.join(hdf5_dir, hdf5_name + \".h5\")\n",
    "    \n",
    "    if not os.path.exists(hdf5_path):\n",
    "        continue\n",
    "        \n",
    "    distances.append(row[\"DistanceKm\"])\n",
    "    dates.append(row[\"Date\"])\n",
    "    channels.append(row[\"Channel\"])\n",
    "    \n",
    "    # Open HDF5\n",
    "    with h5py.File(hdf5_path, 'r') as f:\n",
    "        mid_time = (0.5 * (row[\"Begin.Time..s.\"] + row[\"End.Time..s.\"])) % (15*60)\n",
    "        mid_col = int(round(mid_time * (sr/hop_length)))\n",
    "        start_col = max(0, mid_col - half_clip_n_cols)\n",
    "        stop_col = min(start_col + 2*half_clip_n_cols, f[\"stft\"].shape[1])\n",
    "        start_col = stop_col - 2*half_clip_n_cols\n",
    "        X_stft = f[\"stft\"][low_freq_bin:, start_col:stop_col, channel_id-1]\n",
    "        X_stft_sf = np.maximum(0, np.diff(np.log1p(X_stft), axis=1))\n",
    "        stft_max_list.append(np.max(X_stft_sf))\n",
    "        stft_avg_list.append(np.max(np.mean(X_stft_sf, axis=0)))\n",
    "        X_pcen = f[\"pcen\"][low_freq_bin:, start_col:stop_col, channel_id-1]\n",
    "        log_pcen_max_list.append(np.max(np.log1p(X_pcen)))\n",
    "        log_pcen_avg_list.append(np.max(np.mean(np.log1p(X_pcen), axis=0)))\n",
    "\n",
    "feature_dict = {\n",
    "    \"stft_avg\": stft_avg_list,\n",
    "    \"stft_max\": stft_max_list,\n",
    "    \"log_pcen_avg\": log_pcen_avg_list,\n",
    "    \"log_pcen_max\": log_pcen_max_list,\n",
    "    \"distance\": distances,\n",
    "    \"date\": dates,\n",
    "    \"channel\": channels\n",
    "}\n",
    "feature_df = pd.DataFrame(feature_dict)\n",
    "\n",
    "feature_df.to_csv(\n",
    "    \"/beegfs/vl1019/waspaa2019_data/ccb18/ccb18_features_v-2bis\" +\\\n",
    "    pcen_version_str + \".csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import librosa\n",
    "import soundfile as sf\n",
    "\n",
    "ships_dir = '/beegfs/vl1019/waspaa2019_data/shipsEar_AUDIOS'\n",
    "settings = {\n",
    "    \"T\": 1.0,\n",
    "    \"alpha\": 1.0,\n",
    "    \"delta\": 0.0,\n",
    "    \"r\": 1.0,\n",
    "    \"eps\": 1e-6,\n",
    "    \"n_fft\": 2**8,\n",
    "    \"win_length\": 2**8,\n",
    "    \"hop_length\": 2**7,\n",
    "    \"sr\": 2000\n",
    "}\n",
    "\n",
    "stft_avg_neg = []\n",
    "stft_max_neg = []\n",
    "log_pcen_avg_neg = []\n",
    "log_pcen_max_neg = []\n",
    "\n",
    "\n",
    "for ship_path in tqdm.tqdm(glob.glob(os.path.join(ships_dir, \"*.wav\"))):\n",
    "\n",
    "    ship_waveform, orig_sr = sf.read(ship_path)\n",
    "\n",
    "    ship_waveform = librosa.resample(ship_waveform, orig_sr, settings[\"sr\"])\n",
    "\n",
    "    stft = librosa.stft(\n",
    "        ship_waveform * (2**31),\n",
    "        n_fft=settings[\"n_fft\"],\n",
    "        win_length=settings[\"win_length\"],\n",
    "        hop_length=settings[\"hop_length\"],\n",
    "        window=\"hann\")[low_freq_bin:, :]\n",
    "\n",
    "    logspec = np.log1p(np.abs(stft))\n",
    "\n",
    "    logspec_flux = np.maximum(0, np.diff(logspec, axis=1))\n",
    "    stft_avg_neg.append(np.mean(logspec_flux, axis=0))\n",
    "    stft_max_neg.append(np.max(logspec_flux, axis=0))\n",
    "\n",
    "    pcen = librosa.pcen(np.abs(stft),\n",
    "        sr=settings[\"sr\"],\n",
    "        hop_length=settings[\"hop_length\"],\n",
    "        gain=settings[\"alpha\"],\n",
    "        bias=settings[\"delta\"],\n",
    "        power=settings[\"r\"],\n",
    "        time_constant=settings[\"T\"],\n",
    "        eps=settings[\"eps\"])[low_freq_bin:, :]\n",
    "\n",
    "    log_pcen_avg_neg.append(np.mean(np.log1p(pcen), axis=0))\n",
    "    log_pcen_max_neg.append(np.max(np.log1p(pcen), axis=0))\n",
    "    \n",
    "    \n",
    "neg_dict = {\n",
    "    \"log_pcen_avg\": np.concatenate([x[10:-10] for x in log_pcen_avg_neg]),\n",
    "    \"log_pcen_max\": np.concatenate([x[10:-10] for x in log_pcen_max_neg]),\n",
    "    \"stft_avg\": np.concatenate([x[10:-10] for x in stft_avg_neg]),\n",
    "    \"stft_max\": np.concatenate([x[10:-10] for x in stft_max_neg])\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('/beegfs/vl1019/waspaa2019_data/ccb18/ccb18_negatives_v-' + pcen_version_str + \"bis.hdf5\", 'w') as f:\n",
    "    for k in neg_dict.keys():\n",
    "        f[k] = neg_dict[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "neg_dict = {}\n",
    "\n",
    "with h5py.File('/beegfs/vl1019/waspaa2019_data/ccb18/ccb18_negatives_v-' + pcen_version_str + \"bis.hdf5\", 'r') as f:\n",
    "    for k in f.keys():\n",
    "        neg_dict[k] = f[k][()]\n",
    "        \n",
    "        \n",
    "feature_df = pd.DataFrame.from_csv(\n",
    "    \"/beegfs/vl1019/waspaa2019_data/ccb18/ccb18_features_v-2bis\" +\\\n",
    "    pcen_version_str + \".csv\")\n",
    "\n",
    "max_distance = 21\n",
    "command_tpr = 0.5\n",
    "dates = [219, 220, 221, 222, 417]\n",
    "feature_mtbfs = []\n",
    "\n",
    "feature_strs = [\"log_pcen_max\", \"stft_avg\", \"stft_max\", \"log_pcen_avg\"]\n",
    "\n",
    "for feature_str in feature_strs:\n",
    "\n",
    "    distance_mtbfs = []\n",
    "\n",
    "    for min_distance in range(30):\n",
    "        dist_rows =\\\n",
    "                (feature_df[\"distance\"] >= min_distance) &\\\n",
    "                (feature_df[\"distance\"] < (min_distance+1)) &\\\n",
    "                (feature_df[\"channel\"] != 3)\n",
    "\n",
    "        mtbfs = []\n",
    "        for date in dates:\n",
    "            date_dist_rows = dist_rows & (feature_df[\"date\"] == date)\n",
    "            sorted_values = np.sort(feature_df[date_dist_rows][feature_str])[::-1]\n",
    "            n_positives = len(sorted_values)\n",
    "            command_threshold = sorted_values[int(command_tpr*n_positives)]\n",
    "            command_fpr = np.mean(neg_dict[feature_str] > command_threshold)\n",
    "            mtbf = (hop_length/sr) / command_fpr\n",
    "            mtbfs.append(mtbf)\n",
    "\n",
    "        distance_mtbfs.append(np.stack(mtbfs))\n",
    "        \n",
    "    feature_mtbfs.append(np.stack(distance_mtbfs))\n",
    "    \n",
    "feature_mtbfs = np.array(feature_mtbfs)\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "plt.title('Bioacoustic detection of North Atlantic Right Whale', size=12)\n",
    "\n",
    "plt.plot(\n",
    "    np.log10(np.median(feature_mtbfs[0, :max_distance, :], axis=1)),\n",
    "    '-o', markersize=10.0, linewidth=2.0, color=\"#009900\",\n",
    "    label=\"Max-pooled PCEN\")\n",
    "plt.fill_between(\n",
    "    range(max_distance),\n",
    "    np.log10(np.ravel(np.quantile(feature_mtbfs[0, :max_distance, :], [0.25], axis=1))),\n",
    "    np.log10(np.ravel(np.quantile(feature_mtbfs[0, :max_distance, :], [0.75], axis=1))),\n",
    "    alpha = 0.5,\n",
    "    color=\"#009900\",\n",
    ")\n",
    "\n",
    "\n",
    "plt.plot(\n",
    "    np.log10(np.median(feature_mtbfs[2, :max_distance, :], axis=1)),\n",
    "    '-s', color=\"#0000B2\",\n",
    "    markersize=10.0, linewidth=2.0,\n",
    "    label=\"Max-pooled spectral flux\")\n",
    "plt.fill_between(\n",
    "    range(max_distance),\n",
    "    np.log10(np.ravel(np.quantile(feature_mtbfs[2, :max_distance, :], [0.25], axis=1))),\n",
    "    np.log10(np.ravel(np.quantile(feature_mtbfs[2, :max_distance, :], [0.75], axis=1))),\n",
    "    alpha = 0.5,\n",
    "    color=\"#0000B2\"\n",
    ")\n",
    "\n",
    "\n",
    "plt.plot(\n",
    "    np.log10(np.median(feature_mtbfs[1, :max_distance, :], axis=1)), '-v', markersize=10.0, linewidth=2.0,\n",
    "    color=\"#E67300\", label=\"Averaged spectral flux\")\n",
    "plt.fill_between(\n",
    "    range(max_distance),\n",
    "    np.log10(np.ravel(np.quantile(feature_mtbfs[1, :max_distance, :], [0.25], axis=1))),\n",
    "    np.log10(np.ravel(np.quantile(feature_mtbfs[1, :max_distance, :], [0.75], axis=1))),\n",
    "    alpha = 0.33, color=\"#E67300\"\n",
    ")\n",
    "\n",
    "yticks = np.array([0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0])\n",
    "\n",
    "\n",
    "plt.rcParams[\"font.family\"] = \"serif\"\n",
    "plt.legend(prop={'size': 11})\n",
    "plt.xlabel(\"CCB18 dataset: distance between sensor and source (km)\", size=12)\n",
    "plt.ylabel(\"ShipsEar dataset: mean time\\nbetween false alarms at half recall (s)\", size=12)\n",
    "plt.gca().set_xticks(np.arange(0, 21, 2), minor=True)\n",
    "plt.gca().set_xticks(np.arange(0, 21, 4), minor=False)\n",
    "plt.gca().set_xticklabels(np.arange(0, 21, 4), size=11)\n",
    "plt.yticks(np.log10(yticks))\n",
    "plt.gca().set_yticklabels(yticks, size=11)\n",
    "plt.xlim(0, 20)\n",
    "plt.grid(linestyle='--', alpha=1.0, which=\"minor\")\n",
    "plt.grid(linestyle='--', alpha=1.0, which=\"major\")\n",
    "ax = plt.gca()\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "plt.savefig('lostanlen_waspaa2019_ccb_mtbfa-50_semilogy.eps', bbox_inches='tight')\n",
    "plt.savefig('lostanlen_waspaa2019_ccb_mtbfa-50_semilogy.png', bbox_inches='tight', dpi=500)"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}