{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Trajectory\n",
    "``noctiluca`` is intended to facilitate handling and analysis of single (or multi-) particle tracking trajectories. The centerpiece is the [Trajectory](../noctiluca.rst#noctiluca.trajectory.Trajectory) class; let's see how to use it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import noctiluca as nl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "At its very core, a `Trajectory` object is similar to a numpy array. As such, it can be initialized simply from raw data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory([1, 0.5, np.nan, 5.4, -2])\n",
    "nl.plot.vstime(traj)\n",
    "\n",
    "plt.title('Example trajectory')\n",
    "plt.xlabel('time [frames]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we specified a single list of data points. What if we have a 2D trajectory?\n",
    "\n",
    "Usual SPT trajectories will consist of multiple spatial (xyz) coordinates for a particle. In some cases there might also be multiple particles tracked simultaneously (e.g. multi-color experiments). So in its most general form, a `Trajectory` has\n",
    "\n",
    "+ $N$ particles,\n",
    "+ $T$ frames,\n",
    "+ $d$ spatial dimensions,\n",
    "\n",
    "which is represented internally as a numpy array of shape ``(N, T, d)``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(traj.data.shape)\n",
    "print(f\"N = {traj.N}, T = {traj.T}, d = {traj.d}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, `Trajectory()` accepts arrays of shapes ``(T,)``, ``(T, d)``, or ``(N, T, d)``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj0 = nl.Trajectory(np.random.normal(size=(10,)))\n",
    "traj1 = nl.Trajectory(np.random.normal(size=(10, 3)))\n",
    "traj2 = nl.Trajectory(np.random.normal(size=(2, 10, 2),\n",
    "                                       loc=[[[0]], [[5]]], # <-- offset to distinguish loci\n",
    "                                      ))\n",
    "\n",
    "########## Plotting ############\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=[15, 4])\n",
    "\n",
    "for i, traj in enumerate([traj0, traj1, traj2]):\n",
    "    ax = axs[i]\n",
    "    nl.plot.vstime(traj, ax)\n",
    "    ax.set_title(f\"traj{i}: N = {traj.N}, T = {traj.T}, d = {traj.d}\")\n",
    "    ax.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metadata\n",
    "Trajectories often carry some annotation, metadata, or other additional information that might be useful. These are stored in the `Trajectory.meta` dict:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory([1, 2, 3], comments=\"This is an example trajectory; not for consumption\")\n",
    "traj.meta['date'] = \"2306/07/13\"\n",
    "traj.meta['spot intensity'] = np.array([1.02, 1.5, 0.7])\n",
    "traj.meta['localization_error'] = np.array([0.5])\n",
    "\n",
    "print(\"traj.meta = {\")\n",
    "for key, val in traj.meta.items():\n",
    "    print(f\"    {repr(key):>20s} : {repr(val)}\")\n",
    "print(\"}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entries in this ``meta``-dict can be pretty much anything, whatever you deem useful to know about each trajectory. Note that instead of writing to ``traj.meta`` explicitly, you can also provide metadata as keyword arguments to the ``Trajectory(...)`` constructor.\n",
    "\n",
    "Beyond the ``meta``-dict, a ``Trajectory`` has two (as of 2022/09/20) additional attributes: ``localization_error`` and ``parity``. These are hardcoded to allow propagating them through modifiers; see below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying trajectories\n",
    "There are a few common operations that one might want to apply to a trajectory, such as taking an absolute value or calculating an increment trajectory. These are implemented as so-called modifiers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory([[1.3, -2], [8.9, 3], [-0.5, 4], [3.2, 1.7]])\n",
    "abstraj = traj.abs() # takes Euclidean norm --> d = 1\n",
    "print(f\"abstraj.d = {abstraj.d}\")\n",
    "difftraj = traj.diff() # returns increments (like numpy.diff())\n",
    "traj2 = traj.rescale(1e3) # useful for changing units (e.g. nm <--> μm)\n",
    "\n",
    "traj = nl.Trajectory(np.random.normal(size=(2, 10, 3)))\n",
    "reltraj = traj.relative() # returns trajectory of distances between particle 1 and 2\n",
    "print(f\"reltraj.N = {reltraj.N}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A modifier always returns a new trajectory with the given modification applied. By default, metadata is not copied (since it might not apply anymore). To override this, specify the ``keepmeta`` keyword argument for any of the above modifiers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory([1, -2, 3], mymeta=5)\n",
    "abstraj = traj.abs()\n",
    "print(abstraj.meta) # {}\n",
    "\n",
    "abstraj = traj.abs(keepmeta=['mymeta'])\n",
    "print(abstraj.meta) # {'mymeta': 5}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Missing data\n",
    "Real data is not perfect; specifically, often there are missing data points. These are specified by ``np.nan`` and analysis methods should generally be able to deal with this. It is often useful to know how many frames of a given trajectory are actually valid, i.e. not ``np.nan``; use ``Trajectory.count_valid_frames()`` or its alias ``Trajectory.F`` (F for \"(valid) Frame\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory([1, 2, np.nan, 3, np.nan, np.nan, 5])\n",
    "print(f\"valid frames in the trajectory: {traj.count_valid_frames()}\")\n",
    "print(f\"                       traj.F = {traj.F}\")\n",
    "print(f\"     for reference: len(traj) = {len(traj)} = traj.T\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next: many trajectories\n",
    "Check out the [next example](02_TaggedSet.ipynb) for ``TaggedSet`` to learn how ``noctiluca`` represents complete datasets containing many individual trajectories. Alternatively, see below for a few more tips and tricks on ``Trajectory`` handling. For full reference, always refer to the [documentation](../noctiluca.rst#API-reference)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing\n",
    "Trajectories can be indexed along the time dimension using the ``[]`` operator, which returns a numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory(np.random.normal(size=(2, 10, 3)))\n",
    "print(f\"traj[:].shape  = {traj[:].shape}\")\n",
    "print(f\"traj[:7].shape = {traj[:7].shape}\")\n",
    "print(type(traj[:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The three dimensions of the internal array ``(N, T, d)`` may be present or absent in the output, following these rules:\n",
    "\n",
    "+ ``N`` is squeezed, i.e. removed if ``N == 1``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a trajectory with N = 1\n",
    "print(f\"traj.relative()[:].shape = {traj.relative()[:].shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ ``T`` follows numpy convention: whether it's present or not depends on the format of the key:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"traj[:5].shape  = {traj[:5].shape}\")\n",
    "print(f\"traj[5].shape   = {traj[5].shape}\")\n",
    "print(f\"traj[[5]].shape = {traj[[5]].shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ `d` is always present in the output, even if there is only one spatial dimension. This is such that analysis code can be fully agnostic of spatial dimension; squeezing this last dimension can be done by hand if necessary, as illustrated below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abstraj = traj.abs() # get a trajectory with d = 1\n",
    "print(f\"abstraj[:].shape                  = {abstraj[:].shape}\")\n",
    "print(f\"abstraj[:][..., 0].shape          = {abstraj[:][..., 0].shape}\")\n",
    "print(f\"abstraj.relative()[:].shape       = {abstraj.relative()[:].shape}\")\n",
    "print(f\"abstraj.relative()[:][:, 0].shape = {abstraj.relative()[:][:, 0].shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "We already saw the ``nl.plot.vstime()`` method for plotting a trajectory over time. If you have a complicated trajectory (with multiple particles and/or spatial dimensions), you can use modifiers to simplify plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory(np.random.normal(size=(2, 10, 3)))\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=[15, 4])\n",
    "\n",
    "ax = axs[0]\n",
    "nl.plot.vstime(traj, ax=ax)\n",
    "ax.set_title('Complicated plot; difficult to understand')\n",
    "\n",
    "ax = axs[1]\n",
    "nl.plot.vstime(traj.relative(), ax=ax)\n",
    "ax.set_title('Plotting only difference between particles')\n",
    "\n",
    "ax = axs[2]\n",
    "nl.plot.vstime(traj.relative().abs(), ax=ax)\n",
    "ax.set_title('Absolute distance between particles')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there is no useful simplification to be made to the data before plotting, you can beautify the plot after the fact to make it more digestible. To that end, ``nl.plot.vstime()`` returns a list of the lines it adds to the axes, so we can simply adjust their properties. The returned list is sorted by ``N * d``, i.e. ``[:d]`` are the spatial dimensions of particle 0, ``[d:(2*d)]`` are particle 1, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory(np.random.normal(size=(2, 10, 3)))\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=[10, 4])\n",
    "\n",
    "ax = axs[0]\n",
    "nl.plot.vstime(traj, ax=ax)\n",
    "ax.set_title('Complicated plot; difficult to understand')\n",
    "\n",
    "ax = axs[1]\n",
    "lines = nl.plot.vstime(traj, ax=ax)\n",
    "for i in range(traj.d):\n",
    "    lines[traj.d+i].set_color(lines[i].get_color())\n",
    "    lines[traj.d+i].set_linestyle('--')\n",
    "    \n",
    "    lines[i].set_label(f'd={i}')\n",
    "    lines[traj.d+i].set_label('')\n",
    "\n",
    "# Add some fake legend entries\n",
    "ax.plot(np.nan, np.nan, linestyle='-',  color='k', label='N=0')\n",
    "ax.plot(np.nan, np.nan, linestyle='--', color='k', label='N=1')\n",
    "\n",
    "ax.set_title('(somewhat) more parseable plot')\n",
    "ax.legend(loc=(1.02, 0.5))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is another method for visualizing trajectories: ``nl.plot.spatial()`` plots a given trajectory in the xy-coordinate plane. For trajectories with more than two spatial dimensions you can specify which dimensions you want to plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = nl.Trajectory(np.random.normal(size=(2, 10, 3)))\n",
    "\n",
    "fig, axs = plt.subplots(1, 3,\n",
    "                        figsize=[15, 4],\n",
    "                        gridspec_kw={'wspace' : 0.3},\n",
    "                       )\n",
    "\n",
    "for i, (ax, dims) in enumerate(zip(axs, [(0, 1), (1, 2), (2, 0)])):\n",
    "    lines = nl.plot.spatial(traj, ax=ax, dims=dims)\n",
    "    lines[1].set_color(lines[0].get_color())\n",
    "    lines[1].set_linestyle('--')\n",
    "    \n",
    "    dims_chr = tuple(chr(ord('x')+d) for d in dims)\n",
    "    ax.set_title(dims_chr[0]+dims_chr[1]+'-plane')\n",
    "    ax.set_xlabel(dims_chr[0])\n",
    "    ax.set_ylabel(dims_chr[1])\n",
    "    \n",
    "axs[-1].legend(loc=(1.02, 0.5))\n",
    "plt.show()"
   ]
  }
 ],
 "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
}