index.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reproducing Watts and Strogatz ideas\n",
    "\n",
    "We've all said at some point \"it's such a small world!\" when facing an unexpected connection&mdash;a friend or a relative in common&mdash;with someone you just met. It just feels like an incredible coincidence! In their paper, Watts and Strogatz try to quantify that.\n",
    "\n",
    "The authors show that some networked systems showcase what they call \"small-worlds\" features. These involve a short path length and high clustering coefficient. In plain terms, what this means is that some systems have their components very well connected to their close neighbors while not too far from distant ones. If you got completely lost with the jargon of the paper, [Veritasium and Neil deGrasse Tyson](https://www.youtube.com/watch?v=TcxZSmzPw8k) have a nice video to illustrate the core idea.\n",
    "\n",
    "Here, we will reproduce the main result from Watts and Strogatz. We will build a clustered network that we will progressively randomize to generate a gradient from \"structured\" to \"randomized\". Then, we will use network properties to show that there is a region in this gradient that shows \"small-world\" properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ring Lattice\n",
    "\n",
    "Let's start by building the following network that the authors call \"Ring Lattice\" in the paper:\n",
    "\n",
    "![](RingLattice.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this network, everyone is connected to their nearest neighbors. Can you imagine how this would look like in the form of an adjacency matrix? In this case, the network has N=20 nodes, and K=4 edges or links per node."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "<font color='green'>**Exercise 1:** Manually transform the following identify matrix into a Ring Lattice with N=5 and K=4:<font>\n",
    "\n",
    "```{r}\n",
    "rbind(\n",
    "        c(1,0,0,0,0),\n",
    "        c(0,1,0,0,0),\n",
    "        c(0,0,1,0,0),\n",
    "        c(0,0,0,1,0),\n",
    "        c(0,0,0,0,1)\n",
    "        )\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "```{r}\n",
    "rbind(\n",
    "        c(0,1,1,1,1),\n",
    "        c(1,0,1,1,1),\n",
    "        c(1,1,0,1,1),\n",
    "        c(1,1,1,0,1),\n",
    "        c(1,1,1,1,0)\n",
    "        )\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can actually automatically generate an adjacency matrix like this one for any N and K using the function below (don't worry too much about the details of the function). Try different combinations of N nodes and K edges and try to see if you understand the shape of the adjacency matrix. <font color='red'>Notice that N>K, and the number of edges K will always have to be an even number if we want to have a Ring Lattice!</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=10\n",
    "K=4\n",
    "ring_lattice <- function(N, K){\n",
    "    if(K%%2!=0){\n",
    "     stop(\"The number of edges K must be an even number\", call. = FALSE)\n",
    "    }\n",
    "    K <- K*0.5\n",
    "    return(toeplitz(c(0, rep(1, K), rep(0, N - 2 * K - 1), rep(1, K))))\n",
    "}\n",
    "\n",
    "print(ring_lattice(N=N,K=K))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the network\n",
    "\n",
    "There is a package called \"igraph\" that will help you visualize the network. The first thing we need to do is import the library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(igraph)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Then we need to convert our Ring lattice---which is represented as an adjacency matrix---to an igraph network. To do this, we can use the function 'graph_from_adjacency_matrix()', making sure we define the network as 'indirected' (there is no directionality in this network's links, as opposed to a food web where they represent the transfer of energy from one node to another). We can visualize this network using the 'plot()' function in R:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat <- ring_lattice(N,K)\n",
    "g <- graph_from_adjacency_matrix(mat, mode=\"undirected\")\n",
    "plot(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 2:** Try to reproduce the exact ring lattice shown in the Watts and Strogatz paper. You can tune the graph visualization using the plot options shown in the following [tutorial](https://kateto.net/netscix2016.html) (don't spend too much time on this, just give it a quick try; e.g. give some curvature to the links).<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert the code for the visualization here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 3:** Calculate the connectance for this graph.<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Insert the code for calculating the connectance here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "<font color='green'>_The connectance for the Watts and Strogatz ring lattice is 0.2_<font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Randomize ring lattice\n",
    "\n",
    "Randomizing the ring lattice it's conceptually very easy but a bit tricky to code. Don't worry about the details of how to code this. Instead, try to reread the caption of Figure 1 in the paper, and understand what are the conditions for this randomization process. What is the probability $p$ actually doing? What are the key steps of the process? Discuss this with the members of your group."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the following code to produce a ring lattice that has been randomized by a probability $p$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p <- 0.01\n",
    "watts_random <- function(N,K,p){\n",
    "  d <- diag(N)\n",
    "  mat <- ring_lattice(N,K) + d\n",
    "  for (i in 1:K){\n",
    "    for (j in 1:N){\n",
    "      if(sample(c(FALSE,TRUE),size = 1, prob = c(1-p,p))){\n",
    "        idx <- sample(1:(N-1),1)\n",
    "        k <- (i+j)%%(N+1)\n",
    "        k <- N*(k==0)+k\n",
    "        k_ <- (i+j+idx)%%(N+1)\n",
    "        k_ <- N*(k_==0)+k_\n",
    "        if(j==k){next}\n",
    "        if(mat[j, k_]==0){\n",
    "          mat[j,k] <- 0\n",
    "          mat[k,j] <- 0\n",
    "          mat[j,k_] <- 1\n",
    "          mat[k_,j] <- 1          \n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "  mat <- mat-diag(N)\n",
    "  return(mat)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 4:** Will a randomized graph with probability $p$, $N=20$ and $K=6$ have the same connectance as the ring lattice counterpart?<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert the code for doing this comparison here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 5:** Play around with the plot function to see how the graph changes when modifying the probability $p$<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g <- graph_from_adjacency_matrix(watts_random(N=15,K=4,p=0.1), mode=\"undirected\")\n",
    "plot(g, layout=layout_in_circle(g), edge.color=\"#525252\", edge.width=0.6,\n",
    "            vertex.color=\"#dfc27d\", edge.curved=0, vertex.label=NA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Calculate the clustering coefficient\n",
    "\n",
    "Watts and Strogatz define the clustering coefficient as follows:\n",
    "> The clustering coefficient $C(p)$ is defined as follows. Suppose that a vertex $v$ has $k_v$ neighbors; then at most $k_v(k_v-1)/2$ edges can exist between them (this occurs when every neighbor of $v$ is connected to every other neighbor of $v$). Let $C_v$ denote the fraction of these allowable edges that actually exist. Define $C$ as the average of $C_v$ over all $v$. \n",
    "\n",
    "We will try to write a piece of code that calculates that step by step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 6:** We will start by generating a fully random graph with the function 'watts_random()' for $N=20$, $K=4$, and $p=1$. For a given node $idx=1$, try to write a line of code that finds all of its neighbors. You can represent this as a Boolean array.<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Insert the line of code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "```{r}\n",
    "mat <- watts_random(20,4,1)\n",
    "idx <- 1\n",
    "print((mat[idx,]==1 | mat[,idx]==1))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 7:** So now we know what the neighbors are for a particular node $idx$. Find the way to generate the adjacency matrix $M$ for all neighbors of $idx$<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Insert the line of code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "```{r}\n",
    "mat <- watts_random(20,4,1)\n",
    "idx <- 1\n",
    "condition <- (mat[idx,]==1 | mat[,idx]==1)\n",
    "M <- mat[condition,condition]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 8:** Use this matrix to calculate the contribution of the node $idx$ to the clustering coefficient.<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Insert the line of code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "solution2": "hidden"
   },
   "source": [
    "```{r}\n",
    "contribution <- sum(M)/(nrow(M)**2-nrow(M))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 9:** Put all the steps together to build a function called 'clust_coeff()' that calculates the clustering coefficient for a given matrix 'mat'. Hint: you will have to loop across all nodes and average the different contributions to the clustering coefficient<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Insert the line of code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "clust_coeff <- function(mat){\n",
    "  cc <- 0\n",
    "  for (i in 1:nrow(mat)){\n",
    "    j <- (mat[i,]==1 | mat[,i]==1)\n",
    "    k <- sum(j)\n",
    "    cc <- cc + sum(mat[j,j])/(k**2-k)\n",
    "  }\n",
    "  return(cc/nrow(mat))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Calculate the Path length\n",
    "\n",
    "We won't do this step by step, but we will find online a function that does it for us."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<font color='green'>**Exercise 10:** Try to find online a function online that could do this for the same random graph that we used above. Hint: the \"igraph\" library that we used to visualize the network has plenty of functions to calculate different network properties, maybe you should try googling that.<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden",
    "solution2_first": true
   },
   "outputs": [],
   "source": [
    "# Insert the line of code here\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "# Since the matrix is a random matrix, your value and my value will most likely be different!\n",
    "mat <- watts_random(20,4,1)\n",
    "g <- graph_from_adjacency_matrix(mat, mode=\"undirected\")\n",
    "average.path.length(g, directed=FALSE, unconnected=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Reproducing Figure 2 from Watts and Strogatz \n",
    "\n",
    "To reproduce figure 2 let's first read what are their instructions:\n",
    "> The figure are averages over 20 random realizations of the rewiring process described in Fig.1, and have been normalized by the values $L(0)$, $C(0)$ for a regular lattice. All the graphs have $N=1000$ vertices and an average degree of $k=10$ edges per vertex. We note that a logarithmic horizontal scale has been used to resolve the rapid drop in $L(p)$, corresponding to the onset of the small-world phenomenon.\n",
    "\n",
    "The first thing we need to do is define the $x$ axis. This is the randomization probability p, and they choose a logarithmic scale to be able to visualize the \"small-world\" property"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x <- exp(seq(log(0.0001), log(1), length.out = 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now we can randomly generate matrices and calculate the path length and clustering coefficient. Make sure that you've compiled all functions necessary for this. This might take some time, be patient!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N <- 1000\n",
    "K <- 10\n",
    "replicates <- 20\n",
    "results <- c()\n",
    "for (i in x){\n",
    "  for (j in 1:replicates){\n",
    "    rmat <- watts_random(N,K,i)\n",
    "    g <- graph_from_adjacency_matrix(rmat, mode=\"undirected\")\n",
    "    results <- rbind(results, c(i,average.path.length(g, directed=FALSE, unconnected=FALSE), clust_coeff(rmat)))\n",
    "  }\n",
    "}\n",
    "results <- data.frame(p=results[,1], l=results[,2], c=results[,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now we just need to average across all probabilities so that for every $p$ we have the average clustering coefficient and path length across replicates the 20 replicates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat <- aggregate(results[,2:3], list(results$p), mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now we can calculate the quantities $L(p)/L(0)$ and $C(p)/C(0)$. To do so, we need to calculate L(0) and C(0), and normalize our dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat <- watts_random(N,K,0)\n",
    "g <- graph_from_adjacency_matrix(mat, mode=\"undirected\")\n",
    "\n",
    "l0 <- average.path.length(g, directed=FALSE, unconnected=FALSE)\n",
    "c0 <- clust_coeff(mat)\n",
    "\n",
    "dat$l <- dat$l/l0\n",
    "dat$c <- dat$c/c0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Finally, let's visualize the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(as.numeric(as.character(dat$Group.1)), dat$l, pch=16, type=\"p\", ylim=c(0,1), log=\"x\", ylab = \"\", xlab = \"p\")\n",
    "points(as.numeric(as.character(dat$Group.1)), dat$c, pch=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='purple'>Notice that there is a region (a particular range of probabilities), where the clustering coefficient is very high and the path length is very short. The networks falling within such a range of probabilities show small-world features (i.e. high clustering coefficient and short path lengths)!<font>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Investigate whether or not empirical ecological networks show \"small-world\" features \n",
    "\n",
    "Watts and Strogatz show that some real-world networks present small-world structures. However, they don't look at any ecological network. <font color='red'>**Everyone has to do this part individually in your own computer.** <font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, download a food web from a database such as [web of life](http://www.web-of-life.es), [globalwebdb](https://www.globalwebdb.com) or [icon.colorado.edu](https://icon.colorado.edu/#!/networks), and investigate $L$ and $C$. Here, I will show you an example, but you should download your own network and do the same analyses. For the sake of simplicity, we are going to consider all networks undirected.\n",
    "\n",
    "I will be using the food web with identifier \"FW_008\" from [web of life](http://www.web-of-life.es). Once you have downloaded the network, you can read it in you R session (you don't need to use this notebook for it). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat <- as.matrix(read.table(\"FW_008.csv\", sep=\",\"))\n",
    "\n",
    "#Since this network is weighted, I will convert the matrix to a binary form and remove any NA:\n",
    "mat <- 1*(mat>0)\n",
    "mat[is.na(mat)] <- 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the number of species and links:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n <- nrow(mat)\n",
    "k <- sum(mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's generate a random matrix with these parameters. You can do this as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmat <- matrix(sample(c(rep(1,k), rep(0,n*n-k))), n, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the path length and clustering coefficient for both \"mat\" and \"rmat\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Clustering coefficient\n",
    "cc_random <- clust_coeff(rmat)\n",
    "cc_empirical <- clust_coeff(mat)\n",
    "\n",
    "#Path length\n",
    "g_random <- graph_from_adjacency_matrix(rmat, mode=\"undirected\")\n",
    "g_empirical <- graph_from_adjacency_matrix(mat, mode=\"undirected\")\n",
    "pl_random <- average.path.length(g_random, directed=FALSE, unconnected=FALSE)\n",
    "pl_empirical <- average.path.length(g_empirical, directed=FALSE, unconnected=FALSE)\n",
    "\n",
    "print(data.frame(row.names =c(\"random\",\"empirical\"), cc = c(cc_random, cc_empirical), pl=c(pl_random, pl_empirical)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The network shows similar path length than its random counterpart but much higher clustering coefficient. It's such a small world!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill [this form](https://ee.kobotoolbox.org/x/#kfyNyy1V) to report the results you see regarding your food web. <font color='green'>**Run the code below to see everyone's results:**<font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(httr)\n",
    "rawdata<-httr::GET(\"https://kc.kobotoolbox.org/api/v1/data/422894\")\n",
    "form_metadata <- httr::content(rawdata)\n",
    "dat <- sapply(form_metadata, function(x) c(x$What_s_your_name, x$What_s_the_network_identifier,\n",
    "                                   x$What_s_the_random_path_length, x$What_s_the_measured_path_length,\n",
    "                                   x$What_s_the_random_clustering_coefficient, x$What_s_the_measured_ustering_coefficient))\n",
    "print(data.frame(who=dat[1,], network=dat[2,], random.pl=dat[3,], real.pl=dat[4,], random.cc=dat[5,], real.cc=dat[6,]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}