5_DissimilarityMatrix.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e52d338c-6919-4b73-98d8-0158c2f482f4",
   "metadata": {},
   "source": [
    "# Compute Fused Dissimilarity Matrix\n",
    "Compute geodesic fused angle-delay profile (ADP) / timestamp dissimilarity matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0eeabd91-dbae-419c-9007-64c201325893",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-02-06 09:32:04.085052: I tensorflow/core/util/port.cc:110] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.\n",
      "2024-02-06 09:32:04.112754: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.\n",
      "To enable the following instructions: AVX2 AVX_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
      "2024-02-06 09:32:04.555105: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT\n",
      "2024-02-06 09:32:05.147590: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.162686: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.162803: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.164623: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.164725: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.164799: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.487342: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.487453: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.487526: I tensorflow/compiler/xla/stream_executor/cuda/cuda_gpu_executor.cc:995] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero. See more at https://github.com/torvalds/linux/blob/v6.0/Documentation/ABI/testing/sysfs-bus-pci#L344-L355\n",
      "2024-02-06 09:32:05.487585: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1639] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 224 MB memory:  -> device: 0, name: NVIDIA GeForce GTX 1080 Ti, pci bus id: 0000:01:00.0, compute capability: 6.1\n",
      "2024-02-06 09:32:05.500447: I tensorflow/compiler/xla/stream_executor/cuda/cuda_driver.cc:753] failed to allocate 224.06MiB (234946560 bytes) from device: CUDA_ERROR_OUT_OF_MEMORY: out of memory\n",
      "2024-02-06 09:32:05.893766: I tensorflow/core/grappler/optimizers/data/replicate_on_split.cc:32] Running replicate on split optimization\n",
      "2024-02-06 09:32:05.914956: I tensorflow/core/grappler/optimizers/data/replicate_on_split.cc:32] Running replicate on split optimization\n"
     ]
    }
   ],
   "source": [
    "from dichasus_cf0x import training_set\n",
    "import multiprocessing as mp\n",
    "import scipy.sparse.csgraph\n",
    "import sklearn.neighbors\n",
    "import tensorflow as tf\n",
    "import scipy.spatial\n",
    "import numpy as np\n",
    "import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "16c6a961-92c3-4173-a474-9c2d358dff70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def csi_time_domain(csi, pos, time):\n",
    "    csi = tf.signal.fftshift(tf.signal.ifft(tf.signal.fftshift(csi, axes = -1)), axes = -1)\n",
    "    return csi, pos, time\n",
    "\n",
    "def cut_out_taps(tap_start, tap_stop):\n",
    "    def cut_out_taps_func(csi, pos, time):\n",
    "        return csi[:,:,:,tap_start:tap_stop], pos, time\n",
    "\n",
    "    return cut_out_taps_func\n",
    "\n",
    "\n",
    "training_set = training_set.map(csi_time_domain, num_parallel_calls = tf.data.AUTOTUNE)\n",
    "training_set = training_set.map(cut_out_taps(507, 520), num_parallel_calls = tf.data.AUTOTUNE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84a36010-50d0-46bf-80b8-b5f68d8d8329",
   "metadata": {},
   "source": [
    "### Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0013363e-49b5-4292-abee-dbace51faa67",
   "metadata": {},
   "outputs": [],
   "source": [
    "classical_estimated_positions = np.load(\"results/estimated_positions.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5cb15d7-9a6e-4015-9010-76624de34d11",
   "metadata": {},
   "source": [
    "# From TensorFlow to NumPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6c378b60-fce8-42ea-8fcf-ca326e98ac72",
   "metadata": {},
   "outputs": [],
   "source": [
    "groundtruth_positions = []\n",
    "csi_time_domain = []\n",
    "timestamps = []\n",
    "\n",
    "for csi, pos, time in training_set.prefetch(tf.data.AUTOTUNE).batch(1000):\n",
    "    csi_time_domain.append(csi.numpy())\n",
    "    groundtruth_positions.append(pos.numpy())\n",
    "    timestamps.append(time.numpy())\n",
    "\n",
    "csi_time_domain = np.concatenate(csi_time_domain)\n",
    "groundtruth_positions = np.concatenate(groundtruth_positions)\n",
    "timestamps = np.concatenate(timestamps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2323add4-59f8-422d-b396-39eedf539f2d",
   "metadata": {},
   "source": [
    "# Dissimilarity Metrics and Dissimilarity Matrix Computation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd029ddd-29ef-4aea-a936-f98225c8ba35",
   "metadata": {},
   "source": [
    "### Step 1: ADP-based dissimilarity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b18f625e-972e-4ba5-8000-b1cf1f75f71c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 439866729/439866729 [00:58<00:00, 7493696.88it/s] \n"
     ]
    }
   ],
   "source": [
    "adp_dissimilarity_matrix = np.zeros((csi_time_domain.shape[0], csi_time_domain.shape[0]), dtype=np.float32)\n",
    "\n",
    "def adp_dissimilarities_worker(todo_queue, output_queue):\n",
    "    def adp_dissimilarities(index):\n",
    "        # h has shape (arrays, antenna rows, antenna columns, taps), w has shape (datapoints, arrays, antenna rows, antenna columns, taps)\n",
    "        h = csi_time_domain[index]\n",
    "        w = csi_time_domain[index:]\n",
    "\n",
    "        dotproducts = np.abs(np.einsum(\"brmt,lbrmt->lbt\", np.conj(h), w))**2\n",
    "        norms = np.real(np.einsum(\"brmt,brmt->bt\", h, np.conj(h)) * np.einsum(\"lbrmt,lbrmt->lbt\", w, np.conj(w)))\n",
    "        \n",
    "        return np.sum(1 - dotproducts / norms, axis = (1, 2))\n",
    "\n",
    "    while True:\n",
    "        index = todo_queue.get()\n",
    "\n",
    "        if index == -1:\n",
    "            output_queue.put((-1, None))\n",
    "            break\n",
    "        \n",
    "        output_queue.put((index, adp_dissimilarities(index)))\n",
    "\n",
    "with tqdm.tqdm(total = csi_time_domain.shape[0]**2) as bar:\n",
    "    todo_queue = mp.Queue()\n",
    "    output_queue = mp.Queue()\n",
    "\n",
    "    for i in range(csi_time_domain.shape[0]):\n",
    "        todo_queue.put(i)\n",
    "    \n",
    "    for i in range(mp.cpu_count()):\n",
    "        todo_queue.put(-1)\n",
    "        p = mp.Process(target = adp_dissimilarities_worker, args = (todo_queue, output_queue))\n",
    "        p.start()\n",
    "\n",
    "    finished_processes = 0\n",
    "    while finished_processes != mp.cpu_count():\n",
    "        i, d = output_queue.get()\n",
    "\n",
    "        if i == -1:\n",
    "            finished_processes = finished_processes + 1\n",
    "        else:\n",
    "            adp_dissimilarity_matrix[i,i:] = d\n",
    "            adp_dissimilarity_matrix[i:,i] = d\n",
    "            bar.update(2 * len(d) - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e90c51df-55f5-4438-bb5d-848a522163ec",
   "metadata": {},
   "source": [
    "### Step 2: Timestamp-based dissimilarity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "54929cab-7bf5-40cb-8cd3-a6ccc49bbb85",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute timestamp-based dissimilarity matrix\n",
    "timestamp_dissimilarity_matrix = np.abs(np.subtract.outer(timestamps, timestamps))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67c37817-a709-4fbd-b97a-a05f7f3c151a",
   "metadata": {},
   "source": [
    "### Step 3: Fusing $\\mathbf D_\\mathrm{ADP}$ with $\\mathbf D_\\mathrm{time}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d44a378a-a58f-4867-8e14-539d21b5d360",
   "metadata": {},
   "outputs": [],
   "source": [
    "TIME_THRESHOLD = 2\n",
    "small_time_dissimilarity_indices = np.logical_and(timestamp_dissimilarity_matrix < TIME_THRESHOLD, timestamp_dissimilarity_matrix > 0)\n",
    "small_time_dissimilarities = timestamp_dissimilarity_matrix[small_time_dissimilarity_indices]\n",
    "small_adp_dissimilarities = adp_dissimilarity_matrix[small_time_dissimilarity_indices]\n",
    "\n",
    "occurences, edges = np.histogram(small_adp_dissimilarities / small_time_dissimilarities, range = (0, 50), bins = 1500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2659b2fc-394a-403f-8790-209b91d961d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_centers = edges[:-1] + np.diff(edges) / 2.\n",
    "max_bin = np.argmax(occurences)\n",
    "min_threshold = np.quantile(occurences[:max_bin], 0.5)\n",
    "\n",
    "for threshold_bin in range(max_bin - 1, -1, -1):\n",
    "    if occurences[threshold_bin] < min_threshold:\n",
    "        break\n",
    "\n",
    "scaling_factor = bin_centers[threshold_bin]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6e6b3780-238d-4f2b-b88b-4b162feabfcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fuse ADP-based and time-based dissimilarity matrices\n",
    "dissimilarity_matrix_fused = np.minimum(adp_dissimilarity_matrix, timestamp_dissimilarity_matrix * scaling_factor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdca4a35-7f73-48e3-8f40-c62afa873b42",
   "metadata": {},
   "source": [
    "### Step 4: Geodesic Dissimilarity Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "32b34306-0c96-46ef-94e9-e7bed2d7beb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_neighbors = 20\n",
    "\n",
    "nbrs_alg = sklearn.neighbors.NearestNeighbors(n_neighbors = n_neighbors, metric=\"precomputed\", n_jobs = -1)\n",
    "nbrs = nbrs_alg.fit(dissimilarity_matrix_fused)\n",
    "nbg = sklearn.neighbors.kneighbors_graph(nbrs, n_neighbors, metric = \"precomputed\", mode=\"distance\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a6150477-9dcc-4abc-b667-03bb0cad53bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 439866729/439866729 [00:15<00:00, 28272258.03it/s] \n"
     ]
    }
   ],
   "source": [
    "dissimilarity_matrix_geodesic = np.zeros((nbg.shape[0], nbg.shape[1]), dtype = np.float32)\n",
    "\n",
    "def shortest_path_worker(todo_queue, output_queue):\n",
    "    while True:\n",
    "        index = todo_queue.get()\n",
    "\n",
    "        if index == -1:\n",
    "            output_queue.put((-1, None))\n",
    "            break\n",
    "\n",
    "        d = scipy.sparse.csgraph.dijkstra(nbg, directed=False, indices=index)\n",
    "        output_queue.put((index, d))\n",
    "\n",
    "with tqdm.tqdm(total = nbg.shape[0]**2) as bar:\n",
    "    todo_queue = mp.Queue()\n",
    "    output_queue = mp.Queue()\n",
    "\n",
    "    for i in range(nbg.shape[0]):\n",
    "        todo_queue.put(i)\n",
    "    \n",
    "    for i in range(mp.cpu_count()):\n",
    "        todo_queue.put(-1)\n",
    "        p = mp.Process(target = shortest_path_worker, args = (todo_queue, output_queue))\n",
    "        p.start()\n",
    "\n",
    "    finished_processes = 0\n",
    "    while finished_processes != mp.cpu_count():\n",
    "        i, d = output_queue.get()\n",
    "\n",
    "        if i == -1:\n",
    "            finished_processes = finished_processes + 1\n",
    "        else:\n",
    "            dissimilarity_matrix_geodesic[i,:] = d\n",
    "            bar.update(len(d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbbf01d2-8e63-44cd-a78b-79f028861f2f",
   "metadata": {},
   "source": [
    "# Scaling the Dissilimarity Matrix\n",
    "All values should be in meters..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "65fdf4cc-362c-4aec-8740-69a10130c4ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaling_nth_reduction = 30\n",
    "classical_positions_reduced = classical_estimated_positions[::scaling_nth_reduction]\n",
    "dissimilarity_matrix_reduced = dissimilarity_matrix_geodesic[::scaling_nth_reduction, ::scaling_nth_reduction]\n",
    "classical_distance_matrix = scipy.spatial.distance_matrix(classical_positions_reduced, classical_positions_reduced)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "e5a8fbef-5d63-4ec6-b09c-b5f0ffd979ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "dissimilarity_unit_meters = np.full_like(dissimilarity_matrix_reduced, np.nan)\n",
    "diff = np.divide(dissimilarity_matrix_reduced, classical_distance_matrix, out = dissimilarity_unit_meters, where = classical_distance_matrix != 0)\n",
    "dissimilarity_unit_meters = dissimilarity_unit_meters.flatten()\n",
    "scaling_factor_meters = np.median(dissimilarity_unit_meters[np.isfinite(dissimilarity_unit_meters)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d0dbb0dc-e1f1-44dc-bdf3-ead3989b6719",
   "metadata": {},
   "outputs": [],
   "source": [
    "dissimilarity_matrix_geodesic_meters = dissimilarity_matrix_geodesic / scaling_factor_meters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "17ea3577-6d31-430c-8799-bb97fdb2780e",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(\"results/dissimilarity_matrix_geodesic_meters.npy\", dissimilarity_matrix_geodesic_meters)"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}