{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "4c9f2995", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "--- PROCESSING SESSION: 085315 ---\n", "Running track association...\n", "Executing spectral fusion...\n", "Session 085315 successfully clustered and saved.\n", "\n", "--- PROCESSING SESSION: 081356 ---\n", "Running track association...\n", "Executing spectral fusion...\n", "Session 081356 successfully clustered and saved.\n", "\n", "--- PROCESSING SESSION: 085452 ---\n", "Running track association...\n", "Executing spectral fusion...\n", "Session 085452 successfully clustered and saved.\n", "\n", "--- BATCH CLUSTERING COMPLETE ---\n", "Note: output columns are 'fragment_id' (pre-fusion) and 'cluster_merged' (post-fusion).\n" ] } ], "source": [ "\"\"\"\n", "\n", "Module 03: Click train spectral clustering\n", "Journal of the Acoustical Society of America (JASA) - Data Availability\n", "\n", "(c) Lara Berkenbaum and Hervé Glotin (2025, 2026)\n", "Original code produced for the CIAN Tool Kit, associated with the manuscript: \n", "Berkenbaum, L., Glotin, H., Sarano, F., Zimmer, W.M.X., Sarano, V., Adam, O., Giraudet, P. (2026). \"Audiovisual diarization of overlapping click trains in sperm whale (Physeter macrocephalus) vocal sparring using a three-hydrophone array,\" J. Acoust. Soc. Am.\n", "\n", "Description:\n", "This script performs unsupervised clustering on the extracted CQT patches to group sperm whale clicks into coherent trains. \n", "It utilizes a tracking algorithm with short-term and long-term template memory, compensating for temporal jitter via cross-correlation shifts. \n", "Processes multiple recording sessions dynamically.\n", "\n", "Inputs: \n", "- CQT patches (Pickle): `../01_data/ground_truth_acoustic/*_CQT_patch_H0.pkl`\n", "- Spectral features (CSV): `../01_data/ground_truth_acoustic/*_CQT_feat_H0.csv`\n", "\n", "Outputs: \n", "- Click clustered: `../01_data/ground_truth_acoustic/*_click_clustering.csv`\n", "\n", "\"\"\"\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import pickle\n", "import os\n", "import glob\n", "import re\n", "from scipy.spatial.distance import cosine\n", "\n", "# ======================================================================\n", "# 1. CONFIGURATION & BATCH SETTINGS\n", "# ======================================================================\n", "csv_dir = \"../01_data/ground_truth_acoustic/\"\n", "\n", "# Tracking thresholds\n", "CORR_THRESH_HIGH = 0.60\n", "CORR_THRESH_LOW = 0.35\n", "ICI_TOLERANCE = 0.25\n", "MAX_ICI = 2.50\n", "MIN_TRACK_LEN = 2\n", "\n", "# Template and parameters\n", "TEMPLATE_FREEZE_LEN = 4\n", "WEIGHT_COSINE = 0.3 \n", "\n", "# ======================================================================\n", "# 2. UTILITY FUNCTION: ROBUST JITTER CORRELATION\n", "# ======================================================================\n", "def robust_correlation(sig_a, sig_b):\n", " \"\"\"\n", " Computes cross-correlation while accounting for minor temporal jitter (-1, 0, +1 shifts).\n", " Prevents track fragmentation due to sub-sample misalignment.\n", " \"\"\"\n", " score_0 = np.corrcoef(sig_a, sig_b)[0,1]\n", " \n", " sig_b_left = np.roll(sig_b, -1)\n", " score_L = np.corrcoef(sig_a, sig_b_left)[0,1]\n", " \n", " sig_b_right = np.roll(sig_b, 1)\n", " score_R = np.corrcoef(sig_a, sig_b_right)[0,1]\n", " \n", " return max(score_0, score_L, score_R)\n", "\n", "# ======================================================================\n", "# 3. CLASS: SPERM WHALE TRACKER\n", "# ======================================================================\n", "class SpermWhaleTracker:\n", " def __init__(self, track_id, start_idx, start_time, start_sig):\n", " self.track_id = track_id\n", " self.indices = [start_idx]\n", " self.times = [start_time]\n", " self.signatures = [start_sig]\n", " self.active = True\n", " self.last_ici = 0.5 \n", " self.ici_history = [] \n", " self.template = start_sig \n", "\n", " def predict_window(self):\n", " last_t = self.times[-1]\n", " if len(self.times) == 1:\n", " return last_t + 0.04, last_t + 1.5\n", " \n", " current_ici = np.median(self.ici_history[-5:]) if self.ici_history else self.last_ici\n", " t_min = last_t + (current_ici * (1 - ICI_TOLERANCE))\n", " t_max = last_t + (current_ici * (1 + ICI_TOLERANCE))\n", " return t_min, t_max\n", "\n", " def score_candidate(self, cand_time, cand_sig):\n", " # A. Long-term memory comparison (template)\n", " corr_template = robust_correlation(self.template, cand_sig)\n", " \n", " # B. Short-term memory comparison (last click)\n", " corr_last = robust_correlation(self.signatures[-1], cand_sig)\n", " \n", " corr_val = max(corr_template, corr_last)\n", " cos_sim = 1 - cosine(self.template, cand_sig)\n", " \n", " vis_score = (1 - WEIGHT_COSINE) * corr_val + WEIGHT_COSINE * cos_sim\n", " \n", " # Temporal penalty formulation\n", " time_penalty = 0\n", " if len(self.times) > 1:\n", " current_ici = np.median(self.ici_history[-5:])\n", " expected = self.times[-1] + current_ici\n", " dt_ratio = abs(cand_time - expected) / current_ici\n", " time_penalty = min(dt_ratio**2, 1.0) \n", " \n", " return vis_score, time_penalty\n", "\n", " def add_click(self, idx, time, sig, vis_score):\n", " if len(self.times) > 0:\n", " self.ici_history.append(time - self.times[-1])\n", " \n", " self.indices.append(idx)\n", " self.times.append(time)\n", " self.signatures.append(sig)\n", " \n", " if len(self.indices) > TEMPLATE_FREEZE_LEN:\n", " if vis_score > 0.6:\n", " self.template = (0.8 * self.template) + (0.2 * sig)\n", "\n", "# ======================================================================\n", "# 4. BATCH PROCESSING PIPELINE\n", "# ======================================================================\n", "feature_files = glob.glob(os.path.join(csv_dir, \"*_CQT_feat_H0.csv\"))\n", "\n", "if not feature_files:\n", " print(f\"Error: No feature datasets found in {csv_dir}.\")\n", "\n", "for csv_filename in feature_files:\n", " \n", " # Session extraction: anchoring the regex to the date format YYYYMMDD_HHMMSS\n", " base_name = os.path.basename(csv_filename).replace(\"_CQT_feat_H0.csv\", \"\")\n", " match = re.search(r'\\d{8}_(\\d{6})', base_name)\n", " \n", " if not match:\n", " print(f\"Warning: could not extract session ID from {base_name}. Skipping.\")\n", " continue\n", " \n", " session_id = match.group(1)\n", " \n", " print(f\"\\n--- PROCESSING SESSION: {session_id} ---\")\n", " \n", " # Identify corresponding pickle file dynamically\n", " pkl_candidates = glob.glob(os.path.join(csv_dir, f\"*{session_id}*_CQT_patch_H0.pkl\"))\n", " if not pkl_candidates:\n", " print(f\"Error: missing PKL patch file for session {session_id}. Skipping.\")\n", " continue\n", " file_pkl = pkl_candidates[0]\n", "\n", " # Data loading and preparation\n", " df = pd.read_csv(csv_filename, sep=';')\n", " \n", " if 'T_H0' not in df.columns:\n", " print(f\"Error: temporal column 'T_H0' not found in features {csv_filename}. Check Module 02 export. Skipping.\")\n", " continue\n", "\n", " df = df.sort_values('T_H0').reset_index(drop=True)\n", "\n", " with open(file_pkl, \"rb\") as f:\n", " data_pkl = pickle.load(f)\n", " \n", " pmap = {p[\"index\"]: p[\"patch\"] for p in data_pkl[\"patches\"]}\n", "\n", " signatures = {}\n", " for idx_csv, row in df.iterrows():\n", " orig_idx = int(row['index']) if pd.notna(row['index']) else -1\n", " \n", " if orig_idx in pmap:\n", " p_db = 10 * np.log10(np.maximum(pmap[orig_idx], 1e-10))\n", " p_flat = p_db.flatten()\n", " if np.std(p_flat) > 0:\n", " signatures[idx_csv] = (p_flat - np.mean(p_flat)) / np.std(p_flat)\n", "\n", " # ------------------------------------------------------------------\n", " # Tracking execution\n", " # ------------------------------------------------------------------\n", " print(\"Running track association...\")\n", " active_tracks = []\n", " completed_tracks = []\n", " assigned_indices = set()\n", " global_track_counter = 1\n", "\n", " for i, row in df.iterrows():\n", " if i not in signatures: continue\n", " current_time = row['T_H0']\n", " current_sig = signatures[i]\n", " \n", " candidates = []\n", " for track in active_tracks:\n", " t_min, t_max = track.predict_window()\n", " \n", " if current_time > t_max + 1.5:\n", " track.active = False\n", " continue\n", " \n", " if current_time >= t_min and current_time <= t_max:\n", " vis, pen = track.score_candidate(current_time, current_sig)\n", " score = vis - (0.5 * pen)\n", " \n", " is_valid = False\n", " if vis > CORR_THRESH_HIGH: is_valid = True\n", " elif vis > CORR_THRESH_LOW and pen < 0.1: is_valid = True\n", " \n", " if is_valid:\n", " candidates.append((score, track, vis))\n", " \n", " if candidates:\n", " candidates.sort(key=lambda x: x[0], reverse=True)\n", " best_score, best_track, best_vis = candidates[0]\n", " best_track.add_click(i, current_time, current_sig, best_vis)\n", " assigned_indices.add(i)\n", " else:\n", " if i not in assigned_indices:\n", " new_track = SpermWhaleTracker(global_track_counter, i, current_time, current_sig)\n", " active_tracks.append(new_track)\n", " global_track_counter += 1\n", " assigned_indices.add(i)\n", "\n", " # Track maintenance\n", " still_active = []\n", " for t in active_tracks:\n", " if t.active: still_active.append(t)\n", " else:\n", " if len(t.indices) >= MIN_TRACK_LEN: completed_tracks.append(t)\n", " active_tracks = still_active\n", "\n", " for t in active_tracks:\n", " if len(t.indices) >= MIN_TRACK_LEN: completed_tracks.append(t)\n", "\n", " # ------------------------------------------------------------------\n", " # Advanced spectral fusion\n", " # ------------------------------------------------------------------\n", " print(\"Executing spectral fusion...\")\n", " tracks_meta = []\n", " for t in completed_tracks:\n", " ici_med = np.median(t.ici_history) if t.ici_history else 0\n", " tracks_meta.append({\n", " 'id': t.track_id,\n", " 'start': t.times[0],\n", " 'end': t.times[-1],\n", " 'ici': ici_med,\n", " 'end_sig': t.signatures[-1],\n", " 'start_sig': t.signatures[0]\n", " })\n", "\n", " meta_df = pd.DataFrame(tracks_meta).sort_values('start')\n", " merged_map = {row['id']: row['id'] for _, row in meta_df.iterrows()}\n", "\n", " MAX_GAP = 5.0\n", " MAX_ICI_DIFF = 0.04\n", " MIN_LINK_CORR = 0.55 \n", "\n", " changed = True\n", " while changed:\n", " changed = False\n", " curr_ids = sorted(list(set(merged_map.values())))\n", " \n", " for i in range(len(curr_ids)):\n", " id_a = curr_ids[i]\n", " constituents_a = [k for k,v in merged_map.items() if v == id_a]\n", " last_meta_a = meta_df[meta_df['id'].isin(constituents_a)].sort_values('end').iloc[-1]\n", " \n", " for j in range(i+1, len(curr_ids)):\n", " id_b = curr_ids[j]\n", " constituents_b = [k for k,v in merged_map.items() if v == id_b]\n", " first_meta_b = meta_df[meta_df['id'].isin(constituents_b)].sort_values('start').iloc[0]\n", " \n", " gap = first_meta_b['start'] - last_meta_a['end']\n", " \n", " if gap < 0.06 or gap > MAX_GAP: continue\n", " if abs(last_meta_a['ici'] - first_meta_b['ici']) > MAX_ICI_DIFF: continue\n", " \n", " sig_a = last_meta_a['end_sig']\n", " sig_b = first_meta_b['start_sig']\n", " \n", " link_corr = robust_correlation(sig_a, sig_b)\n", " \n", " if link_corr > MIN_LINK_CORR:\n", " for k, v in merged_map.items():\n", " if v == id_b: merged_map[k] = id_a\n", " changed = True\n", " break\n", " if changed: break\n", "\n", " # Output preparation\n", " df['fragment_id'] = -1\n", " for t in completed_tracks:\n", " df.loc[t.indices, 'fragment_id'] = t.track_id\n", " \n", " df['cluster_merged'] = df['fragment_id'].map(merged_map)\n", " df['cluster_merged'] = df['cluster_merged'].fillna(-1)\n", "\n", " out_file = os.path.join(csv_dir, f\"20230504_{session_id}_click_clustering.csv\")\n", " df.to_csv(out_file, sep=';', index=False)\n", " print(f\"Session {session_id} successfully clustered and saved.\")\n", "\n", "print(\"\\n--- BATCH CLUSTERING COMPLETE ---\")\n", "print(\"Note: output columns are 'fragment_id' (pre-fusion) and 'cluster_merged' (post-fusion).\")" ] }, { "cell_type": "code", "execution_count": null, "id": "94659ea8-a782-4524-8dba-9226f8b63c99", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.4" } }, "nbformat": 4, "nbformat_minor": 5 }