[demo]detect_duplicates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1ea3ca61",
   "metadata": {},
   "source": [
    "Adapted from methods in dynamic-vision pipeline.\n",
    "\n",
    "https://github.com/cajal/dynamic-vision/blob/7ab1ae8530a7681645b4f56766b6359cbb6fc0b3/scans/scan.py#L529"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "a62a5b63",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "from utils.signal import CenteredSpline\n",
    "\n",
    "import datajoint as dj\n",
    "from pipeline import meso\n",
    "from stimulus import stimulus\n",
    "stats = dj.create_virtual_module('neurodata_static_stats', 'neurodata_static_stats')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "58cb8925",
   "metadata": {},
   "outputs": [],
   "source": [
    "scan_key  = {'animal_id': 21067, 'session': 10, 'scan_idx': 18} \n",
    "max_distance = 80 # microns\n",
    "min_correlation = 0.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edbcd19b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# restrict to masks classified as soma\n",
    "soma = (meso.ScanSet.Unit & scan_key & (meso.MaskClassification.Type() & dict(type=\"soma\"))).proj()\n",
    "\n",
    "# get the activity traces, unit positions, field membership, ms_delays\n",
    "fetch_strs = ['KEY','um_x','um_y','um_z','ms_delay','trace','field','pearson']\n",
    "joint_rel = meso.ScanSet.UnitInfo * meso.Activity.Trace * soma * stats.Oracle.UnitScores\n",
    "unit_keys,xs,ys,zs,ms_delays,traces,fields,oracles = joint_rel.fetch(*fetch_strs)\n",
    "delays = ms_delays / 1000 # correct from ms to s\n",
    "\n",
    "positions = np.stack([xs, ys, zs], 1)\n",
    "\n",
    "df = pd.DataFrame(unit_keys)\n",
    "df['oracle'] = oracles\n",
    "df['position'] = list(map(np.array,positions))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e580e78a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# downsample from scan depth times (nframes * ndepths) to scan frame times \n",
    "ndepth = len(dj.U('z') & (meso.ScanInfo.Field & scan_key))\n",
    "frame_times = (stimulus.Sync & scan_key).fetch1(\"frame_times\").squeeze()[::ndepth]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc14d3f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create hanning filter for desired downsample frequency\n",
    "unique_hz = 1\n",
    "source_hz = 1/ np.nanmedian(np.diff(frame_times))\n",
    "\n",
    "M = 2*round(source_hz/unique_hz)+1\n",
    "n = np.arange(M)    \n",
    "h = 0.5 - 0.5 * np.cos(2.0 * np.pi * n / (M - 1))\n",
    "h /= h.sum()\n",
    "\n",
    "#  apply low pass filter to traces\n",
    "traces = [np.convolve(t,h,mode='same') for t in traces]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "23fa23ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clip scan times and spikes to the same length\n",
    "assert abs(len(frame_times) - len(traces[0])) < 2, \"length of times differs from spikes by more than 1\"\n",
    "min_length = min(len(frame_times), len(traces[0]))\n",
    "frame_times = frame_times[:min_length]\n",
    "traces = [t[:min_length] for t in traces]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "651231de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample at specified hz\n",
    "scan_start = frame_times[0] + delays.max()\n",
    "n_samples = math.floor((frame_times[-1] - scan_start) * unique_hz)\n",
    "sample_times = scan_start + np.arange(n_samples) / unique_hz\n",
    "\n",
    "traces = np.stack(\n",
    "    [\n",
    "        CenteredSpline(frame_times + d, f)(sample_times)\n",
    "        for d, f in zip(delays, traces)\n",
    "    ],\n",
    "    0,\n",
    ")\n",
    "norm_traces = (traces - traces.mean(-1, keepdims=True)) / (traces.std(-1, keepdims=True) + 1e-8)\n",
    "df['norm_trace'] = list(map(np.array,norm_traces))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d5a88183",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Undirected graph\n",
    "G = nx.Graph()\n",
    "G.add_nodes_from(range(len(traces)))\n",
    "for node, (norm_trace, position) in enumerate(zip(norm_traces, positions)):\n",
    "    # distance threshold\n",
    "    distances = np.linalg.norm(position[None] - positions, 2, 1)\n",
    "    distance_idx = np.nonzero(distances < max_distance)[0]\n",
    "\n",
    "    # correlation threshold\n",
    "    corrs = (norm_trace[None] * norm_traces[distance_idx]).mean(-1)\n",
    "    corr_idx = np.nonzero(corrs > min_correlation)[0]\n",
    "\n",
    "    # add edges\n",
    "    G.add_edges_from([(node, n) for n in distance_idx[corr_idx].tolist()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "b1aa8858",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select representative unit from that group using max oracle\n",
    "neuron_keys = []\n",
    "rep_units = []\n",
    "unit_keys = []\n",
    "distances = []\n",
    "corrs = []\n",
    "for neuron_id, node_set in enumerate(nx.connected_components(G)):\n",
    "    _df = df.loc[list(node_set)]\n",
    "    unit_id = _df.loc[_df[\"oracle\"].idxmax(), \"unit_id\"]\n",
    "    neuron_keys += [dict(neuron_id=neuron_id, unit_id=unit_id)]\n",
    "    rep_units += [unit_id]\n",
    "    unit_keys += [dict(neuron_id=neuron_id, unit_id=uid) for uid in _df.unit_id]\n",
    "    if len(_df)>1:\n",
    "        com = np.mean(np.vstack(_df['position']),axis=0)\n",
    "        distances.append(np.linalg.norm(com - np.vstack(_df['position']),2,1))\n",
    "        corrs.append(np.vstack((np.array(_df.loc[_df['unit_id']==unit_id,'norm_trace']) * \n",
    "                                np.array(_df.loc[_df['unit_id']!=unit_id,'norm_trace']))).mean(-1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "65fb2372",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.00e+00, 2.75e+03, 6.10e+02, 5.00e+02, 3.43e+02, 2.49e+02,\n",
       "        2.00e+01, 4.00e+00, 3.00e+00, 2.00e+00]),\n",
       " array([-0.5,  0.6,  1.7,  2.8,  3.9,  5. ,  6.1,  7.2,  8.3,  9.4, 10.5]),\n",
       " <BarContainer object of 10 artists>)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAPJ0lEQVR4nO3dXawcd3nH8e+vNtDyosaRXSvYVh1VbitTiSRahbSpqrSUxEmrOtxEQSpYUSRzkbShQqoMN0FwwwUvLRKNZIiLUdNEEQTFQhHBdZG4IngdoiSOiXIUEmzXiQ81DbRI0MDTix03i7F93vbs2Of//UirnX1mducZ7dFv5/xndjZVhSSpDb/WdwOSpOkx9CWpIYa+JDXE0Jekhhj6ktSQ1X03cD5r166tzZs3992GJF1UDh069IOqWne2eRd06G/evJnhcNh3G5J0UUny4rnmObwjSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNuaC/kXuxStLbuv1RHEnn456+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWrInKGfZFOSbyR5JsnhJHd19Y8kOZ7kie5209hzPpRkJsmzSW4Yq2/rajNJdi3PJkmSzmU+19N/FfhgVT2e5C3AoST7u3mfrqpPjC+cZCtwK/A24K3AvyX53W72Z4F3AceAg0n2VdUzk9gQSdLc5gz9qjoBnOimf5zkCLDhPE/ZDjxQVT8FvpdkBri6mzdTVc8DJHmgW9bQl6QpWdCYfpLNwJXAY13pziRPJtmTZE1X2wAcHXvasa52rrokaUrmHfpJ3gx8GfhAVf0IuAf4HeAKRv8JfHISDSXZmWSYZDg7OzuJl5QkdeYV+klexyjw76uqhwCq6uWq+nlV/QL4HK8N4RwHNo09fWNXO1f9l1TV7qoaVNVg3bp1C90eSdJ5zOfsnQD3Akeq6lNj9cvGFns38HQ3vQ+4NckbklwObAG+DRwEtiS5PMnrGR3s3TeZzZAkzcd8zt65Fngv8FSSJ7rah4H3JLkCKOAF4P0AVXU4yYOMDtC+CtxRVT8HSHIn8CiwCthTVYcntiWSpDmlqvru4ZwGg0ENh8O+21iw0T9H/biQ309J05HkUFUNzjbPb+RKUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDVkztBPsinJN5I8k+Rwkru6+qVJ9id5rrtf09WT5DNJZpI8meSqsdfa0S3/XJIdy7dZkqSzmc+e/qvAB6tqK3ANcEeSrcAu4EBVbQEOdI8BbgS2dLedwD0w+pAA7gbeAVwN3H36g0KSNB1zhn5Vnaiqx7vpHwNHgA3AdmBvt9he4OZuejvwxRr5FnBJksuAG4D9VXWqqn4I7Ae2TXJjJEnnt6Ax/SSbgSuBx4D1VXWim/USsL6b3gAcHXvasa52rvqZ69iZZJhkODs7u5D2JElzmHfoJ3kz8GXgA1X1o/F5VVVATaKhqtpdVYOqGqxbt24SLylJ6swr9JO8jlHg31dVD3Xll7thG7r7k139OLBp7Okbu9q56pKkKZnP2TsB7gWOVNWnxmbtA06fgbMDeHis/r7uLJ5rgFe6YaBHgeuTrOkO4F7f1SRJU7J6HstcC7wXeCrJE13tw8DHgQeT3A68CNzSzXsEuAmYAX4C3AZQVaeSfAw42C330ao6NYmNkCTNT0bD8RemwWBQw+Gw7zYWbPTPUT8u5PdT0nQkOVRVg7PN8xu5ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNmTP0k+xJcjLJ02O1jyQ5nuSJ7nbT2LwPJZlJ8mySG8bq27raTJJdk98USdJc5rOn/wVg21nqn66qK7rbIwBJtgK3Am/rnvNPSVYlWQV8FrgR2Aq8p1tWkjRFq+daoKq+mWTzPF9vO/BAVf0U+F6SGeDqbt5MVT0PkOSBbtlnFt6yJGmxljKmf2eSJ7vhnzVdbQNwdGyZY13tXPVfkWRnkmGS4ezs7BLakySdabGhfw/wO8AVwAngk5NqqKp2V9Wgqgbr1q2b1MtKkpjH8M7ZVNXLp6eTfA74avfwOLBpbNGNXY3z1CVJU7KoPf0kl409fDdw+syefcCtSd6Q5HJgC/Bt4CCwJcnlSV7P6GDvvsW3LUlajDn39JPcD1wHrE1yDLgbuC7JFUABLwDvB6iqw0keZHSA9lXgjqr6efc6dwKPAquAPVV1eNIbI0k6v1RV3z2c02AwqOFw2HcbC5akt3VfyO+npOlIcqiqBmeb5zdyJakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1JA5Qz/JniQnkzw9Vrs0yf4kz3X3a7p6knwmyUySJ5NcNfacHd3yzyXZsTybI0k6n/ns6X8B2HZGbRdwoKq2AAe6xwA3Alu6207gHhh9SAB3A+8ArgbuPv1BIUmanjlDv6q+CZw6o7wd2NtN7wVuHqt/sUa+BVyS5DLgBmB/VZ2qqh8C+/nVDxJJ0jJb7Jj++qo60U2/BKzvpjcAR8eWO9bVzlX/FUl2JhkmGc7Ozi6yPUnS2Sz5QG5VFVAT6OX06+2uqkFVDdatWzepl5UksfjQf7kbtqG7P9nVjwObxpbb2NXOVZckTdFiQ38fcPoMnB3Aw2P193Vn8VwDvNINAz0KXJ9kTXcA9/quJkmaotVzLZDkfuA6YG2SY4zOwvk48GCS24EXgVu6xR8BbgJmgJ8AtwFU1akkHwMOdst9tKrOPDgsSVpmGQ3JX5gGg0ENh8O+21iwJL2t+0J+PyVNR5JDVTU42zy/kStJDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWrI6r4b0GQl6WW9VdXLeiUtjHv6ktQQQ1+SGmLoS1JDDH1JasiSQj/JC0meSvJEkmFXuzTJ/iTPdfdrunqSfCbJTJInk1w1iQ2QJM3fJPb0/7SqrqiqQfd4F3CgqrYAB7rHADcCW7rbTuCeCaxbkrQAyzG8sx3Y203vBW4eq3+xRr4FXJLksmVYvyTpHJYa+gV8PcmhJDu72vqqOtFNvwSs76Y3AEfHnnusq/2SJDuTDJMMZ2dnl9iepiVJLzdJC7PUL2f9cVUdT/JbwP4k3x2fWVWVZEHf2qmq3cBugMFg4Dd+JGmClrSnX1XHu/uTwFeAq4GXTw/bdPcnu8WPA5vGnr6xq0mSpmTRoZ/kTUnecnoauB54GtgH7OgW2wE83E3vA97XncVzDfDK2DCQJGkKljK8sx74Sjeuuhr416r6WpKDwINJbgdeBG7pln8EuAmYAX4C3LaEdUuSFmHRoV9VzwNvP0v9P4F3nqVewB2LXZ8kaen8Rq4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiL+Rq4tan9ff8XeBdTFyT1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xFM2pUXq63RRTxXVUrinL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkOmHvpJtiV5NslMkl3TXr8ktWyqP6KSZBXwWeBdwDHgYJJ9VfXMNPuQLmZ9/XgL+AMuK8G0fznramCmqp4HSPIAsB0w9KWLgL8WdvGbduhvAI6OPT4GvGN8gSQ7gZ3dw/9O8uyUehu3FvhBD+udBrft4rWSt++829bnfzcT0Mf79tvnmnHB/UZuVe0GdvfZQ5JhVQ367GG5uG0Xr5W8fW7b9Ez7QO5xYNPY441dTZI0BdMO/YPAliSXJ3k9cCuwb8o9SFKzpjq8U1WvJrkTeBRYBeypqsPT7GGeeh1eWmZu28VrJW+f2zYl8ai4JLXDb+RKUkMMfUlqiKE/ZqVeIiLJpiTfSPJMksNJ7uq7p0lLsirJd5J8te9eJi3JJUm+lOS7SY4k+cO+e5qUJH/X/U0+neT+JL/ed09LkWRPkpNJnh6rXZpkf5Lnuvs1ffZo6HfGLhFxI7AVeE+Srf12NTGvAh+sqq3ANcAdK2jbTrsLONJ3E8vkH4GvVdXvA29nhWxnkg3A3wKDqvoDRid33NpvV0v2BWDbGbVdwIGq2gIc6B73xtB/zf9fIqKqfgacvkTERa+qTlTV4930jxmFxoZ+u5qcJBuBvwA+33cvk5bkN4E/Ae4FqKqfVdV/9drUZK0GfiPJauCNwH/03M+SVNU3gVNnlLcDe7vpvcDN0+zpTIb+a852iYgVE4ynJdkMXAk81nMrk/QPwN8Dv+i5j+VwOTAL/HM3fPX5JG/qu6lJqKrjwCeA7wMngFeq6uv9drUs1lfViW76JWB9n80Y+g1J8mbgy8AHqupHffczCUn+EjhZVYf67mWZrAauAu6pqiuB/6Hn4YFJ6ca2tzP6YHsr8KYkf91vV8urRufI93qevKH/mhV9iYgkr2MU+PdV1UN99zNB1wJ/leQFRkNyf5bkX/ptaaKOAceq6vR/Zl9i9CGwEvw58L2qmq2q/wUeAv6o556Ww8tJLgPo7k/22Yyh/5oVe4mIjC5ReC9wpKo+1Xc/k1RVH6qqjVW1mdF79u9VtWL2FqvqJeBokt/rSu9k5VyK/PvANUne2P2NvpMVcpD6DPuAHd30DuDhHnu58K6y2ZeL6BIRi3Et8F7gqSRPdLUPV9Uj/bWkBfgb4L5uZ+R54Lae+5mIqnosyZeAxxmdYfYdLrBLFixUkvuB64C1SY4BdwMfBx5McjvwInBLfx16GQZJaorDO5LUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNeT/ANIY2HhFvE8hAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist([len(x) for x in nx.connected_components(G)],bins = np.linspace(-0.5,10.5,11),color='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "id": "42f90d0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([6.060e+02, 9.600e+01, 1.985e+03, 5.300e+01, 0.000e+00, 1.412e+03,\n",
       "        1.000e+01, 7.860e+02, 4.400e+01, 2.000e+00, 4.070e+02, 5.000e+00,\n",
       "        1.270e+02, 1.900e+01, 0.000e+00, 4.800e+01, 0.000e+00, 8.000e+00,\n",
       "        4.000e+00, 0.000e+00, 6.000e+00, 0.000e+00, 3.000e+00, 1.000e+00,\n",
       "        0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00]),\n",
       " array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,\n",
       "        13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,\n",
       "        26., 27., 28., 29., 30.]),\n",
       " <BarContainer object of 30 artists>)"
      ]
     },
     "execution_count": 240,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAATfUlEQVR4nO3df6xk5X3f8fenYLsRtsVSbtFmd+lia50KR+lijzBVbIvWDb8UBVxVFKSajWtlbQUkW47UYvcPqCMkKzV2hZquhWtkkGwIjU1ZVaT2BlmhkQrmLt7w04QFg9jVevemJMbUES3w7R/z3Hiyvrs7987s/fW8X9JoznzPmTnPc8/uZ46ec+acVBWSpD78nZVugCRp+Rj6ktQRQ1+SOmLoS1JHDH1J6sipK92AEznzzDNr69atK90MSVoz9u7d+xdVNbPQvFUf+lu3bmV2dnalmyFJa0aSF441z+EdSeqIoS9JHTH0Jakjhr4kdeSEoZ9kS5LvJnkyyRNJPtnqZyTZk+SZ9ryh1ZPkliT7kzya5D0jn7WjLf9Mkh0nr1uSpIWMs6f/GvA7VXUucAFwbZJzgeuB+6tqG3B/ew1wKbCtPXYCu2D4JQHcALwPOB+4Yf6LQpK0PE4Y+lV1qKoeadM/AZ4CNgGXA7e3xW4HrmjTlwN31NCDwOlJNgIXA3uq6qWq+ktgD3DJNDsjSTq+RY3pJ9kKnAc8BJxVVYfarB8BZ7XpTcCLI2870GrHqi+0np1JZpPMzs3NLaaJkqTjGDv0k7wV+Cbwqap6eXReDS/KP7UL81fVrVU1qKrBzMyCPyqTJC3BWKGf5E0MA//rVfWtVj7chm1oz0da/SCwZeTtm1vtWPU1I8lYD0larcY5eyfAV4GnquqLI7N2A/Nn4OwA7h2pX9PO4rkA+HEbBvo2cFGSDe0A7kWtJklaJuNce+dXgY8AjyXZ12qfBT4P3J3kY8ALwJVt3n3AZcB+4KfARwGq6qUkvws83Jb7XFW9NI1OSJLGk9V+j9zBYFCr5YJr4w7drPa/qaT1LcneqhosNM9f5EpSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHxrlH7m1JjiR5fKT2B0n2tcfz87dRTLI1yV+PzPvyyHvem+SxJPuT3BLvIC5Jy26ce+R+DfhPwB3zhar6l/PTSW4Gfjyy/LNVtX2Bz9kF/BbwEMP76F4C/NGiWyxJWrIT7ulX1QPAgjcwb3vrVwJ3Hu8zkmwE3l5VD9bwBrJ3AFcsurWSpIlMOqb/AeBwVT0zUjsnyfeT/EmSD7TaJuDAyDIHWm1BSXYmmU0yOzc3N2ETJUnzJg39q/nbe/mHgLOr6jzg08A3krx9sR9aVbdW1aCqBjMzMxM2UZI0b5wx/QUlORX458B752tV9Srwapvem+RZ4F3AQWDzyNs3t5okaRlNsqf/z4AfVNXfDNskmUlySpt+B7ANeK6qDgEvJ7mgHQe4Brh3gnVLkpZgnFM27wT+F/BLSQ4k+VibdRU/fwD3g8Cj7RTOPwQ+UVXzB4F/G/gvwH7gWTxzR5KWXYYn06xeg8GgZmdnV7oZAIz704LV/jeVtL4l2VtVg4Xm+YtcSeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6sg4t0u8LcmRJI+P1G5McjDJvva4bGTeZ5LsT/J0kotH6pe02v4k10+/K5KkExlnT/9rwCUL1L9UVdvb4z6AJOcyvHfuu9t7/nOSU9rN0n8fuBQ4F7i6LStJWkannmiBqnogydYxP+9y4K6qehX4YZL9wPlt3v6qeg4gyV1t2ScX32RJ0lJNMqZ/XZJH2/DPhlbbBLw4ssyBVjtWfUFJdiaZTTI7Nzc3QRMlSaOWGvq7gHcC24FDwM3TahBAVd1aVYOqGszMzEzzoyWpaycc3llIVR2en07yFeC/t5cHgS0ji25uNY5TlyQtkyXt6SfZOPLyw8D8mT27gauSvCXJOcA24HvAw8C2JOckeTPDg727l95sSdJSnHBPP8mdwIXAmUkOADcAFybZDhTwPPBxgKp6IsndDA/QvgZcW1Wvt8+5Dvg2cApwW1U9Me3OSJKOL1W10m04rsFgULOzsyvdDACSjLXcav+bSlrfkuytqsFC8/xFriR1xNCXpI4Y+pLUkSWdsqnl5/EESdPgnr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHThj6SW5LciTJ4yO1/5DkB0keTXJPktNbfWuSv06yrz2+PPKe9yZ5LMn+JLdk3MtGSpKmZpw9/a8BlxxV2wP8clX9CvDnwGdG5j1bVdvb4xMj9V3AbzG8Wfq2BT5TknSSnTD0q+oB4KWjat+pqtfayweBzcf7jCQbgbdX1YM1vOD7HcAVS2qxJGnJpjGm/6+BPxp5fU6S7yf5kyQfaLVNwIGRZQ602oKS7Ewym2R2bm5uCk2UJMGEoZ/k3wGvAV9vpUPA2VV1HvBp4BtJ3r7Yz62qW6tqUFWDmZmZSZooSRqx5NslJvlN4NeBD7UhG6rqVeDVNr03ybPAu4CD/O0hoM2tJklaRkva009yCfBvgN+oqp+O1GeSnNKm38HwgO1zVXUIeDnJBe2snWuAeyduvSRpUU64p5/kTuBC4MwkB4AbGJ6t8xZgTzvz8sF2ps4Hgc8l+X/AG8Anqmr+IPBvMzwT6BcYHgMYPQ4gSVoGaSMzq9ZgMKjZ2dmVbgYA4/604GT8TVdy3ZLWliR7q2qw0Dx/kStJHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR0x9CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdGSv0k9yW5EiSx0dqZyTZk+SZ9ryh1ZPkliT7kzya5D0j79nRln8myY7pd0eSdDzj7ul/DbjkqNr1wP1VtQ24v70GuJThDdG3ATuBXTD8kmB4f933AecDN8x/UUiSlsdYoV9VDwAvHVW+HLi9Td8OXDFSv6OGHgROT7IRuBjYU1UvVdVfAnv4+S8SSdJJNMmY/llVdahN/wg4q01vAl4cWe5Aqx2r/nOS7Ewym2R2bm5ugiZKkkZN5UBuVRVQ0/is9nm3VtWgqgYzMzPT+lhJ6t4koX+4DdvQno+0+kFgy8hym1vtWHVJ0jKZJPR3A/Nn4OwA7h2pX9PO4rkA+HEbBvo2cFGSDe0A7kWtJklaJqeOs1CSO4ELgTOTHGB4Fs7ngbuTfAx4AbiyLX4fcBmwH/gp8FGAqnopye8CD7flPldVRx8cliSdRBkOx69eg8GgZmdnV7oZACQZa7mT8TddyXVLWluS7K2qwULz/EWuJHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdWTJoZ/kl5LsG3m8nORTSW5McnCkftnIez6TZH+Sp5NcPJ0uaCmSjPWQtL6MdY/chVTV08B2gCSnAAeBexjeE/dLVfWF0eWTnAtcBbwb+EXgj5O8q6peX2obJEmLM63hnQ8Bz1bVC8dZ5nLgrqp6tap+yPDG6edPaf2SpDFMK/SvAu4ceX1dkkeT3JZkQ6ttAl4cWeZAq0mSlsnEoZ/kzcBvAP+1lXYB72Q49HMIuHkJn7kzyWyS2bm5uUmbKElqprGnfynwSFUdBqiqw1X1elW9AXyFnw3hHAS2jLxvc6v9nKq6taoGVTWYmZmZQhMlSTCd0L+akaGdJBtH5n0YeLxN7wauSvKWJOcA24DvTWH9kqQxLfnsHYAkpwG/Bnx8pPx7SbYDBTw/P6+qnkhyN/Ak8BpwrWfuSNLymij0q+r/AH/vqNpHjrP8TcBNk6xTkrR0/iJXkjpi6EtSRwx9SeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcmugzDajfu7f6q6iS3RJJWB/f0Jakjhr4kdcTQl6SOGPqS1BFDX5I6YuhLUkcmDv0kzyd5LMm+JLOtdkaSPUmeac8bWj1JbkmyP8mjSd4z6folSeOb1p7+P6mq7VU1aK+vB+6vqm3A/e01wKUMb4i+DdgJ7JrS+iVJYzhZwzuXA7e36duBK0bqd9TQg8DpSTaepDZIko4yjdAv4DtJ9ibZ2WpnVdWhNv0j4Kw2vQl4ceS9B1pNkrQMpnEZhvdX1cEkfx/Yk+QHozOrqpIs6joH7ctjJ8DZZ589hSZKkmAKe/pVdbA9HwHuAc4HDs8P27TnI23xg8CWkbdvbrWjP/PWqhpU1WBmZmbSJkqSmolCP8lpSd42Pw1cBDwO7AZ2tMV2APe26d3ANe0snguAH48MA0mSTrJJh3fOAu5pV7M8FfhGVf2PJA8Ddyf5GPACcGVb/j7gMmA/8FPgoxOuX5K0CBOFflU9B/yjBer/G/jQAvUCrp1knZKkpfMXuZLUEUNfkjpi6EtSRwx9SerIur5HrlYf71ssrSz39CWpI4a+JHXE0Jekjhj6ktQRQ1+SOmLoS1JHDH1J6oihL0kdMfQlqSOGviR1xNCXpI4Y+pLUEUNfkjqy5NBPsiXJd5M8meSJJJ9s9RuTHEyyrz0uG3nPZ5LsT/J0koun0QFJ0vgmubTya8DvVNUjSd4G7E2yp837UlV9YXThJOcCVwHvBn4R+OMk76qq1ydogyRpEZa8p19Vh6rqkTb9E+ApYNNx3nI5cFdVvVpVPwT2A+cvdf2SpMWbyph+kq3AecBDrXRdkkeT3JZkQ6ttAl4cedsBjvElkWRnktkks3Nzc9NooiSJKYR+krcC3wQ+VVUvA7uAdwLbgUPAzYv9zKq6taoGVTWYmZmZtImSpGai0E/yJoaB//Wq+hZAVR2uqter6g3gK/xsCOcgsGXk7ZtbTZK0TCY5eyfAV4GnquqLI/WNI4t9GHi8Te8GrkryliTnANuA7y11/ZKkxZvk7J1fBT4CPJZkX6t9Frg6yXaggOeBjwNU1RNJ7gaeZHjmz7WeuSNJy2vJoV9VfwpkgVn3Hec9NwE3LXWdkqTJ+ItcSeqIoS9JHTH0Jakjhr4kdcTQl6SOGPqS1BFDX5I6MsmPs6RVYfjj8BOrqpPcEmn1c09fkjpi6EtSRxzeYfzhAUla69zTl6SOGPqS1BFDX5I64pj+SeAphJJWK0Nfx+UX2LH5t9Fa5PCOJHVk2UM/ySVJnk6yP8n1y71+SerZsoZ+klOA3wcuBc5leD/dc5ezDepXkrEe0nq23GP65wP7q+o5gCR3AZczvFm6pmC9hNZ66cdinIw+j3s8YdrHJzzesXotd+hvAl4ceX0AeN/RCyXZCexsL19J8vQS13cm8BdLfO9Jt8j/5OulL6u6HzD9vqzkF9hq78tJ+Nus+n9fizBJX/7BsWasyrN3qupW4NZJPyfJbFUNptCkFbde+rJe+gH2ZTVaL/2Ak9eX5T6QexDYMvJ6c6tJkpbBcof+w8C2JOckeTNwFbB7mdsgSd1a1uGdqnotyXXAt4FTgNuq6omTuMqJh4hWkfXSl/XSD7Avq9F66QecpL7Eo+eS1A9/kStJHTH0Jakj6zL019OlHpI8n+SxJPuSzK50exYjyW1JjiR5fKR2RpI9SZ5pzxtWso3jOkZfbkxysG2bfUkuW8k2jiPJliTfTfJkkieSfLLV19x2OU5f1uJ2+btJvpfkz1pf/n2rn5PkoZZlf9BOgJlsXettTL9d6uHPgV9j+OOvh4Grq2pN/uo3yfPAoKrW3A9OknwQeAW4o6p+udV+D3ipqj7fvpA3VNW/Xcl2juMYfbkReKWqvrCSbVuMJBuBjVX1SJK3AXuBK4DfZI1tl+P05UrW3nYJcFpVvZLkTcCfAp8EPg18q6ruSvJl4M+qatck61qPe/p/c6mHqvq/wPylHrTMquoB4KWjypcDt7fp2xn+J131jtGXNaeqDlXVI236J8BTDH8pv+a2y3H6subU0Cvt5Zvao4B/Cvxhq09lu6zH0F/oUg9r8h9CU8B3kuxtl6dY686qqkNt+kfAWSvZmCm4Lsmjbfhn1Q+JjEqyFTgPeIg1vl2O6guswe2S5JQk+4AjwB7gWeCvquq1tshUsmw9hv568/6qeg/DK5Ne24YZ1oUaji2u5fHFXcA7ge3AIeDmFW3NIiR5K/BN4FNV9fLovLW2XRboy5rcLlX1elVtZ3ilgvOBf3gy1rMeQ39dXeqhqg625yPAPQz/Maxlh9tY7PyY7JEVbs+SVdXh9h/1DeArrJFt08aMvwl8vaq+1cprcrss1Je1ul3mVdVfAd8F/jFwepL5H9FOJcvWY+ivm0s9JDmtHaAiyWnARcDjx3/Xqrcb2NGmdwD3rmBbJjIfks2HWQPbph0w/CrwVFV9cWTWmtsux+rLGt0uM0lOb9O/wPBElKcYhv+/aItNZbusu7N3ANopWv+Rn13q4aaVbdHSJHkHw717GF4y4xtrqS9J7gQuZHiJ2MPADcB/A+4GzgZeAK6sqlV/gPQYfbmQ4RBCAc8DHx8ZF1+Vkrwf+J/AY8AbrfxZhmPha2q7HKcvV7P2tsuvMDxQewrDnfG7q+pzLQPuAs4Avg/8q6p6daJ1rcfQlyQtbD0O70iSjsHQl6SOGPqS1BFDX5I6YuhLUkcMfUnqiKEvSR35/9fp7zCB14xOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(np.concatenate(distances),bins=np.linspace(0,30,31),color='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "id": "59da0d5a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 4480/4480 [13:03<00:00,  5.72it/s]\n"
     ]
    }
   ],
   "source": [
    "rep_traces = df.loc[[u in rep_units for u in df['unit_id']],'norm_trace']\n",
    "rep_corrs = []\n",
    "for i in tqdm(np.arange(len(rep_traces)-1)):\n",
    "    \n",
    "    rep_corrs.append((np.vstack(rep_traces)[i][None] * \n",
    "                       np.vstack(rep_traces)[i+1:][None]).mean(-1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "id": "c27819e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 823/823 [00:21<00:00, 37.66it/s]\n",
      "100%|██████████| 822/822 [00:21<00:00, 37.84it/s]\n",
      "100%|██████████| 810/810 [00:22<00:00, 36.32it/s]\n",
      "100%|██████████| 837/837 [00:23<00:00, 35.74it/s]\n",
      "100%|██████████| 848/848 [00:23<00:00, 35.34it/s]\n",
      "100%|██████████| 859/859 [00:24<00:00, 34.59it/s]\n",
      "100%|██████████| 857/857 [00:24<00:00, 34.86it/s]\n",
      "100%|██████████| 839/839 [00:23<00:00, 35.59it/s]\n",
      "100%|██████████| 837/837 [00:23<00:00, 35.54it/s]\n",
      "100%|██████████| 830/830 [00:23<00:00, 35.95it/s]\n"
     ]
    }
   ],
   "source": [
    "field_corrs = []\n",
    "\n",
    "for field in np.arange(10)+1:\n",
    "    field_traces = df.loc[df['field'] == field,'norm_trace']\n",
    "    for i in tqdm(np.arange(len(field_traces)-1)):\n",
    "\n",
    "        field_corrs.append((np.vstack(field_traces)[i][None] * \n",
    "                           np.vstack(field_traces)[i+1:][None]).mean(-1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "id": "0e66bd82",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAa1ElEQVR4nO3de3RU9d3v8feXEIOIBY/Ep1Iwg30QA0kIENRUi+FWWMqCpce6pGLlYqlYLO2peKmwwJ521WfpsVZLPUVteXzqBY1aWb1qISlKVUhQNAZQjyQKooQIVFCQwPf8kTFPQm6TmZ2ZbPN5rTWLuezZv+8vEz7Z89t7/7a5OyIiEj49Ul2AiIjERwEuIhJSCnARkZBSgIuIhJQCXEQkpHoms7H+/ft7JBJJZpMiIqFXXl6+x90zj38+qQEeiUQoKytLZpMiIqFnZtUtPa8hFBGRkFKAi4iElAJcRCSkkjoGLiLtO3LkCDt27ODQoUOpLkWSrFevXgwcOJD09PSYlleAi3QxO3bs4OSTTyYSiWBmqS5HksTdqa2tZceOHQwePDim92gIRaSLOXToEKeeeqrCu5sxM0499dQOffNSgIt0QQrv7qmjn7sCXEQkpBTgIl1dJAJmwd1iOBu6T58+bb5eVVVFTk5Oh7oxa9YsiouLY14+1jY6ut7W2igrK+P73/9+XOu5++67+eSTT+J6byJCsxMz8uKLVB8+HNd7szIyqCosDLgikSSproYgL7yi4ZkWFRQUUFBQENd77777bmbOnEnv3r0DrqptodkCrz58GC8qiusWb/CLdHcHDhxgwoQJjBo1itzcXJ555pmG1+rq6rjyyivJzs7msssua9gCLS8v58ILL2T06NFMnjyZXbt2NVtva8uUl5czYsQIRowYwfLly1usyd1ZsGABQ4cOZeLEiezevbvhtUgkwp49e4D6LeqioiIAli1bxlVXXUVhYSFDhgzh/vvvb7be0tJSpk6d2tDv2bNnk5ubS15eHk8++SQA8+fPp6CggOHDh7N06VIA7rnnHt5//33GjRvHuHHjAHj22WcpLCxk1KhRfPOb3+TAgQMA3HzzzQwbNoy8vDxuuOGGGD+FNrh70m6jR4/2eFFS4kB8t5KSuNsVSbbKysqmT0CwDcSwvpNOOsnd3Y8cOeL79+93d/eamhr/6le/6seOHfPt27c74C+88IK7u8+ePdvvuOMO/+yzz7ywsNB3797t7u6PPfaYz549293dr776an/iiSfaXCY3N9f/8Y9/uLv7DTfc4MOHD29W25NPPukTJ070uro637lzp/ft29efeOIJd3fPysrympoad3ffuHGjX3jhhe7uvnTpUs/Ly/NPPvnEa2pqfODAgb5z507fvn17QxslJSV+8cUXu7v7jTfe6AsXLmxo86OPPnJ399raWnd3r6ur8wsvvNA3b97crN2amhr/+te/7gcOHHB399tvv91vu+0237Nnj5911ll+7Ngxd3ffu3dviz/7Zp+/uwNl3kKmhmYIBer/2MTDSkuDLUSkm3B3fvzjH7Nu3Tp69OjBzp07+fDDDwEYNGgQ559/PgAzZ87knnvuYcqUKVRUVDBp0iQAjh49yumnn95kndu2bWtxmX379rFv3z7Gjh0LwFVXXcVf/vKXZjWtW7eOGTNmkJaWxoABAxg/fnxMfZk+fTonnngiJ554IuPGjWPDhg3k5+e3uOzf//53HnvssYbHp5xyCgCPP/44K1asoK6ujl27dlFZWUleXl6T97700ktUVlY2/Gw+++wzCgsL6du3L7169WLu3LlMnTq1YWs/EaEKcBFJrocffpiamhrKy8tJT08nEok0HKd8/CFvZoa7M3z4cF588cVW19naMvv27Uu43p49e3Ls2DGAZsdTt1RvR2zfvp0777yTjRs3csoppzBr1qwWj9l2dyZNmsSjjz7a7LUNGzawZs0aiouL+dWvfsXatWs7VMPxQjMGLiLJt3//fk477TTS09MpKSmhuvq/ZzV99913G0L4kUce4YILLmDo0KHU1NQ0PH/kyBHeeOONJutsbZl+/frRr18/XnjhBaD+j0dLxo4dy6pVqzh69Ci7du2ipKSk4bVIJEJ5eTlAw7j155555hkOHTpEbW0tpaWljBkzptV+T5o0qckY/N69e/nXv/7FSSedRN++ffnwww+bfDs4+eST+fjjjwE477zzWL9+PW+//TYABw8e5M033+TAgQPs37+fiy66iF/84hds3ry51fZjpQAX6eqysoI9jDArK+amr7zySsrKysjNzeWhhx7i7LPPbnht6NChLF++nOzsbPbu3cv8+fM54YQTKC4u5qabbmLEiBHk5+fzz3/+s8k621rmd7/7Hd/73vfIz89vdcj0kksuYciQIQwbNoxvf/vbFDY6wmzp0qUsXLiQgoIC0tLSmrwvLy+PcePGcd5557FkyRIGDBjQar8XL17M3r17ycnJYcSIEZSUlDBixAhGjhzJ2Wefzbe+9a2GIRKAefPmMWXKFMaNG0dmZiYrV65kxowZ5OXlUVhYyNatW/n444+ZOnUqeXl5XHDBBdx1110xfw6tsXjHleNRUFDg8V7QwUpL8ege5WS+VyTZtmzZQnZ2dqrL+EJZtmwZffr0CebIj07W0udvZuXu3uwYx3a3wM3st2a228wqGj13h5ltNbPXzOxpM+sXROEiIhK7WIZQVgJTjnvuOSDH3fOAN4FbAq5LRCQwy5YtC8XWd0e1G+Duvg746LjnnnX3uujDl4CBnVCbiIi0IYidmHOA5gdrRpnZPDMrM7OympqaAJoTERFIMMDN7FagDmj5eB/A3Ve4e4G7F2RmZibSnIiINBL3iTxmNguYCkzwZB7KIiIiQJxb4GY2BbgRmObuyZ9DUaQb+fzSakHdIl+w6WRjsXr1am6//fZA1tWVtLsFbmaPAkVAfzPbASyl/qiTDOC56OmoL7n7tZ1Yp0i3VV1dHfc8QC3pjlf7mTZtGtOmTUt1GYGL5SiUGe5+urunu/tAd3/Q3f/d3Qe5e370pvAW+QLqitPJNp72FWDBggWsXLkSqP+2snTp0oZ6t27dCsDKlStZsGABUD+nSWFhIbm5uSxevLjh20Zb642lT6mgU+lFpFW9evXi6aefZtOmTZSUlPCjH/2o4dvAtm3buO6669iyZQtf+tKX+PWvf82RI0e4/vrrKS4upry8nDlz5nDrrbc2WWdby8yePZt77703oXlC+vfvz6ZNm5g/fz533nlns9cXLlzI/Pnzef3115vNlNiSWPqUKpqNUERa1RWnk23PpZdeCsDo0aN56qmnmr2+fv36homurrrqKm666aY219davV2BAlxEWtUVp5NtPGUsNJ82NiMjA4C0tDTq6upoSUv7AVpbbyx9ShUNoYhIq7ridLJZWVlUVlZy+PBh9u3bx5o1azrUp/PPP7/hYg2N22htvbH0KVUU4CJdXFZWVqCHEWaFfDrZQYMGcfnll5OTk8Pll1/OyJEjO/Tz/OUvf8ny5cvJzc1l586d7a43lj6liqaTFeliNJ1scvXp06fhosNdQaDTyYqISNekABeRbq0rbX13lAJcRCSkFOAiIiGlABcRCSkFuIhISOlMTJEuLvLii1QfPhzY+rIyMqgqLAxsfUF6/vnnufbaa0lPT+dPf/oTCxcubHMK2qqqKqZOnUpFRUWz14qKirjzzjspKGh29F2Hfe1rX+syx343pgAX6eKqDx8O9DwGKy0NbF1Be/jhh7nllluYOXMmQIfmD+9MXTG8QUMoInKcgwcPcvHFFzNixAhycnJYtWoVAD/5yU8YM2YMOTk5zJs3r+FMyaKiIn74wx9SUFBAdnY2Gzdu5NJLL2XIkCEsXry4Yb2///3vOeecc8jPz+e73/0uR48ebdLuAw88wOOPP86SJUu48sorm1zQ4ejRoyxatIgxY8aQl5fHb37zm2Z1f/rpp1xxxRVkZ2dzySWX8Omnn7bYv0gkwp49ewAoKyujKPrHcdmyZcyZM4eioiLOPPNM7rnnnob3fD7lrLuzYMEChg4dysSJE7nooosa/si0tt6DBw8yZ84czjnnHEaOHNlkSt5EKcBFpIm//vWvDBgwgM2bN1NRUcGUKVOA+vmxN27cSEVFBZ9++il//OMfG95zwgknUFZWxrXXXsv06dNZvnw5FRUVrFy5ktraWrZs2cKqVatYv349r776Kmlpac3mOrnmmmuYNm0ad9xxR7PXHnzwQfr27cvGjRvZuHEj999/P9u3b2+yzH333Ufv3r3ZsmULt912G+Xl5R3u+9atW/nb3/7Ghg0buO222zhy5EiT159++mm2bdtGZWUlDz30UExb5j/72c8YP348GzZsoKSkhEWLFnHw4MEO19YSBbiINJGbm8tzzz3HTTfdxPPPP0/fvn0BKCkp4dxzzyU3N5e1a9c2mdDp86vd5ObmMnz4cE4//XQyMjI488wzee+991izZg3l5eWMGTOG/Px81qxZwzvvvBNzTc8++ywPPfQQ+fn5nHvuudTW1vLWW281WWbdunUNQy95eXnk5eV1uO8XX3wxGRkZ9O/fn9NOO61h6tzGbcyYMYO0tDQGDBjA+PHjY6r99ttvJz8/n6KiIg4dOsS7777b4dpaojFwEWnirLPOYtOmTfz5z39m8eLFTJgwgRtvvJHrrruOsrIyBg0axLJly5pM4/r5FK49evRouP/547q6Otydq6++mp///Odx1eTu3HvvvUyePLnJ81VVVR1eV+NpY1ubihbano62I+t1d5588kmGDh3a4Vrboy1wEWni/fffp3fv3sycOZNFixaxadOmhkDq378/Bw4c6PDOxQkTJlBcXMzu3bsB+Oijj5pMTdueyZMnc9999zUMabz55pvNhiHGjh3LI488AkBFRQWvvfZai+uKRCINwyufX9ghVmPHjmXVqlUcPXqUXbt2UVJS0u56J0+ezL333tuwz+CVV17pUJtt0Ra4SBeXlZER6JEjWY22Mlvy+uuvs2jRInr06EF6ejr33Xcf/fr14zvf+Q45OTl8+ctfZsyYMR1qc9iwYfz0pz/lG9/4BseOHSM9PZ3ly5fHPLXtNddcQ1VVFaNGjcLdyczM5A9/+EOTZebPn8/s2bPJzs4mOzub0aNHt7iupUuXMnfuXJYsWdKwozFWl1xyCWvXrmXYsGGcccYZFDY6HLO19S5ZsoQf/OAH5OXlcezYMQYPHtxk/0EiNJ2sSBej6WTDY9asWUydOpXLLrsssHVqOlkRkW5AQygiInFauXJlSttvdwvczH5rZrvNrKLRc//DzJ4zs7ei/57SuWWKdC/JHNqUrqOjn3ssQygrgSnHPXczsMbdhwBroo9FJAC9evWitrZWId7NuDu1tbX06tUr5ve0O4Ti7uvMLHLc09OBouj9/wRKgZtiblVEWjVw4EB27NhBTU1NqkuRJOvVqxcDBw6Mefl4x8D/zd13Re9/APxbawua2TxgHsAZZ5wRZ3Mi3Ud6ejqDBw9OdRkSAgkfheL13/Na/a7n7ivcvcDdCzIzMxNtTkREouIN8A/N7HSA6L+7gytJRERiEW+Arwaujt6/GghufkQREYlJLIcRPgq8CAw1sx1mNhe4HZhkZm8BE6OPRUQkiWI5CmVGKy9NCLgWERHpAJ1KLyISUgpwEZGQUoCLiISUAlxEJKQU4CIiIaUAFxEJKQW4iHQfkQiYxXeLRFJdfTO6oIOIdB/V1RDvNL1mwdYSAG2Bi4iElAJcRCSkFOAiIiGlABcRCSkFuIhISCnARURCSgEuIhJSCnARkZBSgIuIhJQCXEQkpBTgIiIhpQAXEQkpBbiISEgpwEVEQkoBLiISUgpwEZGQSijAzeyHZvaGmVWY2aNm1iuowkREpG1xB7iZfQX4PlDg7jlAGnBFUIWJiEjbEh1C6QmcaGY9gd7A+4mXJCIisYg7wN19J3An8C6wC9jv7s8ev5yZzTOzMjMrq6mpib9SERFpIpEhlFOA6cBgYABwkpnNPH45d1/h7gXuXpCZmRl/pSIi0kQiQygTge3uXuPuR4CngK8FU5aIiLQnkQB/FzjPzHqbmQETgC3BlCUiIu1JZAz8ZaAY2AS8Hl3XioDqEhGRdvRM5M3uvhRYGlAtIiLSAToTU0QkpBTgIiIhpQAXEQkpBbiISEgpwEVEQkoBLiISUgpwEZGQUoCLiISUAlxEJKQU4CIiIaUAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkFKAi4iElAJcRCSkFOAiIiGlABcRCSkFuIhILLKywCz+WyQSeEkJXdRYRKTbqKpK7P1mgZTRmLbARURCSgEuIhJSCnARkZBKKMDNrJ+ZFZvZVjPbYmaFQRUmIiJtS3Qn5i+Bv7r7ZWZ2AtA7gJpERCQGcQe4mfUFxgKzANz9M+CzYMoSEZH2JDKEMhioAX5nZq+Y2QNmdtLxC5nZPDMrM7OympqaBJoTEZHGEgnwnsAo4D53HwkcBG4+fiF3X+HuBe5ekJmZmUBzIiLSWCIBvgPY4e4vRx8XUx/oIiKSBHEHuLt/ALxnZkOjT00AKgOpSkRE2pXoUSjXAw9Hj0B5B5ideEkiIhKLhALc3V8FCoIpRUREOkJnYoqIhJQCXEQkpBTgIiIhpQAXEQkpBbiISEgpwEVEQkoBLiISUt0mwM0s7lukEy5GKiKSqG5xUeOsjAyqS0rifn/1Bx8EWI2ISDC6RYBXFSZ2oSArLQ2mEBGRAHWbIRQRkS8aBbiISEgpwEVEQkoBLiISUgpwEZGQUoCLiISUAlxEJKQU4CIiIaUAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkEo4wM0szcxeMbM/BlGQiIjEJogt8IXAlgDWIyIiHZBQgJvZQOBi4IFgyhERkVglugV+N3AjcKy1BcxsnpmVmVlZTU1Ngs2JiMjn4g5wM5sK7Hb38raWc/cV7l7g7gWZmZnxNiciIsdJZAv8fGCamVUBjwHjzez3gVQlIiLtijvA3f0Wdx/o7hHgCmCtu88MrDIREWmTjgMXEQmpQK5K7+6lQGkQ6xIRkdhoC1xEJKQU4CIiIaUAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkFKAi4iElAJcRCSkFOAiIiGlABcRCSkFuIhISCnARURCSgEuIhJSCnARCZdIBMziu2Vlpbr6QAUyH7iISNJUV4N7qqvoErQFLiISUgpwEZGQUoCLiISUAlxEJKQU4CIiIaUAFxEJKQW4iEhIKcBFREIq7gA3s0FmVmJmlWb2hpktDLIwERFpWyJnYtYBP3L3TWZ2MlBuZs+5e2VAtYmISBvi3gJ3913uvil6/2NgC/CVoAoTEZG2BTIGbmYRYCTwcguvzTOzMjMrq6mpCaI5EREhgAA3sz7Ak8AP3P1fx7/u7ivcvcDdCzIzMxNtTkREohIKcDNLpz68H3b3p4IpSUREYpHIUSgGPAhscfe7gitJRERikcgW+PnAVcB4M3s1ersooLpERKQdcR9G6O4vABZgLSIi0gE6E1NEJKQU4DEys7hukUgk1aWLyBeUrokZg6yMDKpLSuJ6b/UHHwRcjYhIPQV4DKoKC+N+r5WWBleIiEgjGkIREQkpBbiISEgpwEVEQkoBLiISUgpwEZGQUoCLiISUAlxEJKQU4CIiIaUAF5Hki0TALL5bVlaqq+8ydCamiCRfdTW4p7qK0NMWuIhISCnARURCSgEuIhJSCnARkZBSgIuIhJQCPAl0NR8R6Qw6jLCT6Wo+ItJZFOCdTFfzEZHOoiEUEZGQUoCLSMclciq8TocPTEIBbmZTzGybmb1tZjcHVZSIdHGfnwof762qKtU9+EKIewzczNKA5cAkYAew0cxWu3tlUMV1d2l79sQ9Dp62Zw91l10WbEHyxRKJ1AdxPLQF3SUkshPzHOBtd38HwMweA6YDCvCAJBLAPYuLFf4dkUiYhVVWliaUCrlEAvwrwHuNHu8Azj1+ITObB8yLPjxgZtvibK+/wZ443xtW/UlBn48CluxG/1tK+pxiqelzdXX9eHRqdM/P2SzePrf4lafTDyN09xXAikTXY2Zl7l4QQEmhoT53D+pz99AZfU5kJ+ZOYFCjxwOjz4mISBIkEuAbgSFmNtjMTgCuAFYHU5aIiLQn7iEUd68zswXA34A04Lfu/kZglTWX8DBMCKnP3YP63D0E3mdz7YUWEQklnYkpIhJSCnARkZDqcgHe3un5ZpZhZquir79sZpEUlBmoGPr8v8ys0sxeM7M1Zhb60+BinYbBzP6nmbmZhfqQs1j6a2aXRz/nN8zskWTXGLQYfq/PMLMSM3sl+rt9USrqDJKZ/dbMdptZRSuvm5ndE/2ZvGZmoxJq0N27zI36naH/DzgTOAHYDAw7bpnrgP8bvX8FsCrVdSehz+OA3tH787tDn6PLnQysA14CClJddyd/xkOAV4BToo9PS3XdSejzCmB+9P4woCrVdQfQ77HAKKCildcvAv5C/bly5wEvJ9JeV9sCbzg9390/Az4/Pb+x6cB/Ru8XAxPMUnc6WQDa7bO7l7j7J9GHL1F/zH2YxfI5A/xv4D+AQ8ksrhPE0t/vAMvdfS+Au+9Oco1Bi6XPDnwper8v8H4S6+sU7r4O+KiNRaYDD3m9l4B+ZnZ6vO11tQBv6fT8r7S2jLvXAfuBU5NSXeeIpc+NzaX+L3iYtdvn6FfLQe7+p2QW1kli+YzPAs4ys/Vm9pKZTUladZ0jlj4vA2aa2Q7gz8D1ySktpTr6/71NuiJPiJjZTKAAuDDVtXQmM+sB3AXMSnEpydST+mGUIuq/Ya0zs1x335fKojrZDGClu/8fMysE/svMctz9WKoLC4uutgUey+n5DcuYWU/qv3rVJqW6zhHTlARmNhG4FZjm7oeTVFtnaa/PJwM5QKmZVVE/Vrg6xDsyY/mMdwCr3f2Iu28H3qQ+0MMqlj7PBR4HcPcXgV7UT3L1RRboFCRdLcBjOT1/NXB19P5lwFqP7h0IqXb7bGYjgd9QH95hHxuFdvrs7vvdvb+7R9w9Qv24/zR3L0tNuQmL5ff6D9RvfWNm/akfUnkniTUGLZY+vwtMADCzbOoDvCapVSbfauDb0aNRzgP2u/uuuNeW6r22reylfZP6Pdi3Rp/7CfX/gaH+Q34CeBvYAJyZ6pqT0Oe/Ax8Cr0Zvq1Ndc2f3+bhlSwnxUSgxfsZG/bBRJfA6cEWqa05Cn4cB66k/QuVV4BuprjmAPj8K7AKOUP+tai5wLXBto895efRn8nqiv9c6lV5EJKS62hCKiIjESAEuIhJSCnARkZBSgIuIhJQCXEQkpBTgIiIhpQAXEQmp/w/wSX035yXeNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.linspace(0,1,21)\n",
    "plt.hist(np.concatenate(corrs),bins=bins,color='r',histtype=u'step',\n",
    "         density=True,label='labeled duplicates')\n",
    "plt.hist(np.hstack(rep_corrs).squeeze(),bins=bins,color='k',histtype=u'step',\n",
    "         density=True, label='labeled unique')\n",
    "plt.hist(np.hstack(field_corrs).squeeze(),bins=bins,color='c',histtype=u'step',\n",
    "         density=True, label='same field unique')\n",
    "_ = plt.legend()"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}