demos/relative_heading_angle_explorer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A notebook with a simple example of two datapoints of two hypothetical fish to explore relative angles between the two. For the implementation of algorithms that include relative positions between the fish and their heading trajectories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import atan2, degrees\n",
    "import numpy as np\n",
    "import gridtools as gt\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD8CAYAAABgtYFHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAASTElEQVR4nO3dfYxV9Z3H8fcXGBURHbATnAEaIIKGJWtxpi5G3aLU+LBOwUSNg0lJQzImdcXWSRV3s5Vt9o81Fp+y1jARW9rIrCw1wJithkVou5ss60xtFEWQ9REYZIwC1kco3/3jnmHvjDPOj3PuveecO59XQphz7rn3fj3Bd845c+aOuTsiIsMZlfYAIpIPioWIBFEsRCSIYiEiQRQLEQmiWIhIkGFjYWZPmNlBM9tRtG6imW02s9ejvydE683MHjGzPWb2kpldWM7hRaRyQo4sfgFcPWDdcmCLu88EtkTLANcAM6M/rcBjpRlTRNI2bCzc/XfABwNWLwTWRF+vARYVrf+lF/w3UGtm9SWaVURSNCbm8ya5e0/09QFgUvT1ZODdou32Rut6GMDMWikcfTBu3LjG888/P+YoIhKiu7v7fXevi/v8uLE4wd3dzE76nnF3bwfaAZqamryrqyvpKCLyFczs7STPj/vdkPf6Ti+ivw9G6/cBU4u2mxKtE5GcixuLTcCS6OslwMai9d+NvisyDzhcdLoiIjk27GmImXUA84Gvmdle4F7gn4F1ZrYUeBu4Kdr834FrgT3AJ8D3yjCziKRg2Fi4e8sQDy0YZFsHbks6lIhkj+7gFJEgioWIBFEsRCSIYiEiQRQLEQmiWIhIEMVCRIIoFiISRLEQkSCKhYgEUSxEJIhiISJBFAsRCaJYiEgQxUJEgigWIhJEsRCRIIqFiARRLEQkiGIhIkEUCxEJoliISBDFQkSCKBYiEkSxEJEgioWIBFEsRCSIYiEiQRQLEQmiWIhIEMVCRIIoFiISRLEQkSCKhYgEUSxEJEiiWJjZD83sFTPbYWYdZnaamU03s+1mtsfMnjKzU0o1rIikJ3YszGwysAxocvc5wGjgZuA+4EF3Pxf4EFhaikFFJF1JT0PGAGPNbAxwOtADXAGsjx5fAyxK+B4ikgGxY+Hu+4CfAu9QiMRhoBs45O7Hos32ApMHe76ZtZpZl5l19fb2xh1DRCokyWnIBGAhMB1oAMYBV4c+393b3b3J3Zvq6urijiEiFZLkNOTbwJvu3uvuR4GngUuA2ui0BGAKsC/hjCKSAUli8Q4wz8xONzMDFgCvAluBG6JtlgAbk40oIlmQ5JrFdgoXMv8AvBy9VjtwN3Cnme0BzgZWl2BOEUnZmOE3GZq73wvcO2D1G8BFSV5XRLJHd3CKSBDFQkSCKBYiEkSxEJEgioWIBFEsRCSIYiEiQRQLEQmiWIhIEMVCRIIoFiISRLEQkSCKhYgESfRTpyJZ9soHn/Hb/Z9w5OhxzqwZxbcaTucvJp6W9li5pVhIVXrlg8/4zTt/4pgXlo8cPc5v3vkTgIIRk05DpCr9dv8nJ0LR55gX1ks8ioVUpSNHj5/UehmeYiFVxd154YUX+Pj99wZ9/Mwa/ZOPS3tOcq8vEHfddRczZszgyiuvpPOBH3P82Bf9thtj8K2G01OaMv8UC8m19vZ2ZsyYwS233EJNTQ3z5s3jyJEjvP67Z7nynFM5s2YU7s6ZNaO45utn6OJmAoqF5Nrs2bPZsGEDu3bt4siRI6xbt45ly5Zx+eWX882pZ/P9ORP5fNOjfH/ORIUiIcVCcu3SSy/lggsuYNmyZfzsZz9j/fr1vPXWW9x4441pj1Z1FAvJvdtvv/1EKBYsWMDWrVtZuHDhicebm5tTnK56KBaSa8WhuP766+ns7OSyyy6jtrb2xDaNjY3pDVhFFAvJrYGhADjnnHO4++67+223YsWKFKarPrrdW3JpsFAALFiwIMWpqpuOLCR3hgqFlJdiIbkSJxSzZs0q81Qjg2IhuRH3iGLx4sVlnGrkUCwkF5Kceqxdu7ZMU40sioVkXtJrFLt37y7DVCOPYiGZpouZ2aFYSGYpFNmiWEgmlTIUuimrNBLFwsxqzWy9mb1mZjvN7GIzm2hmm83s9ejvCaUaVkaGUh9RdHd3l2AqSXpk8TDwrLufD1wA7ASWA1vcfSawJVoWCVKOU4/Ozs6SvM5IFzsWZnYW8NfAagB3/8LdDwELgTXRZmuARclGlJFC1yiyLcmRxXSgF/i5mb1oZo+b2Thgkrv3RNscACYN9mQzazWzLjPr6u3tTTCGVAOFIvuSxGIMcCHwmLvPBT5mwCmHuzvggzwXd2939yZ3b6qrq0swhuRduUPR0tJS8tcciZLEYi+w1923R8vrKcTjPTOrB4j+PphsRKlmlTiiaGhoKMvrjjSxY+HuB4B3zey8aNUC4FVgE7AkWrcE2JhoQqlalTr1WLlyZdleeyRJ+nkWtwNPmtkpwBvA9ygEaJ2ZLQXeBm5K+B5ShXSNIn8SxcLd/wg0DfKQPoFEhqRQ5JPu4JSKSiMU+gzO0lAspGLSOqLQp3uXhmIhFZHmqceqVasq+n7VSrGQskv7GkVPT8/wG8mwFAspq7RDIaWjWEjZZCUU48ePT+29q4liIWWRlVAAtLW1pfr+1UKxkJLLUigAtm3blvYIVUGxkJLKWihAsSgVxUJKJouhkNJRLKQkFIrqp1hIYlkPRWtra9ojVAXFQhLJeiikdJL+iLqMMBte3Mf9z+1i/6FPOe34p+x9fkfmQ9He3q5fB1ACioUE2/DiPu55+mU+PfpnAD4dNZZzFrZh0+amPJlUgk5DJNj9z+06EYo+R48b9z+3K6WJpJJ0ZBHDQw89xKFDh04s911Aa29vP7Fu/vz5zJ8/n5UrV/LRRx8BUF9fz6233kpnZ2e/X3zT1tbG/v376ejoOLGuubmZxsbGfofPs2bNYvHixaxdu7bfL/tdsWIF3d3d/X4/RktLCw0NDf0+Uq6xsZHm5mZWrVp14oerxo8fT1tbG9u2bet3P8Jg/037PvvmoPtj36FPTsyZxblra2sHnVtOjhU+gDtdTU1N3tXVlfYYwVasWDEiz4G/8Q+dHDr65YPRybVj+a/lV6QwkZwMM+t298E+2S6ITkMkSEdHB71bVnPqaOu3fmzNaH501XlDPEuqiU5DZFgdHR3ceeedbN68mT1HJ5z4bkhD7Vh+dNV5LJo7Oe0RpQIUixhG0k0+xaGYM2cOc0BxGKF0GiJDGhgKGdkUixiKr7RXK4VCBlIs5EsUChmMYiH9KBQyFMUihvnz56c9QlkoFPJVFIsYqjEWCoUMR7GIodp+K7dCISEUixj6ftajGigUEkqxGMEUCjkZikUM9fX1aY+QmEIhJ0uxiOHWW29Ne4REFAqJQ7GIofjzF/JGoZC4FIsYij+4Jk8UCklCsRghFApJKnEszGy0mb1oZs9Ey9PNbLuZ7TGzp8zslORjShIKhZRCKY4s7gB2Fi3fBzzo7ucCHwJLS/AemZKn38qtUEipJIqFmU0B/gZ4PFo24ApgfbTJGmBRkvfIov3796c9QhCFQkop6ZHFQ8BdwPFo+WzgkLsfi5b3AoN+rJKZtZpZl5l19fb2Jhyjsoo/hTurFAoptdixMLPrgIPuHutbA+7e7u5N7t5UV1cXdwwZhEIh5ZDkMzgvAb5jZtcCpwFnAg8DtWY2Jjq6mALsSz6mhFIopFxiH1m4+z3uPsXdpwE3A8+7+y3AVuCGaLMlwMbEU2ZMc3Nz2iMMSqGQcirHfRZ3A3ea2R4K1zBWl+E9UtXY2Jj2CF+iUEi5lSQW7r7N3a+Lvn7D3S9y93Pd/UZ3/7wU75ElWfttZAqFVILu4Mw5hUIqRbHIMYVCKkmxiGHWrFlpj6BQSMUpFjEsXrw41fdXKCQNikUMa9euTe29FQpJi2IRw+7du1N5X4VC0qRY5IRCIWlTLHJAoZAsUCxiqORNWQqFZIViEUOlPoNToZAsUSxiqMSneysUkjWKRQYpFJJFikXGKBSSVYpFDC0tLWV5XYVCskyxiKGhoaHkr6lQSNYpFjGsXLmypK+nUEgeKBYpUygkLxSLFCkUkieKRQyl+AxOhULyRrGIIemneysUkkeKRQyrVq2K/VyFQvJKsYihp6cn1vMUCskzxaJCFArJO8UihvHjx5/U9gqFVAPFIoa2trbgbRUKqRaKRQzbtm0L2k6hkGqiWMQQEguFQqqNYlEGCoVUI8WixBQKqVaKRQytra2DrlcopJopFiWiUEi1UyxiaG9v77esUMhIoFgkpFDISKFYJKBQyEgyJu4TzWwq8EtgEuBAu7s/bGYTgaeAacBbwE3u/mHyUTPgpXWw5Sfcy14+/qc1/P7Zj9i8eZtCISNCkiOLY0Cbu88G5gG3mdlsYDmwxd1nAlui5fx7aR10LoPD72I44459wL9cNYo5x19NezKRiogdC3fvcfc/RF9/BOwEJgMLgTXRZmuARQlnzIYtP4Gjn/ZbNerPnxfWi4wAJblmYWbTgLnAdmCSu/d94MMBCqcpgz2n1cy6zKyrt7e3FGOU1+G9J7depMokjoWZnQH8GviBux8pfszdncL1jC9x93Z3b3L3prq6uqRjlN9ZU05uvUiVSRQLM6uhEIon3f3paPV7ZlYfPV4PHEw2YkYs+DHUjO2/rmZsYb3ICBA7FmZmwGpgp7s/UPTQJmBJ9PUSYGP88TLkL2+C5kfgrKmAFf5ufqSwXmQEsMKZQownml0K/B54GTgerf47Ctct1gFfB96m8K3TD77qtZqamryrqyvWHCISxsy63b0p7vNj32fh7v8J2BAPL4j7uiKSTbqDU0SCKBYiEkSxEJEgioWIBFEsRCSIYiEiQRQLEQmiWIhIEMVCRIIoFiISRLEQkSCKhYgEUSxEJIhiISJBFAsRCaJYiEgQxUJEgigWIhJEsRCRIIqFiARRLEQkiGIhIkEUCxEJoliISBDFQkSCKBYiEkSxEJEgioWIBFEsRCSIYiEiQRQLEQmiWIhIEMVCRIIoFiISRLEQkSBliYWZXW1mu8xsj5ktL8d7iEhllTwWZjYaeBS4BpgNtJjZ7FK/j4hUVjmOLC4C9rj7G+7+BfCvwMIyvI+IVNCYMrzmZODdouW9wF8N3MjMWoHWaPFzM9tRhlnK5WvA+2kPEShPs0K+5s3TrADnJXlyOWIRxN3bgXYAM+ty96a0ZjlZeZo3T7NCvubN06xQmDfJ88txGrIPmFq0PCVaJyI5Vo5YvADMNLPpZnYKcDOwqQzvIyIVVPLTEHc/ZmZ/CzwHjAaecPdXhnlae6nnKLM8zZunWSFf8+ZpVkg4r7l7qQYRkSqmOzhFJIhiISJBUo9Flm8NN7OpZrbVzF41s1fM7I5o/UQz22xmr0d/T0h71j5mNtrMXjSzZ6Ll6Wa2Pdq/T0UXnTPBzGrNbL2ZvWZmO83s4qzuWzP7YfRvYIeZdZjZaVnat2b2hJkdLL5faah9aQWPRHO/ZGYXhrxHqrHIwa3hx4A2d58NzANui+ZbDmxx95nAlmg5K+4AdhYt3wc86O7nAh8CS1OZanAPA8+6+/nABRTmzty+NbPJwDKgyd3nULhwfzPZ2re/AK4esG6ofXkNMDP60wo8FvQO7p7aH+Bi4Lmi5XuAe9KcaZh5NwJXAruA+mhdPbAr7dmiWaZE/yiuAJ4BjMIdhmMG298pz3oW8CbRRfai9Znbt/z/XckTKXwH8RngqqztW2AasGO4fQmsAloG2+6r/qR9GjLYreGTU5rlK5nZNGAusB2Y5O490UMHgElpzTXAQ8BdwPFo+WzgkLsfi5aztH+nA73Az6PTpsfNbBwZ3Lfuvg/4KfAO0AMcBrrJ7r7tM9S+jPX/XdqxyAUzOwP4NfADdz9S/JgX0pz695/N7DrgoLt3pz1LoDHAhcBj7j4X+JgBpxwZ2rcTKPww5HSgARjHlw/5M60U+zLtWGT+1nAzq6EQiifd/elo9XtmVh89Xg8cTGu+IpcA3zGztyj8pO8VFK4J1JpZ3813Wdq/e4G97r49Wl5PIR5Z3LffBt509153Pwo8TWF/Z3Xf9hlqX8b6/y7tWGT61nAzM2A1sNPdHyh6aBOwJPp6CYVrGaly93vcfYq7T6OwH59391uArcAN0WaZmBXA3Q8A75pZ309CLgBeJYP7lsLpxzwzOz36N9E3ayb3bZGh9uUm4LvRd0XmAYeLTleGloGLR9cCu4H/Bf4+7XkGzHYphUO3l4A/Rn+upXAtYAvwOvAfwMS0Zx0w93zgmejrGcD/AHuAfwNOTXu+ojm/AXRF+3cDMCGr+xb4R+A1YAfwK+DULO1boIPC9ZSjFI7alg61Lylc+H40+n/uZQrf5Rn2PXS7t4gESfs0RERyQrEQkSCKhYgEUSxEJIhiISJBFAsRCaJYiEiQ/wN9h/Q4IjcwbwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_grid(c1, c1t, c2):\n",
    "\n",
    "    # make position of fish 1\n",
    "    x1 = c1[0]\n",
    "    y1 = c1[1]\n",
    "\n",
    "    # make next positon of fish 1\n",
    "    x1t = c1t[0]\n",
    "    y1t = c1t[1]\n",
    "\n",
    "    # make positon of fish 2\n",
    "    x2 = c2[0]\n",
    "    y2 = c2[1]\n",
    "\n",
    "    # plot positions\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    # markers\n",
    "    ax.plot(x1, y1, \"o\", color = \"tab:blue\")\n",
    "    ax.plot(x1t, y1t, \"o\", color = \"skyblue\")\n",
    "    ax.plot(x2, y2, \"o\", color = \"tab:orange\")\n",
    "\n",
    "    # connectors\n",
    "    ax.plot([x1, x1t], [y1, y1t], color = \"black\", alpha = 1, lw = 1, zorder=-100)\n",
    "    ax.plot([x1, x2], [y1, y2], color = \"black\", alpha = 1, lw = 1, zorder=-100)\n",
    "\n",
    "    # helpers \n",
    "    ax.plot([x1-(x1-x2), x1], [y1, y1], color = \"black\", alpha = 0.5, lw = 1, ls = \"dashed\", zorder=-100)\n",
    "    ax.plot([x2, x2], [y1, y1- (y1-y2)], color = \"black\", alpha = 0.5, lw = 1, ls = \"dashed\", zorder=-100)\n",
    "    ax.plot([x1-(x1-x1t), x1], [y1, y1], color = \"black\", alpha = 0.5, lw = 1, ls = \"dashed\", zorder=-100)\n",
    "    ax.plot([x1t, x1t], [y1, y1- (y1-y1t)], color = \"black\", alpha = 0.5, lw = 1, ls = \"dashed\", zorder=-100)\n",
    "\n",
    "    # limits\n",
    "    ax.set_ylim(0,100)\n",
    "    ax.set_xlim(0,100)\n",
    "    ax.set_aspect(\"equal\")\n",
    "\n",
    "    return ax.lines[3]\n",
    "\n",
    "def add_arrow(line, position=None, direction='right', size=15, color=None):\n",
    "    \"\"\"\n",
    "    add an arrow to a line.\n",
    "\n",
    "    line:       Line2D object\n",
    "    position:   x-position of the arrow. If None, mean of xdata is taken\n",
    "    direction:  'left' or 'right'\n",
    "    size:       size of the arrow in fontsize points\n",
    "    color:      if None, line color is taken.\n",
    "    \"\"\"\n",
    "    if color is None:\n",
    "        color = line.get_color()\n",
    "\n",
    "    xdata = line.get_xdata()\n",
    "    ydata = line.get_ydata()\n",
    "\n",
    "    if position is None:\n",
    "        position = xdata.mean()\n",
    "    # find closest index\n",
    "    start_ind = np.argmin(np.absolute(xdata - position))\n",
    "    if direction == 'right':\n",
    "        end_ind = start_ind + 1\n",
    "    else:\n",
    "        end_ind = start_ind - 1\n",
    "\n",
    "    line.axes.annotate('',\n",
    "        xytext=(xdata[start_ind], ydata[start_ind]),\n",
    "        xy=(xdata[end_ind], ydata[end_ind]),\n",
    "        arrowprops=dict(arrowstyle=\"->\", color=color),\n",
    "        size=size, zorder = -100\n",
    "    )\n",
    "\n",
    "c1 = [50,50]\n",
    "c1t = [70, 70]\n",
    "c2 = [20,20]\n",
    "traj = plot_grid(c1, c1t, c2)\n",
    "add_arrow(traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7853981633974483\n",
      "45.0\n"
     ]
    }
   ],
   "source": [
    "# example values for fish trajectory\n",
    "adj = c1t[0] - c1[0]\n",
    "opp = c1t[1] - c1[1]\n",
    "a = atan2(opp, adj)\n",
    "print(a)\n",
    "print(degrees(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.356194490192345\n",
      "-135.0\n"
     ]
    }
   ],
   "source": [
    "# example values for angle between fish\n",
    "adj  = c2[0] - c1[0]\n",
    "opp = c2[1] - c1[1]\n",
    "b = atan2(opp, adj)\n",
    "print(b)\n",
    "print(degrees(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "radians of trajectory: 0.7853981633974483\n",
      "degrees of trajectory 45.0\n",
      "radians to second fish: 3.9269908169872414\n",
      "degrees to second fish 225.0\n",
      "abs relative angle in radians: 3.141592653589793\n",
      "abs relative angle in degrees: 180.0\n"
     ]
    }
   ],
   "source": [
    "# first transform two radians to keep them positive\n",
    "a = 2*np.pi + a if a < 0 else a\n",
    "b = 2*np.pi + b if b < 0 else b \n",
    "\n",
    "# get the absolute of the relative angle\n",
    "# (because its easier and for now we dont care about left or right)\n",
    "r = np.abs(b - a)\n",
    "\n",
    "print(f\"radians of trajectory: {a}\")\n",
    "print(f\"degrees of trajectory {degrees(a)}\")\n",
    "\n",
    "print(f\"radians to second fish: {b}\")\n",
    "print(f\"degrees to second fish {degrees(b)}\")\n",
    "\n",
    "print(f\"abs relative angle in radians: {r}\")\n",
    "print(f\"abs relative angle in degrees: {degrees(r)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "aim index: 0.0\n"
     ]
    }
   ],
   "source": [
    "# create an index for the \"aim\" of fish 1 towards fish 2\n",
    "\n",
    "# transform r to degrees\n",
    "rd = degrees(r)\n",
    "\n",
    "# norm right side of unit circle to 0\n",
    "if rd <= 180:\n",
    "    aim = 1-rd/180\n",
    "elif rd > 180:\n",
    "    aim = (rd-180)/180\n",
    "\n",
    "print(f\"aim index: {aim}\")\n",
    "\n",
    "# this works: The index is 0 when the relative angle is 180 i.e. fish 1 swims into the opposite direction of fish 2, and 1 if the angle is 0, i.e. fish 1 swims directly towards fish 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'plotstyle'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32m/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb Cell 8\u001b[0m in \u001b[0;36m<cell line: 5>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39m# now we need to incorporate the aim index with the velocity! If we multiply the velocity of a fish at a given time with the aim index, we can get and idea about who approaches whom and how \"vigorously\".\u001b[39;00m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m \n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m \u001b[39m# lets make some functions:\u001b[39;00m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m \u001b[39mfrom\u001b[39;00m \u001b[39mscipy\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39msignal\u001b[39;00m \u001b[39mimport\u001b[39;00m find_peaks\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=4'>5</a>\u001b[0m \u001b[39mfrom\u001b[39;00m \u001b[39mplotstyle\u001b[39;00m \u001b[39mimport\u001b[39;00m PlotStyle\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=5'>6</a>\u001b[0m get_ipython()\u001b[39m.\u001b[39mrun_line_magic(\u001b[39m'\u001b[39m\u001b[39mmatplotlib\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mqt\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/utilities/relative_heading_angle_explorer.ipynb#X10sZmlsZQ%3D%3D?line=7'>8</a>\u001b[0m s \u001b[39m=\u001b[39m PlotStyle()\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'plotstyle'"
     ]
    }
   ],
   "source": [
    "# now we need to incorporate the aim index with the velocity! If we multiply the velocity of a fish at a given time with the aim index, we can get and idea about who approaches whom and how \"vigorously\".\n",
    "\n",
    "# lets make some functions:\n",
    "from scipy.signal import find_peaks\n",
    "from plotstyle import PlotStyle\n",
    "%matplotlib qt\n",
    "\n",
    "s = PlotStyle()\n",
    "\n",
    "# get some nice data:\n",
    "path = \"/home/weygoldt/Data/uni/efish/output/2016-04-20-18_49/\"\n",
    "idx = 8\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "events = pd.read_csv(path + \"events.csv\")\n",
    "grid = gt.GridTracks(path, finespec=False, verbose=False)\n",
    "dyad = gt.Dyad(grid, [int(events.id1[idx]), int(events.id2[idx])])\n",
    "\n",
    "start = gt.utils.find_closest(dyad.times, events.start[idx]-600)\n",
    "stop = gt.utils.find_closest(dyad.times, events.stop[idx]+600)\n",
    "\n",
    "start = 0 \n",
    "stop = len(dyad.times)\n",
    "\n",
    "def find_interactions(dyad, start, stop, maxd, peakprom, plot=False):\n",
    "    \n",
    "    # make arrays for aim index\n",
    "    c1 = [dyad.xpos_smth_id1[start:stop], dyad.ypos_smth_id1[start:stop]]\n",
    "    c2 = [dyad.xpos_smth_id2[start:stop], dyad.ypos_smth_id2[start:stop]]\n",
    "    t = dyad.times[start:stop]\n",
    "\n",
    "    # extract aim index\n",
    "    aims1, relangles1 = gt.utils.aim_index(c1, c2)\n",
    "    aims2, relangles2 = gt.utils.aim_index(c2, c1)\n",
    "\n",
    "    # extract distance index (0 at max dist, 1 at min dist)\n",
    "    dist_index = []\n",
    "    for dpos in dyad.dpos[start:stop]:\n",
    "        # if dpos < maxd:\n",
    "        #     dist_index.append(1-dpos/maxd)\n",
    "        # else:\n",
    "        #     dist_index.append(0)\n",
    "        if dpos < maxd:\n",
    "            dist_index.append(1)\n",
    "        else:\n",
    "            dist_index.append(1)\n",
    "\n",
    "    # extract velocities\n",
    "    v1 = gt.utils.velocity2d(t, c1[0], c1[1])\n",
    "    v2 = gt.utils.velocity2d(t, c2[0], c2[1])\n",
    "\n",
    "    # compute relative velocities\n",
    "    vr = gt.utils.veloc(dyad.times[start:stop], dyad.dpos[start:stop])\n",
    "\n",
    "    interact_index1 = v1 * aims1 * dist_index\n",
    "    interact_index2 = v2 * aims2 * dist_index\n",
    "\n",
    "    # find peaks, i.e. attacks\n",
    "    peaks1 = find_peaks(interact_index1, prominence=peakprom)[0]\n",
    "    peaks2 = find_peaks(interact_index2, prominence=peakprom)[0]\n",
    "\n",
    "    if plot:\n",
    "        fig, ax = plt.subplots(6, 1, figsize=(6,12), sharex= True)\n",
    "\n",
    "        ax[0].set_title(\"Fundamental frequencies\")\n",
    "        ax[0].plot(dyad.times[start:stop], dyad.fund_id1[start:stop])\n",
    "        ax[0].plot(dyad.times[start:stop], dyad.fund_id2[start:stop])\n",
    "\n",
    "        ax[1].set_title(\"Velocities of both fish\")\n",
    "        ax[1].plot(t, v1)\n",
    "        ax[1].plot(t, v2)\n",
    "\n",
    "        ax[2].set_title(\"Aim index of both fish, should have peaks where attack is\")\n",
    "        ax[2].plot(t, aims1)\n",
    "        ax[2].plot(t, aims2)\n",
    "\n",
    "        ax[3].set_title(\"Distance index\")\n",
    "        ax[3].plot(t, dist_index)\n",
    "\n",
    "\n",
    "        ax[4].set_title(\"Attack index, i.e. v * aim_index * dist_index\")\n",
    "        ax[4].plot(t, interact_index1)\n",
    "        ax[4].plot(t, interact_index2)\n",
    "\n",
    "        peakbool1 = np.full(len(interact_index1), False, dtype=bool)\n",
    "        peakbool2 = np.full(len(interact_index2), False, dtype=bool)\n",
    "        peakbool1[peaks1] = True\n",
    "        peakbool2[peaks2] = True\n",
    "\n",
    "        ax[5].set_title(\"Attack index masked for approach phases only\")\n",
    "        ax[5].plot(t, interact_index1)\n",
    "        ax[5].plot(t, interact_index2)\n",
    "        ax[5].plot(t[peakbool1], interact_index1[peakbool1], \".\")\n",
    "        ax[5].plot(t[peakbool2], interact_index2[peakbool2], \".\")\n",
    "\n",
    "    return peaks1, peaks2\n",
    "\n",
    "peaks1, peaks2 = gt.utils.find_interactions(dyad, start, stop, maxd=400, peakprom=10, plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 350.0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(dyad.xpos_smth_id1[peaks1[0]-30:peaks1[0]+2], dyad.ypos_smth_id1[peaks1[0]-30:peaks1[0]+2])\n",
    "ax.plot(dyad.xpos_smth_id2[peaks1[0]-30:peaks1[0]+2], dyad.ypos_smth_id2[peaks1[0]-30:peaks1[0]+2])\n",
    "ax.set_xlim(0,350)\n",
    "ax.set_ylim(0,350)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'relangles1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m/home/weygoldt/Data/uni/efish/code/trigtest.ipynb Cell 9'\u001b[0m in \u001b[0;36m<cell line: 4>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/trigtest.ipynb#ch0000008?line=0'>1</a>\u001b[0m \u001b[39m#start = gt.utils.find_closest(dyad.times, events.start[idx])\u001b[39;00m\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/trigtest.ipynb#ch0000008?line=1'>2</a>\u001b[0m \u001b[39m#stop = gt.utils.find_closest(dyad.times, events.stop[idx])\u001b[39;00m\n\u001b[0;32m----> <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/trigtest.ipynb#ch0000008?line=3'>4</a>\u001b[0m eventra1 \u001b[39m=\u001b[39m relangles1[start:stop]\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/trigtest.ipynb#ch0000008?line=4'>5</a>\u001b[0m eventra2 \u001b[39m=\u001b[39m relangles2[start:stop]\n\u001b[1;32m      <a href='vscode-notebook-cell:/home/weygoldt/Data/uni/efish/code/trigtest.ipynb#ch0000008?line=6'>7</a>\u001b[0m distribution1 \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mhistogram(eventra1,bins\u001b[39m=\u001b[39m\u001b[39m120\u001b[39m, \u001b[39mrange\u001b[39m\u001b[39m=\u001b[39m(\u001b[39m0\u001b[39m,\u001b[39m360\u001b[39m))[\u001b[39m0\u001b[39m]\n",
      "\u001b[0;31mNameError\u001b[0m: name 'relangles1' is not defined"
     ]
    }
   ],
   "source": [
    "#start = gt.utils.find_closest(dyad.times, events.start[idx])\n",
    "#stop = gt.utils.find_closest(dyad.times, events.stop[idx])\n",
    "\n",
    "eventra1 = relangles1[start:stop]\n",
    "eventra2 = relangles2[start:stop]\n",
    "\n",
    "distribution1 = np.histogram(eventra1,bins=120, range=(0,360))[0]\n",
    "distribution2 = np.histogram(eventra2,bins=120, range=(0,360))[0]\n",
    "theta = np.linspace(0,360,120, endpoint=False)\n",
    "\n",
    "f = plt.figure()\n",
    "ax = f.add_subplot(polar=True)\n",
    "ax.set_theta_zero_location('N'); ax.set_theta_direction(-1)\n",
    "ax.set_ylim(-20,21)\n",
    "ax.bar(theta, distribution1, edgecolor=\"black\", width = 20/360 , alpha=0.5)\n",
    "ax.bar(theta, distribution2, edgecolor=\"black\", width = 20/360 , alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fed5ee9d3c0>]"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# extract distance index (0 at max dist, 1 at min dist, thresh at 50 cm dist)\n",
    "mind, maxd = 0, 300\n",
    "dist_index = []\n",
    "for dpos in dyad.dpos:\n",
    "    if dpos < maxd:\n",
    "        dist_index.append(1-dpos/maxd)\n",
    "    else:\n",
    "        dist_index.append(0)\n",
    "\n",
    "\n",
    "# dist_index = [(1-dpos/maxd) if dpos < maxd else 0 for dpos in dyad.dpos]\n",
    "plt.plot(dist_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "119033"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(dist_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2063, 2196])"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "peaks1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.32769999999999966"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dyad.times[1]-dyad.times[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.5 ('env': venv)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "fc22036bd5413fbeea744e8aabe54594d1863936e7bd1686ce73b6aa81314412"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}