nanshe-benchmarking.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Provide ability to parse markdown into html\n",
      "# Borrowed code from < http://stackoverflow.com/questions/21364102/expand-variables-in-markdown-cells-of-ipython-notebook >.\n",
      "\n",
      "from __future__ import absolute_import\n",
      "\n",
      "from IPython.core.getipython import get_ipython\n",
      "from IPython.core.magic import Magics, magics_class, cell_magic, line_magic\n",
      "from IPython.core.display import HTML as _HTML\n",
      "\n",
      "from markdown import markdown, markdownFromFile\n",
      "from StringIO import StringIO\n",
      "\n",
      "\n",
      "class HTML(_HTML):\n",
      "    def __repr__(self):\n",
      "        return \"\"\n",
      "\n",
      "@magics_class\n",
      "class MarkdownMagics(Magics):\n",
      "    @cell_magic\n",
      "    def markdown(self, line, cell):\n",
      "        vars = line.split()\n",
      "\n",
      "        d = {}\n",
      "        for k, v in self.shell.user_ns.items():\n",
      "            if k in vars:\n",
      "                d[k] = v\n",
      "\n",
      "        return HTML(\"<p>{}</p>\".format(markdown(cell.format(**d))))\n",
      "    \n",
      "    @line_magic\n",
      "    def markdown_file(self, filename):\n",
      "        \"\"\"\n",
      "            Read a markdown file and convert it to HTML for display.\n",
      "        \"\"\"\n",
      "\n",
      "        # Provide buffer to read and convert file.\n",
      "        sb = StringIO()\n",
      "        markdownFromFile(filename, sb)\n",
      "        s = sb.getvalue()\n",
      "        \n",
      "        return HTML(s)\n",
      "    \n",
      "\n",
      "get_ipython().register_magics(MarkdownMagics)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%markdown_file Readme.md"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<h1>Goal</h1>\n",
        "<ul>\n",
        "<li>produce a label matrix with output</li>\n",
        "<li>evaluate efficacy</li>\n",
        "</ul>\n",
        "<h2>Ground Truth</h2>\n",
        "<ul>\n",
        "<li>Ground truth will be stored in files named anXXXXX_YYYY_MM_DD_fov_NNNNN_rois.mat ( where all repeated capital letters are decimal digits with XXXXX indicates the animal, YYYY_MM_DD indicates the date in year month day format, and NNNNN represents a field of view ).</li>\n",
        "<li>These files will contain a group called <code>roi</code>. Inside of that group will be the following fields.<ul>\n",
        "<li><code>mask</code> - An array where each index represents a different ROI. All other indices represent spatial coordinates.</li>\n",
        "<li><code>centroid_xy</code> - A matrix where the first index will represent the relevant ROI. The last will represent which axis of the centroid is being represented.</li>\n",
        "<li><code>burst_rate</code> - A vector where the index represents the relevant ROI.</li>\n",
        "</ul>\n",
        "</li>\n",
        "<li>In addition, the raw data will be supplied in TIFFs named Image_Registration_4_anXXXXX_YYYY_MM_DD_main_FFF.tif ( where all repeated capital letters are decimal digits with XXXXX indicates the animal, YYYY_MM_DD indicates the date in year month day format, and FFF represents the frame )</li>\n",
        "</ul>\n",
        "<h2>Evaluating Efficacy</h2>\n",
        "<p>Two metrics will be used.</p>\n",
        "<ol>\n",
        "<li>Distance between centroids of detected and ground truth ROIs (i.e. as measured by the L<sub>2</sub> or Euclidean norm).</li>\n",
        "<li>Percentage overlap of detected and ground truth ROIs (i.e. mask intersection divided by union).</li>\n",
        "</ol>\n",
        "<p>The results will be broked down into 3 categories. The determination of which category they fall into will dependent on a chosen acceptance threshold for the relevant metric.</p>\n",
        "<ol>\n",
        "<li>True Positives - How many matches are there between detected and ground truth ROIs</li>\n",
        "<li>False Negatives - How many ground truth ROIs are missed?</li>\n",
        "<li>False Positives - How many detected ROIs are incorrect (i.e. do not correspond to one in the ground truth).</li>\n",
        "</ol>\n",
        "<p>After assessing both metrics, the next relevant question is how does activity/noise affect detectability? This will require plots of these categories against the burst rate.</p>\n",
        "<p>Finally, these metrics must be stored in some readable format for future use.</p>\n",
        "<h2>Notes</h2>\n",
        "<p>Data is almost always 2 channel. Only use the green channel (1st one).</p>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": []
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Just some imports to get us going\n",
      "\n",
      "import os\n",
      "\n",
      "import operator\n",
      "import operator as op\n",
      "\n",
      "import itertools\n",
      "import itertools as it\n",
      "\n",
      "import numpy\n",
      "import numpy as np\n",
      "\n",
      "import scipy\n",
      "import scipy as sp\n",
      "\n",
      "import scipy.spatial\n",
      "import scipy.spatial as sps\n",
      "\n",
      "import scipy.spatial.distance\n",
      "import scipy.spatial.distance as spsd\n",
      "\n",
      "import h5py\n",
      "import h5py as hp\n",
      "\n",
      "import matplotlib\n",
      "import matplotlib as mpl\n",
      "\n",
      "import matplotlib.pyplot\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "import matplotlib.widgets\n",
      "import matplotlib.widgets as wgt\n",
      "\n",
      "#import nanshe.nanshe.expanded_numpy as enp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Borrowed from nanshe to make it easier to make the notebook self contained.\n",
      "\n",
      "\n",
      "def add_singleton_axis_pos(a_array, new_axis = 0):\n",
      "    \"\"\"\n",
      "        Adds a singleton axis to the given position.\n",
      "        Allows negative values for new_axis.\n",
      "        Also, automatically bounds new_axis in an acceptable regime if it is not already.\n",
      "        \n",
      "        Args:\n",
      "            a_array(numpy.ndarray):            array to add the singleton axis to.\n",
      "            new_axis(int):                     position for the axis to be in the final array (defaults to zero).\n",
      "        \n",
      "        Returns:\n",
      "            (numpy.ndarray):                   a numpy array with the singleton axis added (should be a view).\n",
      "        \n",
      "        Examples:\n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6))).shape\n",
      "            (1, 7, 9, 6)\n",
      "            \n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6)), 0).shape\n",
      "            (1, 7, 9, 6)\n",
      "            \n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6)), new_axis = 0).shape\n",
      "            (1, 7, 9, 6)\n",
      "            \n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6)), new_axis = 1).shape\n",
      "            (7, 1, 9, 6)\n",
      "            \n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6)), new_axis = 2).shape\n",
      "            (7, 9, 1, 6)\n",
      "            \n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6)), new_axis = 3).shape\n",
      "            (7, 9, 6, 1)\n",
      "            \n",
      "            >>> add_singleton_axis_pos(numpy.ones((7,9,6)), new_axis = -1).shape\n",
      "            (7, 9, 6, 1)\n",
      "            \n",
      "    \"\"\"\n",
      "\n",
      "    # Clean up new_axis to be within the allowable range.\n",
      "    new_axis %= (a_array.ndim + 1)\n",
      "    if new_axis < 0:\n",
      "        new_axis += a_array.ndim + 1\n",
      "\n",
      "    # Constructing the current ordering of axis and the singleton dime\n",
      "    new_array_shape = range(1, a_array.ndim + 1)\n",
      "    new_array_shape.insert(new_axis, 0)\n",
      "    new_array_shape = tuple(new_array_shape)\n",
      "\n",
      "    # Adds singleton dimension at front.\n",
      "    # Then changes the order so it is elsewhere.\n",
      "    new_array = a_array[None]\n",
      "    new_array = new_array.transpose(new_array_shape)\n",
      "\n",
      "    return( new_array )\n",
      "\n",
      "def add_singleton_op(op, new_array, axis):\n",
      "    \"\"\"\n",
      "        Performs an operation on the given array on the specified axis, which otherwise would have eliminated the axis\n",
      "        in question. This function will instead ensure that the given axis remains after the operation as a singleton.\n",
      "\n",
      "        Note:\n",
      "            The operation must be able to take only two arguments where the first is the array and the second is the\n",
      "            axis to apply the operation along.\n",
      "\n",
      "        Args:\n",
      "            op(callable):                 callable that takes a numpy.ndarray and an int in order.\n",
      "            new_array(numpy.ndarray):     array to perform operation on and add singleton axis too.\n",
      "            axis(int):                    the axis to apply the operation along and turn into a singleton.\n",
      "\n",
      "        Returns:\n",
      "            (numpy.ndarray):              the array with the operation performed.\n",
      "\n",
      "        Examples:\n",
      "            >>> add_singleton_op(numpy.max, numpy.ones((7,9,6)), 0).shape\n",
      "            (1, 9, 6)\n",
      "\n",
      "            >>> add_singleton_op(numpy.max, numpy.ones((7,9,6)), -1).shape\n",
      "            (7, 9, 1)\n",
      "\n",
      "            >>> add_singleton_op(numpy.min, numpy.ones((7,9,6)), 1).shape\n",
      "            (7, 1, 6)\n",
      "\n",
      "            >>> add_singleton_op(numpy.mean, numpy.ones((7,9,6)), 1).shape\n",
      "            (7, 1, 6)\n",
      "    \"\"\"\n",
      "\n",
      "    return(add_singleton_axis_pos(op(new_array, axis), axis))\n",
      "\n",
      "\n",
      "def threshold_array(an_array, threshold, include_below=True, is_closed=True):\n",
      "    \"\"\"\n",
      "        Given a threshold, this function compares the given array to it to see which entries match.\n",
      "\n",
      "        Args:\n",
      "            an_array(numpy.ndarray):                          an array to threshold.\n",
      "            threshold(int or float or numpy.ndarray):         something to compare to.\n",
      "            include_below(bool):                              whether values below the threshold count or ones above it.\n",
      "            is_closed(bool):                                  whether to include values equal to the threshold\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):                               a mask of entries reflecting the threshold.\n",
      "\n",
      "        Examples:\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 0, include_below=True, is_closed=False)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 0, include_below=True, is_closed=True)\n",
      "            array([[ True, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 0, include_below=False, is_closed=True)\n",
      "            array([[ True,  True,  True],\n",
      "                   [ True,  True,  True]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 0, include_below=False, is_closed=False)\n",
      "            array([[False,  True,  True],\n",
      "                   [ True,  True,  True]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 5, include_below=True, is_closed=True)\n",
      "            array([[ True,  True,  True],\n",
      "                   [ True,  True,  True]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 6, include_below=True, is_closed=False)\n",
      "            array([[ True,  True,  True],\n",
      "                   [ True,  True,  True]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3), 2*numpy.ones((2,3)), include_below=True, is_closed=True)\n",
      "            array([[ True,  True,  True],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3) % 2, 1, include_below=False, is_closed=True)\n",
      "            array([[False,  True, False],\n",
      "                   [ True, False,  True]], dtype=bool)\n",
      "\n",
      "            >>> threshold_array(numpy.arange(6).reshape(2,3) % 2, 1, include_below=True, is_closed=False)\n",
      "            array([[ True, False,  True],\n",
      "                   [False,  True, False]], dtype=bool)\n",
      "    \"\"\"\n",
      "\n",
      "    accepted = numpy.zeros(an_array.shape, dtype=bool)\n",
      "\n",
      "    if include_below and is_closed:\n",
      "        accepted[:] = (an_array <= threshold)\n",
      "    elif not include_below and is_closed:\n",
      "        accepted[:] = (an_array >= threshold)\n",
      "    elif include_below and not is_closed:\n",
      "        accepted[:] = (an_array < threshold)\n",
      "    else:\n",
      "        accepted[:] = (an_array > threshold)\n",
      "\n",
      "    return(accepted)\n",
      "\n",
      "\n",
      "def unique_mapping(mapping, out=None):\n",
      "    \"\"\"\n",
      "        Take a binary mapping between two sets and excludes portions of the mapping that are not one-to-one.\n",
      "\n",
      "        Args:\n",
      "            mapping(numpy.ndarray):      bool array mapping between to sets.\n",
      "            out(numpy.ndarray):          where to store the results.\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):          the results returned.\n",
      "\n",
      "        Examples:\n",
      "            >>> unique_mapping(numpy.zeros((2,2), dtype=bool))\n",
      "            array([[False, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "\n",
      "            >>> unique_mapping(numpy.ones((2,2), dtype=bool))\n",
      "            array([[False, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "\n",
      "            >>> unique_mapping(numpy.eye(2, dtype=bool))\n",
      "            array([[ True, False],\n",
      "                   [False,  True]], dtype=bool)\n",
      "\n",
      "            >>> unique_mapping(\n",
      "            ...     numpy.array([[ True,  True],\n",
      "            ...                  [ True, False]], dtype=bool)\n",
      "            ... )\n",
      "            array([[False, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "\n",
      "            >>> unique_mapping(\n",
      "            ...     numpy.array([[ True, False],\n",
      "            ...                  [False, False]], dtype=bool)\n",
      "            ... )\n",
      "            array([[ True, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "\n",
      "            >>> a = numpy.ones((2,2), dtype=bool); a\n",
      "            array([[ True,  True],\n",
      "                   [ True,  True]], dtype=bool)\n",
      "\n",
      "            >>> unique_mapping(a, out=a)\n",
      "            array([[False, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "\n",
      "            >>> a\n",
      "            array([[False, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "    \"\"\"\n",
      "\n",
      "    assert(issubclass(mapping.dtype.type, numpy.bool_))\n",
      "    assert(mapping.ndim == 2)\n",
      "\n",
      "    if out is None:\n",
      "        out = mapping.copy()\n",
      "    elif id(mapping) != id(out):\n",
      "        out[:] = mapping\n",
      "\n",
      "    injective_into = list()\n",
      "    for i in xrange(mapping.ndim):\n",
      "        injective_into_i = (add_singleton_op(numpy.sum, mapping, i) == 1)\n",
      "        injective_into.append(injective_into_i)\n",
      "\n",
      "    for i in xrange(mapping.ndim):\n",
      "        out *= injective_into[i]\n",
      "\n",
      "    return(out)\n",
      "\n",
      "\n",
      "def threshold_metric(a_metric, threshold, include_below=True, is_closed=True):\n",
      "    \"\"\"\n",
      "        Given a threshold, this function finds which entries uniquely match given the threshold.\n",
      "\n",
      "        Args:\n",
      "            a_metric(numpy.ndarray):                          an array to threshold.\n",
      "            threshold(int or float or numpy.ndarray):         something to compare to.\n",
      "            include_below(bool):                              whether values below the threshold count or ones above it.\n",
      "            is_closed(bool):                                  whether to include values equal to the threshold\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):                               a mapping of unique matches.\n",
      "\n",
      "        Examples:\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 0, include_below=True, is_closed=False)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 0, include_below=True, is_closed=True)\n",
      "            array([[ True, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 0, include_below=False, is_closed=True)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 0, include_below=False, is_closed=False)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 5, include_below=True, is_closed=True)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 6, include_below=True, is_closed=False)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), 2*numpy.ones((2,3)), include_below=True, is_closed=True)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3), numpy.arange(0, 12, 2).reshape(2,3), include_below=True, is_closed=True)\n",
      "            array([[False, False, False],\n",
      "                   [False, False, False]], dtype=bool)\n",
      "\n",
      "            >>> threshold_metric(numpy.arange(6).reshape(2,3) % 2, 1, include_below=True, is_closed=False)\n",
      "            array([[False, False, False],\n",
      "                   [False,  True, False]], dtype=bool)\n",
      "    \"\"\"\n",
      "\n",
      "    assert(a_metric.ndim == 2)\n",
      "\n",
      "    accepted = unique_mapping(threshold_array(a_metric, threshold, include_below=include_below, is_closed=is_closed))\n",
      "\n",
      "    return(accepted)\n",
      "\n",
      "\n",
      "def compute_mapping_matches(mapping):\n",
      "    \"\"\"\n",
      "        Given a mapping this function computes number of matches and mismatches found.\n",
      "\n",
      "        In the array returned, first value is the number of true positives (or matches) and then mismatches along each\n",
      "        dimension of the mapping in order.\n",
      "\n",
      "        If axis 0 is the ground truth, then the second and third values are the number of false negatives (or misses)\n",
      "        and false positives (or false alarm)\n",
      "\n",
      "        Args:\n",
      "            mapping(numpy.ndarray):                           a 2D bool array mapping intersections between 2 groups.\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):                               Counts of the number of matches and mismatches.\n",
      "\n",
      "        Examples:\n",
      "            >>> compute_mapping_matches(numpy.arange(6).reshape(2,3) < 0)\n",
      "            array([0, 2, 3], dtype=uint64)\n",
      "\n",
      "            >>> compute_mapping_matches(numpy.arange(6).reshape(2,3) <= 0)\n",
      "            array([1, 1, 2], dtype=uint64)\n",
      "\n",
      "            >>> compute_mapping_matches((numpy.arange(6).reshape(2,3) % 2) == 1)\n",
      "            array([1, 1, 2], dtype=uint64)\n",
      "\n",
      "            >>> compute_mapping_matches(numpy.eye(2, dtype=bool))\n",
      "            array([2, 0, 0], dtype=uint64)\n",
      "\n",
      "            >>> compute_mapping_matches(numpy.fliplr(numpy.eye(2, dtype=bool)))\n",
      "            array([2, 0, 0], dtype=uint64)\n",
      "    \"\"\"\n",
      "\n",
      "    assert(issubclass(mapping.dtype.type, numpy.bool_))\n",
      "    assert(mapping.ndim == 2)\n",
      "\n",
      "    stats = numpy.zeros((3,), dtype=numpy.uint64)\n",
      "\n",
      "    stats[:] = unique_mapping(mapping).sum()\n",
      "    mapping_shape = numpy.array(mapping.shape, dtype=numpy.uint64)\n",
      "    stats[1:] = mapping_shape - stats[1:]\n",
      "\n",
      "    return(stats)\n",
      "\n",
      "\n",
      "def compute_mapping_relevance(mapping):\n",
      "    \"\"\"\n",
      "        Given a mapping this function computes the recall and precision.\n",
      "\n",
      "        If axis 0 is the ground truth, then the returned values are the recall and precision in order.\n",
      "        If axis 1 is the ground truth, then they are flipped.\n",
      "\n",
      "        Args:\n",
      "            mapping(numpy.ndarray):                       a 2D bool array mapping intersections between 2 groups.\n",
      "\n",
      "        Returns:\n",
      "            relevance(tuple of floats):                   relevance - a combination of recall and precision\n",
      "                                                              recall    -   the ratio of relevant predicted positives\n",
      "                                                                            out of all relevant positives.\n",
      "                                                              precision -   the ratio of relevant predicted positives\n",
      "                                                                            out of all predicted positives.\n",
      "\n",
      "        Examples:\n",
      "            >>> compute_mapping_relevance(numpy.arange(6).reshape(2,3) < 0)\n",
      "            array([ 0.,  0.])\n",
      "\n",
      "            >>> compute_mapping_relevance(numpy.arange(6).reshape(2,3) <= 0)\n",
      "            array([ 0.5       ,  0.33333333])\n",
      "\n",
      "            >>> compute_mapping_relevance((numpy.arange(6).reshape(2,3) % 2) == 1)\n",
      "            array([ 0.5       ,  0.33333333])\n",
      "\n",
      "            >>> compute_mapping_relevance(numpy.eye(2, dtype=bool))\n",
      "            array([ 1.,  1.])\n",
      "\n",
      "            >>> compute_mapping_relevance(numpy.fliplr(numpy.eye(2, dtype=bool)))\n",
      "            array([ 1.,  1.])\n",
      "    \"\"\"\n",
      "\n",
      "    assert(issubclass(mapping.dtype.type, numpy.bool_))\n",
      "    assert(mapping.ndim == 2)\n",
      "\n",
      "    matches = compute_mapping_matches(mapping)\n",
      "    mapping_shape = numpy.array(mapping.shape)\n",
      "\n",
      "    relevance = numpy.zeros((2,), dtype=float)\n",
      "    relevance[:] = matches[0]\n",
      "    relevance /= mapping_shape\n",
      "\n",
      "    return(relevance)\n",
      "\n",
      "\n",
      "def matrix_reduced_op(a, b, op):\n",
      "    \"\"\"\n",
      "        Sort of like numpy.dot. However, it will use the first axis with both arrays.\n",
      "        This means they may not need to be matrices. However, they must have the same number of dimensions. Generally,\n",
      "        though not explicitly required, the operator will likely expect every dimension other than the first to be the\n",
      "        same shape.\n",
      "\n",
      "        Args:\n",
      "            a(numpy.ndarray):      first array.\n",
      "            b(numpy.ndarray):      second array.\n",
      "            op(callable):          an operator that will take a[i] and b[j] as arguments and return a scalar.\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):    an array (matrix) with the shape (len(a), len(b)) with each element out[i, j] the\n",
      "                                   result of op(a[i], b[j]).\n",
      "\n",
      "        Examples:\n",
      "            >>> matrix_reduced_op(numpy.arange(6).reshape(2,3), numpy.arange(0,12,2).reshape(2,3), op=numpy.dot)\n",
      "            array([[ 10,  28],\n",
      "                   [ 28, 100]])\n",
      "\n",
      "            >>> numpy.dot(numpy.arange(6).reshape(2,3), numpy.arange(0,12,2).reshape(2,3).T)\n",
      "            array([[ 10,  28],\n",
      "                   [ 28, 100]])\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.arange(8).reshape(2,4), numpy.arange(0,16,2).reshape(2,4), op=numpy.dot)\n",
      "            array([[ 28,  76],\n",
      "                   [ 76, 252]])\n",
      "\n",
      "            >>> numpy.dot(numpy.arange(8).reshape(2,4), numpy.arange(0,16,2).reshape(2,4).T)\n",
      "            array([[ 28,  76],\n",
      "                   [ 76, 252]])\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.eye(2).astype(bool),\n",
      "            ...                   numpy.ones((2,2), dtype=bool),\n",
      "            ...                   op=lambda _a, _b: (_a & _b).sum())\n",
      "            array([[ True,  True],\n",
      "                   [ True,  True]], dtype=bool)\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.eye(2).astype(bool),\n",
      "            ...                   numpy.ones((2,), dtype=bool),\n",
      "            ...                   op=lambda _a, _b: (_a & _b).sum())\n",
      "            array([[ True],\n",
      "                   [ True]], dtype=bool)\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.ones((2,), dtype=bool),\n",
      "            ...                   numpy.eye(2).astype(bool),\n",
      "            ...                   op=lambda _a, _b: (_a & _b).sum())\n",
      "            array([[ True,  True]], dtype=bool)\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.ones((2,), dtype=bool),\n",
      "            ...                    numpy.ones((2,), dtype=bool),\n",
      "            ...                    op=lambda _a, _b: (_a & _b).sum())\n",
      "            array([[ True]], dtype=bool)\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.eye(2).astype(bool),\n",
      "            ...                   numpy.zeros((2,2), dtype=bool),\n",
      "            ...                   op=lambda _a, _b: (_a & _b).sum())\n",
      "            array([[False, False],\n",
      "                   [False, False]], dtype=bool)\n",
      "\n",
      "            >>> matrix_reduced_op(numpy.eye(2).astype(bool),\n",
      "            ...                   numpy.zeros((2,2), dtype=bool),\n",
      "            ...                   op=lambda _a, _b: (_a | _b).sum())\n",
      "            array([[ True,  True],\n",
      "                   [ True,  True]], dtype=bool)\n",
      "    \"\"\"\n",
      "\n",
      "    if a.ndim == 1:\n",
      "        a = a[None]\n",
      "\n",
      "    if b.ndim == 1:\n",
      "        b = b[None]\n",
      "\n",
      "    assert(a.ndim == b.ndim)\n",
      "\n",
      "    for i in xrange(1, a.ndim):\n",
      "        assert(a.shape[i] == b.shape[i])\n",
      "\n",
      "    out = numpy.empty((len(a), len(b)), dtype=numpy.promote_types(a.dtype, b.dtype))\n",
      "\n",
      "    for i, j in itertools.product(xrange(out.shape[0]), xrange(out.shape[1])):\n",
      "        out[i, j] = op(a[i], b[j])\n",
      "\n",
      "    return(out)\n",
      "\n",
      "\n",
      "def masks_intersection(a, b):\n",
      "    \"\"\"\n",
      "        Generates a mask that contains the points share by both masks.\n",
      "\n",
      "        Args:\n",
      "            a(numpy.ndarray):      first mask.\n",
      "            b(numpy.ndarray):      second mask.\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):    a mask that is only True where both masks are.\n",
      "\n",
      "        Examples:\n",
      "            >>> masks_intersection(numpy.eye(2).astype(bool),\n",
      "            ...                    numpy.ones((2,2), dtype=bool))\n",
      "            array([[1, 1],\n",
      "                   [1, 1]], dtype=uint64)\n",
      "\n",
      "            >>> masks_intersection(numpy.eye(2).astype(bool),\n",
      "            ...                    numpy.ones((2,), dtype=bool))\n",
      "            array([[1],\n",
      "                   [1]], dtype=uint64)\n",
      "\n",
      "            >>> masks_intersection(numpy.ones((2,), dtype=bool),\n",
      "            ...                    numpy.eye(2).astype(bool))\n",
      "            array([[1, 1]], dtype=uint64)\n",
      "\n",
      "            >>> masks_intersection(numpy.ones((2,), dtype=bool),\n",
      "            ...                    numpy.ones((2,), dtype=bool))\n",
      "            array([[2]], dtype=uint64)\n",
      "\n",
      "            >>> masks_intersection(numpy.eye(2).astype(bool),\n",
      "            ...                    numpy.zeros((2,2), dtype=bool))\n",
      "            array([[0, 0],\n",
      "                   [0, 0]], dtype=uint64)\n",
      "\n",
      "            >>> (numpy.arange(6).reshape(2,3) % 2).astype(bool)\n",
      "            array([[False,  True, False],\n",
      "                   [ True, False,  True]], dtype=bool)\n",
      "\n",
      "            >>> numpy.arange(6).reshape(2,3) == 4\n",
      "            array([[False, False, False],\n",
      "                   [False,  True, False]], dtype=bool)\n",
      "\n",
      "            >>> masks_intersection((numpy.arange(6).reshape(2,3) % 2).astype(bool),\n",
      "            ...                    (numpy.arange(6).reshape(2,3) == 4))\n",
      "            array([[0, 1],\n",
      "                   [0, 0]], dtype=uint64)\n",
      "    \"\"\"\n",
      "\n",
      "    assert(issubclass(a.dtype.type, numpy.bool_))\n",
      "    assert(issubclass(b.dtype.type, numpy.bool_))\n",
      "\n",
      "    if a.ndim == 1:\n",
      "        a = a[None]\n",
      "\n",
      "    if b.ndim == 1:\n",
      "        b = b[None]\n",
      "\n",
      "    assert(a.ndim == b.ndim == 2)\n",
      "    assert(a.shape[1] == b.shape[1])\n",
      "\n",
      "    out = numpy.empty((len(a), len(b)), dtype=numpy.uint64)\n",
      "\n",
      "    for i, j in itertools.product(xrange(out.shape[0]), xrange(out.shape[1])):\n",
      "        out[i, j] = (a[i] & b[j]).sum()\n",
      "\n",
      "    return(out)\n",
      "\n",
      "\n",
      "def masks_union(a, b):\n",
      "    \"\"\"\n",
      "        A mask that contains point contained in either mask.\n",
      "\n",
      "        Args:\n",
      "            a(numpy.ndarray):      first mask.\n",
      "            b(numpy.ndarray):      second mask.\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):    a mask that is True where either mask is.\n",
      "\n",
      "        Examples:\n",
      "            >>> masks_union(numpy.eye(2).astype(bool),\n",
      "            ...             numpy.ones((2,2), dtype=bool))\n",
      "            array([[2, 2],\n",
      "                   [2, 2]], dtype=uint64)\n",
      "\n",
      "            >>> masks_union(numpy.eye(2).astype(bool),\n",
      "            ...                    numpy.ones((2,), dtype=bool))\n",
      "            array([[2],\n",
      "                   [2]], dtype=uint64)\n",
      "\n",
      "            >>> masks_union(numpy.ones((2,), dtype=bool),\n",
      "            ...                    numpy.eye(2).astype(bool))\n",
      "            array([[2, 2]], dtype=uint64)\n",
      "\n",
      "            >>> masks_union(numpy.ones((2,), dtype=bool),\n",
      "            ...                    numpy.ones((2,), dtype=bool))\n",
      "            array([[2]], dtype=uint64)\n",
      "\n",
      "            >>> masks_union(numpy.eye(2).astype(bool),\n",
      "            ...             numpy.zeros((2,2), dtype=bool))\n",
      "            array([[1, 1],\n",
      "                   [1, 1]], dtype=uint64)\n",
      "\n",
      "            >>> (numpy.arange(6).reshape(2,3) % 2).astype(bool)\n",
      "            array([[False,  True, False],\n",
      "                   [ True, False,  True]], dtype=bool)\n",
      "\n",
      "            >>> numpy.arange(6).reshape(2,3) == 4\n",
      "            array([[False, False, False],\n",
      "                   [False,  True, False]], dtype=bool)\n",
      "\n",
      "            >>> masks_union((numpy.arange(6).reshape(2,3) % 2).astype(bool),\n",
      "            ...             (numpy.arange(6).reshape(2,3) == 4))\n",
      "            array([[1, 1],\n",
      "                   [2, 3]], dtype=uint64)\n",
      "    \"\"\"\n",
      "\n",
      "    assert(issubclass(a.dtype.type, numpy.bool_))\n",
      "    assert(issubclass(b.dtype.type, numpy.bool_))\n",
      "\n",
      "    if a.ndim == 1:\n",
      "        a = a[None]\n",
      "\n",
      "    if b.ndim == 1:\n",
      "        b = b[None]\n",
      "\n",
      "    assert(a.ndim == b.ndim == 2)\n",
      "    assert(a.shape[1] == b.shape[1])\n",
      "\n",
      "    out = numpy.empty((len(a), len(b)), dtype=numpy.uint64)\n",
      "\n",
      "    for i, j in itertools.product(xrange(out.shape[0]), xrange(out.shape[1])):\n",
      "        out[i, j] = (a[i] | b[j]).sum()\n",
      "\n",
      "    return(out)\n",
      "\n",
      "\n",
      "def masks_overlap_normalized(a, b):\n",
      "    \"\"\"\n",
      "        The area of intersection of the masks divided by the area of their union.\n",
      "\n",
      "        Args:\n",
      "            a(numpy.ndarray):      first mask.\n",
      "            b(numpy.ndarray):      second mask.\n",
      "\n",
      "        Returns:\n",
      "            out(numpy.ndarray):    ratio of the areas of the masks' intersection and union.\n",
      "\n",
      "        Examples:\n",
      "            >>> masks_overlap_normalized(numpy.eye(2).astype(bool),\n",
      "            ...                          numpy.ones((2,2), dtype=bool))\n",
      "            array([[ 0.5,  0.5],\n",
      "                   [ 0.5,  0.5]])\n",
      "\n",
      "            >>> masks_overlap_normalized(numpy.eye(2).astype(bool),\n",
      "            ...                          numpy.ones((2,), dtype=bool))\n",
      "            array([[ 0.5],\n",
      "                   [ 0.5]])\n",
      "\n",
      "            >>> masks_overlap_normalized(numpy.ones((2,), dtype=bool),\n",
      "            ...                          numpy.eye(2).astype(bool))\n",
      "            array([[ 0.5,  0.5]])\n",
      "\n",
      "            >>> masks_overlap_normalized(numpy.ones((2,), dtype=bool),\n",
      "            ...                          numpy.ones((2,), dtype=bool))\n",
      "            array([[ 1.]])\n",
      "\n",
      "            >>> masks_overlap_normalized(numpy.eye(2).astype(bool),\n",
      "            ...                          numpy.zeros((2,2), dtype=bool))\n",
      "            array([[ 0.,  0.],\n",
      "                   [ 0.,  0.]])\n",
      "\n",
      "            >>> masks_overlap_normalized(numpy.zeros((2,2), dtype=bool),\n",
      "            ...                          numpy.zeros((2,2), dtype=bool))\n",
      "            array([[ 0.,  0.],\n",
      "                   [ 0.,  0.]])\n",
      "\n",
      "            >>> (numpy.arange(6).reshape(2,3) % 2).astype(bool)\n",
      "            array([[False,  True, False],\n",
      "                   [ True, False,  True]], dtype=bool)\n",
      "\n",
      "            >>> numpy.arange(6).reshape(2,3) == 4\n",
      "            array([[False, False, False],\n",
      "                   [False,  True, False]], dtype=bool)\n",
      "\n",
      "            >>> masks_overlap_normalized((numpy.arange(6).reshape(2,3) % 2).astype(bool),\n",
      "            ...                          (numpy.arange(6).reshape(2,3) == 4))\n",
      "            array([[ 0.,  1.],\n",
      "                   [ 0.,  0.]])\n",
      "    \"\"\"\n",
      "\n",
      "    assert(issubclass(a.dtype.type, numpy.bool_))\n",
      "    assert(issubclass(b.dtype.type, numpy.bool_))\n",
      "\n",
      "    if a.ndim == 1:\n",
      "        a = a[None]\n",
      "\n",
      "    if b.ndim == 1:\n",
      "        b = b[None]\n",
      "\n",
      "    assert(a.ndim == b.ndim == 2)\n",
      "    assert(a.shape[1] == b.shape[1])\n",
      "\n",
      "    out = masks_intersection(a, b).astype(float)\n",
      "    out /= masks_union(a, b)\n",
      "\n",
      "    out[numpy.isnan(out)] = 0\n",
      "\n",
      "    return(out)\n",
      "\n",
      "\n",
      "def generate_contour(a_image, separation_distance = 1.0, margin = 1.0):\n",
      "    \"\"\"\n",
      "        Takes an image and extracts labeled contours from the mask using some minimum distance from the mask edge\n",
      "        and some margin.\n",
      "\n",
      "        Args:\n",
      "            a_image(numpy.ndarray):            takes an image.\n",
      "            separation_distance(float):        a separation distance from the edge of the mask for the center of the contour.\n",
      "            margin(float):                     the width of contour.\n",
      "\n",
      "        Returns:\n",
      "            (numpy.ndarray):                   an array with the labeled contours.\n",
      "\n",
      "        Examples:\n",
      "            >>> a = numpy.array([[ True,  True, False],\n",
      "            ...                  [False, False, False],\n",
      "            ...                  [ True,  True,  True]], dtype=bool);\n",
      "            >>> generate_contour(a)\n",
      "            array([[ True,  True, False],\n",
      "                   [False, False, False],\n",
      "                   [ True,  True,  True]], dtype=bool)\n",
      "\n",
      "            >>> generate_contour(numpy.eye(3))\n",
      "            array([[ 1.,  0.,  0.],\n",
      "                   [ 0.,  1.,  0.],\n",
      "                   [ 0.,  0.,  1.]])\n",
      "\n",
      "            >>> a = numpy.array([[False, False,  True, False, False, False,  True],\n",
      "            ...                  [ True, False, False, False,  True, False, False],\n",
      "            ...                  [ True,  True, False,  True,  True, False,  True],\n",
      "            ...                  [ True, False, False,  True,  True, False, False],\n",
      "            ...                  [ True, False, False, False, False, False, False],\n",
      "            ...                  [False,  True, False, False, False, False,  True],\n",
      "            ...                  [False,  True,  True, False, False, False, False]], dtype=bool)\n",
      "            >>> generate_contour(a)\n",
      "            array([[False, False,  True, False, False, False,  True],\n",
      "                   [ True, False, False, False,  True, False, False],\n",
      "                   [ True,  True, False,  True,  True, False,  True],\n",
      "                   [ True, False, False,  True,  True, False, False],\n",
      "                   [ True, False, False, False, False, False, False],\n",
      "                   [False,  True, False, False, False, False,  True],\n",
      "                   [False,  True,  True, False, False, False, False]], dtype=bool)\n",
      "    \"\"\"\n",
      "\n",
      "    half_thickness = margin / 2\n",
      "\n",
      "    lower_threshold = separation_distance - half_thickness\n",
      "    upper_threshold = separation_distance + half_thickness\n",
      "\n",
      "    a_mask_transformed = scipy.ndimage.morphology.distance_transform_edt(a_image)\n",
      "\n",
      "    above_lower_threshold = (lower_threshold <= a_mask_transformed)\n",
      "    below_upper_threshold = (a_mask_transformed <= upper_threshold)\n",
      "\n",
      "    a_mask_transformed_thresholded = (above_lower_threshold & below_upper_threshold)\n",
      "\n",
      "    a_image_contours = a_image * a_mask_transformed_thresholded\n",
      "\n",
      "    return(a_image_contours)\n",
      "\n",
      "\n",
      "def generate_labeled_contours(a_mask, separation_distance = 1.0, margin = 1.0):\n",
      "    \"\"\"\n",
      "        Takes a bool mask and extracts labeled contours from the mask using some minimum distance from the mask edge\n",
      "        and some margin.\n",
      "\n",
      "        Args:\n",
      "            a_mask(numpy.ndarray):             takes a bool mask.\n",
      "            separation_distance(float):        a separation distance from the edge of the mask for the center of the contour.\n",
      "            margin(float):                     the width of contour.\n",
      "\n",
      "        Returns:\n",
      "            (numpy.ndarray):                   an array with the labeled contours.\n",
      "\n",
      "        Examples:\n",
      "            >>> a = numpy.array([[ True,  True, False],\n",
      "            ...                  [False, False, False],\n",
      "            ...                  [ True,  True,  True]], dtype=bool);\n",
      "            >>> generate_labeled_contours(a)\n",
      "            array([[1, 1, 0],\n",
      "                   [0, 0, 0],\n",
      "                   [2, 2, 2]], dtype=int32)\n",
      "\n",
      "            >>> generate_labeled_contours(numpy.eye(3))\n",
      "            array([[1, 0, 0],\n",
      "                   [0, 1, 0],\n",
      "                   [0, 0, 1]], dtype=int32)\n",
      "\n",
      "            >>> a = numpy.array([[False, False,  True, False, False, False,  True],\n",
      "            ...                  [ True, False, False, False,  True, False, False],\n",
      "            ...                  [ True,  True, False,  True,  True, False,  True],\n",
      "            ...                  [ True, False, False,  True,  True, False, False],\n",
      "            ...                  [ True, False, False, False, False, False, False],\n",
      "            ...                  [False,  True, False, False, False, False,  True],\n",
      "            ...                  [False,  True,  True, False, False, False, False]], dtype=bool);\n",
      "            >>> generate_labeled_contours(a)\n",
      "            array([[0, 0, 1, 0, 0, 0, 2],\n",
      "                   [3, 0, 0, 0, 4, 0, 0],\n",
      "                   [3, 3, 0, 4, 4, 0, 5],\n",
      "                   [3, 0, 0, 4, 4, 0, 0],\n",
      "                   [3, 0, 0, 0, 0, 0, 0],\n",
      "                   [0, 3, 0, 0, 0, 0, 6],\n",
      "                   [0, 3, 3, 0, 0, 0, 0]], dtype=int32)\n",
      "    \"\"\"\n",
      "\n",
      "    a_mask_contoured = generate_contour(a_mask, separation_distance = separation_distance, margin = margin)\n",
      "\n",
      "    a_mask_contoured_labeled = scipy.ndimage.label(a_mask_contoured, structure = numpy.ones( (3,) * a_mask.ndim ))[0]\n",
      "\n",
      "    return(a_mask_contoured_labeled)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## Convenience for using h5py in iPython.\n",
      "#h5py.enable_ipython_completer()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "raw_data_filename = \"../an271606_2014_12_11_fov_01002.h5\"\n",
      "our_results_filename = \"../an271606_2014_12_11_fov_01002_results.h5\"\n",
      "simons_results_filename = \"../an271606_2014_12_11_fov_01002_rois.mat\"\n",
      "\n",
      "\n",
      "raw_data_filename = os.path.abspath(raw_data_filename)\n",
      "our_results_filename = os.path.abspath(our_results_filename)\n",
      "simons_results_filename = os.path.abspath(simons_results_filename)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read raw data as needed\n",
      "\n",
      "im = None\n",
      "\n",
      "with hp.File(raw_data_filename, \"r\") as f0:\n",
      "    frames = len(f0[\"images\"])\n",
      "    frame_size = np.array(f0[\"images\"].shape[1:])\n",
      "    frame_diag_len = scipy.spatial.distance.cdist(np.zeros_like(frame_size[None]), frame_size[None]).max()\n",
      "    im = f0[\"images\"][()]\n",
      "    im_max = im.max(axis=0)\n",
      "    im_mean = im.mean(axis=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read our results\n",
      "\n",
      "m1 = None\n",
      "c1 = None\n",
      "\n",
      "with hp.File(our_results_filename, \"r\") as f1:\n",
      "    m1 = f1[\"neurons\"][\"mask\"]\n",
      "    c1 = np.array([np.array(_.nonzero()).mean(axis=1) for _ in m1])\n",
      "    #c1 = f1[\"neurons\"][\"centroid\"]\n",
      "\n",
      "    cm1 = np.zeros_like(m1)\n",
      "    for i in xrange(len(m1)):\n",
      "        cm1[i] = generate_contour(m1[i])\n",
      "\n",
      "m1f = m1.reshape(len(m1), -1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read the ground truth\n",
      "\n",
      "m2 = None\n",
      "c2 = None\n",
      "\n",
      "with hp.File(simons_results_filename, \"r\") as f2:\n",
      "    m2 = np.array([f2[_][()].squeeze() for _ in f2[\"roi/mask\"][()].squeeze()])\n",
      "    assert(((m2 == 0) | (m2 == 1)).all())\n",
      "    m2 = m2.astype(bool)\n",
      "    \n",
      "    c2 = np.array([f2[_][()].squeeze() for _ in f2[\"roi/centroid_xy\"][()].squeeze()])\n",
      "    \n",
      "    b2 = np.array([f2[_][()].squeeze() for _ in f2[\"roi/burst_rate\"][()].squeeze()])\n",
      "\n",
      "    cm2 = np.zeros_like(m2)\n",
      "    for i in xrange(len(m2)):\n",
      "        cm2[i] = generate_contour(m2[i])\n",
      "\n",
      "m2f = m2.reshape(len(m2), -1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compute the normalized overlap of the masks and the centroid distance\n",
      "\n",
      "metrics = {}\n",
      "metrics[\"Mask\"] = masks_overlap_normalized(m1f, m2f)\n",
      "metrics[\"Distance\"] = scipy.spatial.distance.cdist(c1, c2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## Make a graph showing the results\n",
      "\n",
      "#fig, (ax1, ax2) = plt.subplots(1, 2)\n",
      "\n",
      "#fig.canvas.set_window_title(os.path.splitext(raw_data_filename)[0] + \"_benchmark\")\n",
      "\n",
      "#ax1.set_title(\"Normalized Mask Overlap\")\n",
      "#im1 = ax1.matshow(metrics[\"Mask\"], vmin=0.0, vmax=1.0, cmap=mpl.cm.Reds)\n",
      "#plt.colorbar(im1, ax=ax1)\n",
      "\n",
      "#ax2.set_title(\"Centroid L_2 Distance\")\n",
      "#im2 = ax2.matshow(metrics[\"Distance\"], vmin=0.0, vmax=frame_diag_len, cmap=mpl.cm.Reds_r)\n",
      "#plt.colorbar(im2, ax=ax2)\n",
      "\n",
      "#plt.tight_layout()\n",
      "#plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make a graph showing the results\n",
      "\n",
      "sliders_max = {}\n",
      "sliders_max[\"Mask\"] = 1.0\n",
      "sliders_max[\"Distance\"] = metrics[\"Distance\"].max()\n",
      "\n",
      "sliders_val = {}\n",
      "sliders_val[\"Mask\"] = metrics[\"Mask\"].mean()\n",
      "sliders_val[\"Distance\"] = metrics[\"Distance\"].mean()\n",
      "\n",
      "accepted = {}\n",
      "accepted[\"Mask\"] = np.zeros(metrics[\"Mask\"].shape, dtype=bool)\n",
      "accepted[\"Distance\"] = np.zeros(metrics[\"Distance\"].shape, dtype=bool)\n",
      "\n",
      "threshold_criteria = {}\n",
      "threshold_criteria[\"Mask\"] = {\"include_below\" : False, \"is_closed\" : True}\n",
      "threshold_criteria[\"Distance\"] = {\"include_below\" : True, \"is_closed\" : True}\n",
      "\n",
      "\n",
      "def compute_mapping_statistics(name):\n",
      "    global metrics\n",
      "    global e\n",
      "    global mm\n",
      "    global accepted\n",
      "    global sliders_val\n",
      "    global threshold_criteria\n",
      "    \n",
      "    name = str(name)\n",
      "    \n",
      "    accepted[name] = threshold_metric(metrics[name], sliders_val[name], **threshold_criteria[name])\n",
      "    \n",
      "    found, missed, unused = compute_mapping_matches(accepted[name].T)\n",
      "    recall, precision = compute_mapping_relevance(accepted[name].T)\n",
      "    \n",
      "    e[0] = found\n",
      "    e[1] = missed\n",
      "    e[2] = unused\n",
      "    mm[0] = recall\n",
      "    mm[1] = precision\n",
      "\n",
      "\n",
      "fig = plt.figure()\n",
      "ax1 = plt.subplot2grid((9,10), (0,0), rowspan=7, colspan=4)\n",
      "ax2 = plt.subplot2grid((9,10), (0,4), rowspan=7, colspan=2)\n",
      "ax3 = plt.subplot2grid((9,10), (0,6), rowspan=7, colspan=2)\n",
      "ax4 = plt.subplot2grid((9,10), (0,8), rowspan=7, colspan=2)\n",
      "rax = plt.subplot2grid((9,10), (7,0), rowspan=2, colspan=4)\n",
      "sax = plt.subplot2grid((9,10), (7,4), rowspan=2, colspan=6)\n",
      "\n",
      "fig.canvas.set_window_title(os.path.splitext(raw_data_filename)[0] + \"_benchmark\")\n",
      "\n",
      "ax1.matshow(im_max, cmap=mpl.cm.gray_r)\n",
      "\n",
      "e = np.zeros((3,), dtype=int)\n",
      "mm = np.zeros((2,), dtype=float)\n",
      "compute_mapping_statistics(\"Mask\")\n",
      "\n",
      "bar_width = 0.5\n",
      "burst_bins = 10\n",
      "\n",
      "ax2.set_title(\"Matches\")\n",
      "ax2.set_ylim(0, max(accepted[\"Mask\"].shape))\n",
      "ax2.set_xlim((0.0, 2*bar_width*(len(e)-0.5)))\n",
      "ax2.bar(np.arange(len(e)), e, width=bar_width, color=[\"purple\", \"red\", \"orange\"])#, bins=len(e), normed=False)\n",
      "ax2.set_xticks(np.arange(len(e)) + bar_width/2.0, minor=False)\n",
      "ax2.set_xticklabels([\"TP\", \"FN\", \"FP\"])\n",
      "\n",
      "ax3.set_title(\"Relevance\")\n",
      "ax3.set_ylim(0.0, 1.0)\n",
      "ax3.set_xlim((0.0, 2*bar_width*(len(mm)-0.5)))\n",
      "ax3.bar(np.arange(len(mm)), mm, width=bar_width, color=[\"blue\", \"green\"])#, bins=len(mm), normed=False)\n",
      "ax3.set_xticks(np.arange(len(mm)) + bar_width/2.0, minor=False)\n",
      "ax3.set_xticklabels([\"Recall\", \"Precision\"])\n",
      "\n",
      "\n",
      "_, burst_bins, _ = ax4.hist(b2, burst_bins, color=[\"red\"])\n",
      "\n",
      "radio = wgt.RadioButtons(rax, (\"Mask\", \"Distance\"))\n",
      "slider = wgt.Slider(sax, \"Mask\", 0.0, sliders_max[\"Mask\"], sliders_val[\"Mask\"])\n",
      "\n",
      "\n",
      "def slider_moving(val, name):\n",
      "    global accepted\n",
      "    global cm1\n",
      "    global cm2\n",
      "    global e\n",
      "    global mm\n",
      "    global sliders_val\n",
      "    global bar_width\n",
      "    global ax1\n",
      "    global ax2\n",
      "    \n",
      "    name = str(name)\n",
      "    \n",
      "    if (name == \"Mask\") or (name == \"Distance\"):\n",
      "        sliders_val[name] = val\n",
      "        compute_mapping_statistics(name)\n",
      "        ax2.cla()\n",
      "        ax2.set_title(\"Matches\")\n",
      "        ax2.set_ylim(0, max(accepted[name].shape))\n",
      "        ax2.set_xlim((0.0, 2*bar_width*(len(e)-0.5)))\n",
      "        \n",
      "        ax1.cla()\n",
      "        \n",
      "        ax2.bar(np.arange(len(e)), e, width=bar_width, color=[\"purple\", \"red\", \"darkorange\"])#, bins=len(e), normed=False)\n",
      "        ax2.set_xticks(np.arange(len(e)) + bar_width/2.0, minor=False)\n",
      "        ax2.set_xticklabels([\"TP\", \"FN\", \"FP\"])\n",
      "\n",
      "        ax3.cla()\n",
      "        ax3.set_title(\"Relevance\")\n",
      "        ax3.set_ylim(0.0, 1.0)\n",
      "        ax3.set_xlim((0.0, 2*bar_width*(len(mm)-0.5)))\n",
      "        \n",
      "        ax3.bar(np.arange(len(mm)), mm, width=bar_width, color=[\"blue\", \"green\"])#, bins=len(mm), normed=False)\n",
      "        ax3.set_xticks(np.arange(len(mm)) + bar_width/2.0, minor=False)\n",
      "        ax3.set_xticklabels([\"Recall\", \"Precision\"])\n",
      "        \n",
      "        tp_cm1 = (accepted[name].max(axis=1)[:, None, None] * cm1).max(axis = 0)\n",
      "        m_tp_cm1 = np.ma.masked_array(tp_cm1, tp_cm1 == 0)\n",
      "        \n",
      "        fn_cm2 = (~(accepted[name].max(axis=0))[:, None, None] * cm2).max(axis = 0)\n",
      "        m_fn_cm2 = np.ma.masked_array(fn_cm2, fn_cm2 == 0)\n",
      "        \n",
      "        fp_cm1 = (~(accepted[name].max(axis=1))[:, None, None] * cm1).max(axis = 0)\n",
      "        m_fp_cm1 = np.ma.masked_array(fp_cm1, fp_cm1 == 0)\n",
      "        \n",
      "        m2_matches = accepted[name].max(axis=0)\n",
      "        fcm2 = (cm2).max(axis = 0)\n",
      "        mfcm2 = np.ma.masked_array(fcm2, fcm2 == 0)\n",
      "        \n",
      "        ax1.matshow(im_max, cmap=mpl.cm.gray_r)\n",
      "        ax1.matshow(mfcm2, cmap=mpl.colors.ListedColormap([\"blue\"], name='from_list', N=None))#, vmin=0, vmax=10)\n",
      "        ax1.matshow(m_tp_cm1, cmap=mpl.colors.ListedColormap([\"purple\"], name='from_list', N=None))#, vmin=0, vmax=10)\n",
      "        ax1.matshow(m_fn_cm2, cmap=mpl.colors.ListedColormap([\"red\"], name='from_list', N=None))#, vmin=0, vmax=10)\n",
      "        ax1.matshow(m_fp_cm1, cmap=mpl.colors.ListedColormap([\"darkorange\"], name='from_list', N=None))#, vmin=0, vmax=10)\n",
      "        \n",
      "        ylim = ax4.get_ylim()\n",
      "        xlim = ax4.get_xlim()\n",
      "        ax4.cla()\n",
      "        ax4.set_title(\"Frequency vs. Burst Rate\")\n",
      "        ax4.set_ylabel(\"Frequency\")\n",
      "        ax4.set_xlabel(\"Burst Rate\")\n",
      "        ax4.set_ylim(ylim)\n",
      "        ax4.set_xlim(xlim)\n",
      "        ax4.hist(np.array([b2 * m2_matches, b2 * ~m2_matches]).T, burst_bins, stacked=True, color=[\"green\", \"red\"])\n",
      "        \n",
      "        plt.draw()\n",
      "    else:\n",
      "        assert(False)\n",
      "\n",
      "cid = None\n",
      "def radio_change(slider_name):\n",
      "    global slider\n",
      "    global sliders_val\n",
      "    global sliders_max\n",
      "    global cid\n",
      "    global sax\n",
      "    \n",
      "    slider_name = str(slider_name)\n",
      "    \n",
      "    if (slider_name == \"Mask\") or (slider_name == \"Distance\"):\n",
      "        slider.disconnect(cid)\n",
      "        \n",
      "        slider.valmax = sliders_max[slider_name]\n",
      "        sax.set_xlim(slider.valmin, slider.valmax)\n",
      "        slider.set_val(sliders_val[slider_name])\n",
      "        slider.label.set_text(slider_name)\n",
      "        \n",
      "        cid = slider.on_changed(lambda val : slider_moving(val, name=slider_name))\n",
      "        \n",
      "        slider_moving(slider.val, name=slider_name)\n",
      "        plt.draw()\n",
      "    else:\n",
      "        assert(False)\n",
      "\n",
      "\n",
      "radio_change(\"Mask\")\n",
      "radio.on_clicked(radio_change)\n",
      "plt.draw()\n",
      "\n",
      "plt.tight_layout()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "WARNING:py.warnings:/xopt_ilastik/local/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib-1.4.2-py2.7-macosx-10.9-x86_64.egg/matplotlib/tight_layout.py:225: UserWarning: tight_layout : falling back to Agg renderer\n",
        "  warnings.warn(\"tight_layout : falling back to Agg renderer\")\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}