SheavesExample.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory sheaves tutorial\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will run through a simple sheaf example using Micahel Robinson's [pysheaf](https://github.com/kb1dds/pysheaf) python package. The code is still under active development so we will demonstrate a simple example from Figure 3 of Blevins and Bassett 2020."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Precompiling PyCall [438e738f-606a-5dbb-bf0a-cddfbfd45ab0]\n",
      "└ @ Base loading.jl:1242\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "PyObject <module 'pysheaf' from '/opt/conda/lib/python3.7/site-packages/pysheaf/__init__.py'>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Import python pysheaf package.\n",
    "using PyCall\n",
    "ps = pyimport(\"pysheaf\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cell signaling example\n",
    "\n",
    "In the cell signaling example from Figure 3a, our system is set up\n",
    "\n",
    "(Cell 1) --------> (Cell 2) --------> (Cell 3)\n",
    "\n",
    "and the effect across edges is\n",
    "\n",
    "(Cell 1) --- 2x ---> (Cell 2) --- 0.5x ---> (Cell 3)\n",
    "\n",
    "so that cell 2 signals at twice the rate of cell 1, and cell 3 signals at half the rate of cell 2.\n",
    "\n",
    "We start to represent this system the nodes of our system with `AddCell`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "shf=ps.Sheaf()\n",
    "\n",
    "shf.AddCell(\"Cell_1\",ps.Cell(\"real\")) # First argument is the name of the Cell, the second describes the stalk\n",
    "shf.AddCell(\"Cell_2\",ps.Cell(\"real\"))\n",
    "shf.AddCell(\"Cell_3\",ps.Cell(\"real\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we have three nodes, one for each cell, and we link them together with restriction maps following the edges in our system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining functions for restriction maps\n",
    "function mult_by_two(x)\n",
    "    return 2.0*x\n",
    "end\n",
    "\n",
    "function mult_by_half(x)\n",
    "    return 0.5*x\n",
    "end\n",
    "\n",
    "\n",
    "# Add these relations to the sheaf\n",
    "shf.AddCoface(\"Cell_1\",\"Cell_2\",ps.Coface(\"real\",\"real\",mult_by_two))\n",
    "shf.AddCoface(\"Cell_2\",\"Cell_3\",ps.Coface(\"real\",\"real\",mult_by_half))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now in order to reproduce the figure, we can assign the rates for each cell:\n",
    "\n",
    "| Cell Number   | Rate (number of molecules per time)        |\n",
    "| ------------- |-------------------------------------------:|\n",
    "| Cell 1        | 4                                          |\n",
    "| Cell 2        | 8                                          |\n",
    "| Cell 3        | 4                                          |\n",
    "\n",
    "\n",
    "and we know these also satisfy our maps. We will add the data assignments with `SetDataAssignment`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "shf.GetCell(\"Cell_1\").SetDataAssignment(ps.Assignment(\"real\",4.0))\n",
    "shf.GetCell(\"Cell_2\").SetDataAssignment(ps.Assignment(\"real\",8.0))\n",
    "shf.GetCell(\"Cell_3\").SetDataAssignment(ps.Assignment(\"real\",4.0))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next few lines of code check the consistency radius of the sheaf. For details see slide 16 of this [presentation](http://www.appliedcategorytheory.org/wp-content/uploads/2018/03/Michael-Robinson-Sheaf-Methods-for-Inference.pdf)\n",
    "\n",
    "Briefly, if we pick one node and propagate that value along the sheaf following the restriction maps, we may get different values at each node than those assigned. The consistency radius is the maximum difference between the actual data assignment and an expected value via propagation, and sets a lower bound on the distance to the nearest global section [1](https://arxiv.org/abs/1603.01446)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The consistency radius is 0.0"
     ]
    }
   ],
   "source": [
    "shf.mPreventRedundantExtendedAssignments = false\n",
    "\n",
    "shf.MaximallyExtendCell(\"Cell_1\")   # Use the data on Cell 1 to decide what values at other nodes *should* be\n",
    "shf.MaximallyExtendCell(\"Cell_2\")\n",
    "shf.MaximallyExtendCell(\"Cell_3\")\n",
    "\n",
    "consistency_radius = shf.ComputeConsistencyRadius()\n",
    "\n",
    "print(\"The consistency radius is $consistency_radius\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New observations\n",
    "\n",
    "If we observe new data from the system, will the new data fit our sheaf?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The consistency radius is 0.0"
     ]
    }
   ],
   "source": [
    "# If we add drug A and then we observe the following rates:\n",
    "\n",
    "shf.GetCell(\"Cell_1\").SetDataAssignment(ps.Assignment(\"real\",2.0))\n",
    "shf.GetCell(\"Cell_2\").SetDataAssignment(ps.Assignment(\"real\",4.0))\n",
    "shf.GetCell(\"Cell_3\").SetDataAssignment(ps.Assignment(\"real\",2.0))\n",
    "\n",
    "# Check to see if the new data is consistent with our original sheaf\n",
    "\n",
    "shf.mPreventRedundantExtendedAssignments = false\n",
    "\n",
    "shf.MaximallyExtendCell(\"Cell_1\")\n",
    "shf.MaximallyExtendCell(\"Cell_2\")\n",
    "shf.MaximallyExtendCell(\"Cell_3\")\n",
    "\n",
    "consistency_radius_A = shf.ComputeConsistencyRadius()\n",
    "\n",
    "print(\"The consistency radius is $consistency_radius_A\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The consistency radius is 4.0"
     ]
    }
   ],
   "source": [
    "# If we add drug B and then we observe the following rates:\n",
    "\n",
    "shf.GetCell(\"Cell_1\").SetDataAssignment(ps.Assignment(\"real\",4.0))\n",
    "shf.GetCell(\"Cell_2\").SetDataAssignment(ps.Assignment(\"real\",4.0))\n",
    "shf.GetCell(\"Cell_3\").SetDataAssignment(ps.Assignment(\"real\",2.0))\n",
    "\n",
    "# Check to see if the new data is consistent with our original sheaf\n",
    "\n",
    "shf.mPreventRedundantExtendedAssignments = false\n",
    "\n",
    "shf.MaximallyExtendCell(\"Cell_1\")\n",
    "shf.MaximallyExtendCell(\"Cell_2\")\n",
    "shf.MaximallyExtendCell(\"Cell_3\")\n",
    "\n",
    "consistency_radius_B = shf.ComputeConsistencyRadius()\n",
    "\n",
    "print(\"The consistency radius is $consistency_radius_B\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see that our new data does not fit the original model!\n",
    "\n",
    "| Cell Number   | Rate original       | Rate after perturbation 1|Rate after perturbation 2|\n",
    "| ------------- |-------------------------------------------:|-------------------------------------------:|-------------------------------------------:|\n",
    "| Cell 1        | 4                                          | 2| 4\n",
    "| Cell 2        | 8                                          | 4 | 4 |\n",
    "| Cell 3        | 4                                          | 2 | 2 |\n",
    "|               |\n",
    "| Consistency radius | 0 | 0 | 2|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}