PersistentHomologyExample.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial for running a subset of the NYC foods from Open Food Facts (https://world.openfoodfacts.org/) through Eirene.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First we load in packages to the environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Recompiling stale cache file /opt/julia/compiled/v1.2/Eirene/mCiG1.ji for Eirene [9c0f25c4-2ca1-5870-89f6-52640788da1d]\n",
      "└ @ Base loading.jl:1240\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<script>\n",
       "// Immediately-invoked-function-expression to avoid global variables.\n",
       "(function() {\n",
       "    var warning_div = document.getElementById(\"webio-warning-1020910187679004178\");\n",
       "    var hide = function () {\n",
       "        var script = document.getElementById(\"webio-setup-9876527207510804383\");\n",
       "        var parent = script && script.parentElement;\n",
       "        var grandparent = parent && parent.parentElement;\n",
       "        if (grandparent) {\n",
       "            grandparent.style.display = \"none\";\n",
       "        }\n",
       "        warning_div.style.display = \"none\";\n",
       "    };\n",
       "    if (typeof Jupyter !== \"undefined\") {\n",
       "        console.log(\"WebIO detected Jupyter notebook environment.\");\n",
       "        // Jupyter notebook.\n",
       "        var extensions = (\n",
       "            Jupyter\n",
       "            && Jupyter.notebook.config.data\n",
       "            && Jupyter.notebook.config.data.load_extensions\n",
       "        );\n",
       "        if (extensions && extensions[\"webio-jupyter-notebook\"]) {\n",
       "            // Extension already loaded.\n",
       "            console.log(\"Jupyter WebIO nbextension detected; not loading ad-hoc.\");\n",
       "            hide();\n",
       "            return;\n",
       "        }\n",
       "    } else if (window.location.pathname.includes(\"/lab\")) {\n",
       "        // Guessing JupyterLa\n",
       "        console.log(\"Jupyter Lab detected; make sure the @webio/jupyter-lab-provider labextension is installed.\");\n",
       "        hide();\n",
       "        return;\n",
       "    }\n",
       "})();\n",
       "\n",
       "</script>\n",
       "<p\n",
       "    id=\"webio-warning-1020910187679004178\"\n",
       "    class=\"output_text output_stderr\"\n",
       "    style=\"padding: 1em; font-weight: bold;\"\n",
       ">\n",
       "    Unable to load WebIO. Please make sure WebIO works for your Jupyter client.\n",
       "    For troubleshooting, please see <a href=\"https://juliagizmos.github.io/WebIO.jl/latest/providers/ijulia/\">\n",
       "    the WebIO/IJulia documentation</a>.\n",
       "    <!-- TODO: link to installation docs. -->\n",
       "</p>\n"
      ],
      "text/plain": [
       "HTML{String}(\"<script>\\n// Immediately-invoked-function-expression to avoid global variables.\\n(function() {\\n    var warning_div = document.getElementById(\\\"webio-warning-1020910187679004178\\\");\\n    var hide = function () {\\n        var script = document.getElementById(\\\"webio-setup-9876527207510804383\\\");\\n        var parent = script && script.parentElement;\\n        var grandparent = parent && parent.parentElement;\\n        if (grandparent) {\\n            grandparent.style.display = \\\"none\\\";\\n        }\\n        warning_div.style.display = \\\"none\\\";\\n    };\\n    if (typeof Jupyter !== \\\"undefined\\\") {\\n        console.log(\\\"WebIO detected Jupyter notebook environment.\\\");\\n        // Jupyter notebook.\\n        var extensions = (\\n            Jupyter\\n            && Jupyter.notebook.config.data\\n            && Jupyter.notebook.config.data.load_extensions\\n        );\\n        if (extensions && extensions[\\\"webio-jupyter-notebook\\\"]) {\\n            // Extension already loaded.\\n            console.log(\\\"Jupyter WebIO nbextension detected; not loading ad-hoc.\\\");\\n            hide();\\n            return;\\n        }\\n    } else if (window.location.pathname.includes(\\\"/lab\\\")) {\\n        // Guessing JupyterLa\\n        console.log(\\\"Jupyter Lab detected; make sure the @webio/jupyter-lab-provider labextension is installed.\\\");\\n        hide();\\n        return;\\n    }\\n})();\\n\\n</script>\\n<p\\n    id=\\\"webio-warning-1020910187679004178\\\"\\n    class=\\\"output_text output_stderr\\\"\\n    style=\\\"padding: 1em; font-weight: bold;\\\"\\n>\\n    Unable to load WebIO. Please make sure WebIO works for your Jupyter client.\\n    For troubleshooting, please see <a href=\\\"https://juliagizmos.github.io/WebIO.jl/latest/providers/ijulia/\\\">\\n    the WebIO/IJulia documentation</a>.\\n    <!-- TODO: link to installation docs. -->\\n</p>\\n\")"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Precompiling Combinatorics [861a8166-3701-5b0c-9a16-15d98fcdc6aa]\n",
      "└ @ Base loading.jl:1242\n",
      "┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]\n",
      "└ @ Base loading.jl:1242\n"
     ]
    }
   ],
   "source": [
    "## Load packages\n",
    "\n",
    "using Pkg\n",
    "using JLD\n",
    "using Eirene\n",
    "using Combinatorics\n",
    "using SparseArrays\n",
    "using Plots\n",
    "theme(:dark)\n",
    "using LightGraphs\n",
    "using GraphPlot\n",
    "\n",
    "include(\"PH_helper_functions.jl\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next load in the NYC food data. \n",
    "\n",
    "We will see a list of food products (products_list), a list of ingredients used in the food (ingredients_list), the incidence matrix indicating which ingredients are in which food products (ingredientsByProducts), and the sugar content of each product. Note for means of this tutorial any product with missing sugar information was given a value of 500."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load in incidence matrix\n",
    "loaded_data = load(\"data/ingredientsByProducts.jld\")\n",
    "\n",
    "products_list = loaded_data[\"products_list\"]\n",
    "ingredients_list = loaded_data[\"ingredients_list\"]\n",
    "ingredientsByProducts = loaded_data[\"ingredientsByProducts\"]\n",
    "sugars = loaded_data[\"sugars\"]\n",
    "\n",
    "# Remove loaded_data to free up memory\n",
    "loaded_data = nothing\n",
    "GC.gc()\n",
    "\n",
    "\n",
    "# Quickly view the productsByIngredients incidence matrix\n",
    "gr()\n",
    "heatmap(ingredientsByProducts)\n",
    "title!(\"Producs by ingredients incidence\")\n",
    "xlabel!(\"Products\")\n",
    "ylabel!(\"Ingredients\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert incidence matrix to the proper format for Eirene\n",
    "\n",
    "Eirene can handle data in the form of an incidence matrix but requires it be converted into a simplex stream. Generally, we need to enumerate and weight each simplex within our final simplicial complex and store this in particular formats for Eirene. We will store the simplices themselves in the array of arrays E, which first has the nodes, then the edges, then the 2-simplices, and so on. Next we need to create the boundary operator D which tells us which simplices are faces of other simplices. Finally we take speciic projections of D and E (called rv, cp, and dv) together with a simplex-weight vector fv and pass these objects to eirene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## First we create lists of simplices of each dimension we care about\n",
    "\n",
    "nIngredients = size(ingredientsByProducts)[1]\n",
    "nProducts = size(ingredientsByProducts)[2]\n",
    "node_list = 1:nProducts\n",
    "\n",
    "\n",
    "edge_list, simp2_list, simp3_list = createSimplexLists(ingredientsByProducts,ingredients_list)\n",
    "# Now edge_list is an array of the edges, simp2_list an array of the 2-simplices, and so on.\n",
    "\n",
    "# Count the number of simplices we have in each dimension\n",
    "nNodes = length(node_list)\n",
    "nEdges = length(edge_list)\n",
    "nsimp2 = length(simp2_list)\n",
    "nsimp3 = length(simp3_list)\n",
    "\n",
    "println(\"Created lists of nodes, 1-simplices, 2-simplices, and 3-simplices\")\n",
    "\n",
    "\n",
    "# Combine these lists into the full cell list E\n",
    "println(\"Creating cell list... \")\n",
    "E = createCellList(nNodes,edge_list,simp2_list)\n",
    "nCells = length(E)\n",
    "\n",
    "println(\"Finished constructing cell list with $nNodes nodes, $nEdges edges, $nsimp2 2-simplices, and $nsimp3 3-simplices.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we compute the boundary matrix D. This matrix is composed of smaller boundary matrices that describe faces of simplices in a specific dimension. For example, the D_2 matrix of size nEdges x n2simp describes which edges make up each 2-simplex. The D_1 matrix is nNodes x nEdges and records which nodes parent each edge. The final boundary matrix D is an nCells x nCells matrix with D_1 and D_2 in the appropriate spots, and 0 elsewhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "println(\"Constructing D matrix...\")\n",
    "D = createDMatrix(E,nNodes,nEdges)\n",
    "println(\"Finished constructing boundary matrix D.\")\n",
    "\n",
    "# D is a sparse matrix due to its size, so we can take the following \n",
    "# slices of the sparse matrix which will be enough to recover all the \n",
    "# information that eirene needs.\n",
    "rv = D.rowval\n",
    "cp = D.colptr\n",
    "dv = length.(E).-1\n",
    "\n",
    "# Assign simplex weights based on the participating food which has \n",
    "# the largest sugar content. This is where we use the node weights\n",
    "# (sugar content) to weight connections.\n",
    "\n",
    "node_weights = sugars\n",
    "fv = determineCellWeight.(E)\n",
    "\n",
    "println(\"Ready for eirene!\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Eirene! \n",
    "\n",
    "We pass our constructed objects that describe the filtration (rv,cp,dv, and fv) to eirene which returns a list C of all the persistent homology information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = eirene(rv=rv,cp=cp,dv=dv,fv=fv)\n",
    "\n",
    "println(\"Finished running persistent homology!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vlsualize outputs\n",
    "\n",
    "### Barcodes and Betti curves\n",
    "First we wil look at the barcode output which shows each persistent cavity as a bar. Then we will plot the Betti curves, which records the number of cavities alive at each filtration value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate Betti_1 from the barcode output\n",
    "bettiCurve_1 = betticurve(C,dim=1)\n",
    "\n",
    "\n",
    "# Plot!\n",
    "p1a = plot(bettiCurve_1[:,1],bettiCurve_1[:,2], label = [\"B_1\"], linewidth = 2,\n",
    "    xlim = (0,100))\n",
    "xlabel!(\"Sugars\")\n",
    "ylabel!(\"Betti_1\")\n",
    "title!(\"Betti curve\")\n",
    "p1b = plotBarcode(barcode(C,dim=1),nNodes)\n",
    "title!(\"Barcode\")\n",
    "xlabel!(\"Sugars\")\n",
    "ylabel!(\"Persistent homology\")\n",
    "\n",
    "plot(p1a,p1b,layout = (2,1), framestyle = :box)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Which food to make?\n",
    "\n",
    "Now that we have discovered many cavities do indeed exist within the NYC food space, let's zero in on a few specific cavities so we can brainstorm exactly what sort of novel food combination is missing (and would be a big hit!).\n",
    "\n",
    "In this next section we'll look at specific representatives and do some food brainstorming!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick persistent cavity number\n",
    "class_number = 45\n",
    "\n",
    "# Extract nodes involved in the cavity and plot cycle\n",
    "rep_i = classrep(C, class = class_number, dim=1)\n",
    "nodes_in_cycle = unique(collect(Iterators.flatten(edge_list[rep_i])))\n",
    "food_in_cycle = products_list[nodes_in_cycle]\n",
    "\n",
    "g = SimpleGraph(nNodes)\n",
    "for edge in edge_list\n",
    "    add_edge!(g, edge[1], edge[2])\n",
    "end\n",
    "\n",
    "sub_graph1,vmap = induced_subgraph(g,nodes_in_cycle)\n",
    "display(gplot(sub_graph1, nodelabel = products_list[nodes_in_cycle]))\n",
    "\n",
    "for product in products_list[nodes_in_cycle]\n",
    "    println(product)\n",
    "end\n",
    "\n"
   ]
  },
  {
   "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": 4
}