{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# PART I: Newtork Basics - Adapted from https://assemblingnetwork.wordpress.com/tutorials/network-basics/\n",
    "# the igraph library is good for network based analyses\n",
    "# Note: NetworkX is a python library for working with networks\n",
    "install.packages(\"igraph\")\n",
    "library(igraph)\n",
    "library(NetIndices)\n",
    "\n",
    "# Build a network using basic formulas:\n",
    "graph.onelink <- graph.formula(A-+B)\n",
    " \n",
    "# This gives us a two species network (A and B) with one link (represented by A-+B). \n",
    "# With this function the (+) sign signifies the \"arrowhead\".\n",
    "plot.igraph(graph.onelink)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Using graph.formula() you have to write out all edges explicitly\n",
    "graph.fournode<-graph.formula(A-+C,B-+C,C-+D)\n",
    "plot.igraph(graph.fournode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#igraph has a function for generating random networks of varying size and connectance. \n",
    "graph.random.gnp<-erdos.renyi.game(n=20,p.or.m=.5,type=\"gnp\",directed=T)\n",
    "plot.igraph(graph.random.gnp)\n",
    " \n",
    "# We can also change the layout of the graph, here we will plot the nodes in a circle\n",
    "plot.igraph(graph.random.gnp,layout=layout.circle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The above random directed graph has 20 species (\"n\")\n",
    "# type='gnp' tells the function to assign links with the probability \"p\" that we sepecifiy \n",
    "# In this case, any two nodes have a 50% probability of being connected.\n",
    "\n",
    "# Can also set the number of links in the system to a value \"m\"\n",
    "graph.random.gnm<-erdos.renyi.game(n=20,p.or.m=100,type=\"gnm\",directed=T)\n",
    "plot.igraph(graph.random.gnm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Rather than being truly random, most real networks exhibit some type of organization. \n",
    "# Many real world networks are scale-free (i.e. the degree distribution follows a power law - \n",
    "# few nodes have many edges, and many nodes have few edges)\n",
    "# To model scale free networks Barabasi and Albert developed the preferential attachment model in 1999. In this model new nodes are more likely to link to nodes with a higher number of links.\n",
    " \n",
    "# In igraph we can use the barabasi.game() function:\n",
    "graph.barabasi.1<-barabasi.game(n=50,power=1)\n",
    " \n",
    "# For this graph I will introduce some new plotting tools to specify layout and vertex/edge properties.\n",
    " \n",
    "plot.igraph(graph.barabasi.1,\n",
    "layout=layout.fruchterman.reingold,\n",
    "vertex.size=10,         # sets size of the vertex, default is 15\n",
    "vertex.label.cex=.5,    # size of the vertex label\n",
    "edge.arrow.size=.5        # sets size of the arrow at the end of the edge\n",
    ")\n",
    " \n",
    "# there are a number of different plotting parameters see\n",
    "#?igraph.plotting\n",
    "#for details\n",
    " \n",
    "plot.igraph(graph.barabasi.1,\n",
    "layout=layout.fruchterman.reingold,\n",
    "vertex.size=10,\n",
    "vertex.label.cex=.5,\n",
    "edge.arrow.size=.5,\n",
    "mark.groups=list(c(1,7,4,13,10,16,15,41,42,29),\n",
    "c(2,48,5,36,43,33,9)), # draws polygon around nodes\n",
    "mark.col=c(\"green\",\"blue\")\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# In the above plot a green and blue polygon are used to highlight nodes 1 and 2 as hubs. The \"mark.groups\" argument allows you to draw a polygon of specified color (\"mark.col\") around the nodes you specify in a list. Because barabasi.game() will give you a different graph each time, the groups I have recorded will not be the same each time.\n",
    " \n",
    "# We can use a community detection algorithm to determine the most densely connected nodes in a graph.\n",
    " \n",
    "barabasi.community<-walktrap.community(graph.barabasi.1)\n",
    " \n",
    "# This algorithm uses random walks to find the most densely connected subgraphs.\n",
    " \n",
    "members<-membership(barabasi.community)\n",
    "# The members() function picks out the membership vector (list of nodes in the most densely connected subgraph) from the communtiy object (e.g., walktrap community).\n",
    " \n",
    "par(mar=c(.1,.1,.1,.1))    # sets the edges of the plotting area\n",
    "plot.igraph(graph.barabasi.1,\n",
    "layout=layout.fruchterman.reingold,\n",
    "vertex.size=10,\n",
    "vertex.label.cex=.5,\n",
    "edge.arrow.size=.5,\n",
    "mark.groups=list(members),\n",
    "mark.col=\"green\"\n",
    ")\n",
    "#barabasicomm1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# With the above plot the group with the green polygon surrounding it is the nodes listed as being a part of the walktrap community.\n",
    " \n",
    "# Now we will play around with the \"power\" argument to see how that impacts the graphs.\n",
    "# We will generate 4 networks with preferential attachment at varying levels.\n",
    "barabasi.game.2<-barabasi.game(n=50,power=.75)\n",
    "barabasi.game.3<-barabasi.game(n=50,power=.5)\n",
    "barabasi.game.4<-barabasi.game(n=50,power=.25)\n",
    "barabasi.game.5<-barabasi.game(n=50,power=0)\n",
    " \n",
    "# These can be organized into a list for convenience.\n",
    "barabasi.graphs<-list(barabasi.game.2,barabasi.game.3,barabasi.game.4,barabasi.game.5)\n",
    " \n",
    "# Now lets use community detection, this time with the walktrap algorithm.\n",
    "bg.community.list<-lapply(barabasi.graphs,walktrap.community)\n",
    "bg.membership.list<-lapply(bg.community.list,membership)\n",
    " \n",
    "txt<-c(\"a\",\"b\",\"c\",\"d\")    # vector for labeling the graphs\n",
    " \n",
    "# Plot these four graphs in one window with:\n",
    "par(mfrow=c(2,2),mar=c(.2,.2,.2,.2))\n",
    "# The for loop here plots each graph in the list one by one into the window prepared by par.\n",
    "for(i in 1:4){\n",
    "plot.igraph(barabasi.graphs[[i]],\n",
    "layout=layout.fruchterman.reingold,\n",
    "vertex.size=10,\n",
    "vertex.label.cex=.5,\n",
    "edge.arrow.size=.5,\n",
    "mark.groups=list(bg.membership.list[[i]]),\n",
    "mark.col=\"green\",\n",
    "frame=T # the frame argument plots a box around the graph\n",
    ")\n",
    "text(1,1,txt[i]) # calls from the vector to label the graph, adds to the graph that was last plotted\n",
    "}\n",
    " \n",
    "# Later we will look at the properties of these graphs to see exactly how they are different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Graphs can also be visualized as an adjacency matrix, usually an NxN matrix of 0s and 1s, \n",
    "# where 1 indicates an interaction, and 0 is no interaction.\n",
    " \n",
    "# The function get.adjacency() converts a graph into a matrix.\n",
    "barabasi.adjacency<-get.adjacency(graph.barabasi.1)\n",
    " \n",
    "# A number of other functions will use the adjacency matrix to calculate different network properties.\n",
    " \n",
    "# And we can get a graph from the adjacency matrix with graph.adjacency(). Often the data for graphs is in matrix form. Using graph.adjacency allows you to convert adjacency matrices into graph objects (that can then be plotted)\n",
    "graph.adjacency(barabasi.adjacency)\n",
    " \n",
    "# Similarly we can arrange the information in a number of other ways\n",
    "# Edgelist\n",
    "barabasi.edgelist<-get.edgelist(graph.barabasi.1)\n",
    " \n",
    "# Dataframe\n",
    "barabasi.data.frame<-get.data.frame(graph.barabasi.1,what=\"edges\")\n",
    "\n",
    "head(barabasi.edgelist)\n",
    "head(barabasi.data.frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# PART II: NETWORK PROPERTIES\n",
    "\n",
    "# Create a network to analyze using the preferential attachment model (power=.5)\n",
    " \n",
    "test.graph<-barabasi.game(100,power=0.5,m=2)\n",
    "# In this case we have set m=2, meaning that\n",
    "# for each new node 2 new links are created\n",
    " \n",
    "par(mar=c(.1,.1,.1,.1))\n",
    "plot.igraph(test.graph,\n",
    "layout=layout.fruchterman.reingold,\n",
    "vertex.size=7,\n",
    "vertex.label.cex=.5,\n",
    "edge.arrow.size=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# How large is the network (I know I set this when we made the network,\n",
    "# but what if I had not?)\n",
    " \n",
    "test.graph      # Tells me that it is an IGRAPH object with 100 nodes and 197 links,\n",
    "                # made with the Barabasi algorithm\n",
    "V(test.graph)   # gives the vertex sequence\n",
    "E(test.graph)   # gives the edge sequence (edge list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The \"GenInd()\" function from the NetIndices package takes an adjacency matrix as input\n",
    "test.graph.adj<-get.adjacency(test.graph,sparse=F)\n",
    "# in older versions of igraph the default was sparse=F,\n",
    "# but now you must specify, other wise you get a matrix of 1s and .s\n",
    " \n",
    "test.graph.properties<-GenInd(test.graph.adj)\n",
    " \n",
    "# The function output consists of 10 network properties.\n",
    "# I will consider five of them here:\n",
    " \n",
    "test.graph.properties$N            #number of nodes\n",
    " \n",
    "test.graph.properties$Ltot        #number of links\n",
    " \n",
    "test.graph.properties$LD        #link density (average # of links per node)\n",
    " \n",
    "test.graph.properties$C            #the connectance of the graph\n",
    "# This function measures connectance as L/(N*(N-1)) where L is links, and N is nodes\n",
    "# Connectance can also be calculated as L/(N^2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The degree of a node refers to the number of links associated with a node.\n",
    "# Degree can be measured as the links going in (\"in degree\"), out (\"out degree\"), or both.\n",
    "# The degree() function takes a graph input and gives the degree of specified nodes.\n",
    "# With the argument \"v=V(graph)\" you tell the function to give the degree of all nodes in the graph,\n",
    "# while the \"mode\" argument specifies in, out, or both.\n",
    " \n",
    "in.deg.testgraph<-degree(test.graph,v=V(test.graph),mode=\"in\")\n",
    "out.deg.testgraph<-degree(test.graph,v=V(test.graph),mode=\"out\")\n",
    "all.deg.testgraph<-degree(test.graph,v=V(test.graph),mode=\"all\")\n",
    " \n",
    "# Degree distribution is the cumulative frequency of nodes with a given degree\n",
    "# this, like degree() can be specified as \"in\", \"out\", or \"all\"\n",
    "deg.distr<-degree.distribution(test.graph,cumulative=T,mode=\"all\")\n",
    " \n",
    "# Using the power.law.fit() function I can fit a power law to the degree distribution\n",
    "power<-power.law.fit(all.deg.testgraph)\n",
    " \n",
    "# The output of the power.law.fit() function gives the exponent of the power law ($alpha)\n",
    "# and the log-likelihood of the parameters used to fit the power law distribution ($logLik)\n",
    "# Also, it performs a Kolmogov-Smirnov test to determine whether the given degree distribution \n",
    "# could have been drawn from the fitted power law distribution.\n",
    "# The function returns a test statistic ($KS.stat) and p-value ($KS.p) for that test\n",
    " \n",
    "# Then I can plot the degree distribution\n",
    "plot(deg.distr,log=\"xy\",\n",
    "ylim=c(.01,10),\n",
    "bg=\"black\",pch=21,\n",
    "xlab=\"Degree\",\n",
    "ylab=\"Cumulative Frequency\")\n",
    " \n",
    "# And the expected power law distribution\n",
    "lines(1:20,10*(1:20)^((-power$alpha)+1))\n",
    " \n",
    "# Graphs typically have a Poisson distribution (if they are random),\n",
    "# power law (preferential attachment), or truncated power law (many real networks) degree distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Diameter is essentially the longest path between two vertices\n",
    "diameter(test.graph) # Gives the length of the diameter while\n",
    "nodes.diameter<-get.diameter(test.graph) # Gives the labels for each node that participate in the diameter\n",
    " \n",
    "# To visualize the diameter\n",
    "# First define the node and edge attributes\n",
    "V(test.graph)$color<-\"skyblue\" # node default color\n",
    "V(test.graph)$size<-7 # node default size\n",
    "V(test.graph)[nodes.diameter]$color<-\"darkgreen\" # diameter node color\n",
    "V(test.graph)[nodes.diameter]$size<-10 # diameter node size\n",
    "V(test.graph)[nodes.diameter]$label.color<-\"white\" # diameter node label color\n",
    "\n",
    "E(test.graph)$color<-\"grey\" # all non-diameter edges will be grey\n",
    "E(test.graph,path=nodes.diameter)$color<-\"darkgreen\" # diameter edges will be dark green\n",
    "E(test.graph,path=nodes.diameter)$width<-2 # diameter edges will be wider\n",
    " \n",
    "# If you do not set the attributes of all of the nodes and edges then it will\n",
    "# default such that you only see what you have defined\n",
    "\n",
    "par(mar=c(.1,.1,.1,.1))\n",
    "plot.igraph(test.graph,\n",
    "layout=layout.fruchterman.reingold,\n",
    "vertex.label.cex=.5,\n",
    "edge.arrow.size=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Clustering coefficient is the proportion of\n",
    "# a nodes neighbors that can be reached by other neighbors\n",
    "# in igraph this property is apparently called \"transitivity\"\n",
    " \n",
    "transitivity(test.graph)\n",
    "# gives the clustering coefficient of the whole network\n",
    " \n",
    "head(transitivity(test.graph,type=\"local\"))\n",
    "# gives the clustering coefficient of each node\n",
    " \n",
    "# Betweenness is the number of shortest paths between two nodes that go through each node of interest\n",
    "graph.betweenness<-betweenness(test.graph,v=V(test.graph))\n",
    "head(graph.betweenness)\n",
    "\n",
    "# The edge betweenness score of an edge measures the number of shortest paths through it\n",
    "graph.edge.betweenness<-edge.betweenness(test.graph,e=E(test.graph))\n",
    "head(graph.edge.betweenness)\n",
    " \n",
    "# Closeness refers to how connected a node is to its neighbors\n",
    "graph.closeness<-closeness(test.graph,vids=V(test.graph))\n",
    "head(graph.closeness)\n",
    "\n",
    "# Clustering coefficient, betweenness, and closeness\n",
    "# all describe the small world properties of the network.\n",
    "# A network with small world properties is one in which\n",
    "# it takes a relatively short path to get from one node to the next\n",
    "# (e.g., six degrees of separation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Every graph can be decomposed into its component n-node subgraphs.\n",
    "# In particular there are 13 unique ways to arrange 3 nodes in directed graphs.\n",
    "# Here are the adjacency matrices for each of the 13 subgraphs\n",
    "s1<-matrix(c(0,1,0,0,0,1,0,0,0),nrow=3,ncol=3)\n",
    "s2<-matrix(c(0,1,1,0,0,1,0,0,0),nrow=3,ncol=3)\n",
    "s3<-matrix(c(0,1,0,0,0,1,1,0,0),nrow=3,ncol=3)\n",
    "s4<-matrix(c(0,0,1,0,0,1,0,0,0),nrow=3,ncol=3)\n",
    "s5<-matrix(c(0,1,1,0,0,0,0,0,0),nrow=3,ncol=3)\n",
    "d2<-matrix(c(0,1,1,1,0,1,0,0,0),nrow=3,ncol=3)\n",
    "d1<-matrix(c(0,1,1,0,0,1,0,1,0),nrow=3,ncol=3)\n",
    "d3<-matrix(c(0,0,1,1,0,0,1,0,0),nrow=3,ncol=3)\n",
    "d4<-matrix(c(0,0,0,1,0,1,0,1,0),nrow=3,ncol=3)\n",
    "d5<-matrix(c(0,1,1,0,0,1,1,0,0),nrow=3,ncol=3)\n",
    "d6<-matrix(c(0,1,1,1,0,1,1,1,0),nrow=3,ncol=3)\n",
    "d7<-matrix(c(0,1,1,1,0,1,1,0,0),nrow=3,ncol=3)\n",
    "d8<-matrix(c(0,1,1,1,0,0,1,0,0),nrow=3,ncol=3)\n",
    " \n",
    "# Make the 13 matrices into a list\n",
    "subgraph3.mat<-list(s1,s2,s3,s4,s5,d1,d2,d3,d4,d5,d6,d7,d8)\n",
    "# And convert the matrices into graph objects\n",
    "subgraph3.graph<-lapply(subgraph3.mat,graph.adjacency)\n",
    " \n",
    "# Loop through the list of subgraphs and count how many times that subgraph appears in the larger test.graph\n",
    "subgraph.count<-c()\n",
    "for(i in 1:13){\n",
    "subgraph.count[i]<- graph.count.subisomorphisms.vf2(test.graph,subgraph3.graph[[i]],vertex.color1=NULL,vertex.color2=NULL,edge.color1=NULL,edge.color2=NULL)\n",
    "}\n",
    "subgraph.count\n",
    "\n",
    "par(mar=c(0,0,0,0),mfrow=c(5,3))\n",
    "plot.igraph(subgraph3.graph[[1]])\n",
    "plot.igraph(subgraph3.graph[[2]])\n",
    "plot.igraph(subgraph3.graph[[3]])\n",
    "plot.igraph(subgraph3.graph[[4]])\n",
    "plot.igraph(subgraph3.graph[[5]])\n",
    "plot.igraph(subgraph3.graph[[6]])\n",
    "plot.igraph(subgraph3.graph[[7]])\n",
    "plot.igraph(subgraph3.graph[[8]])\n",
    "plot.igraph(subgraph3.graph[[9]])\n",
    "plot.igraph(subgraph3.graph[[10]])\n",
    "plot.igraph(subgraph3.graph[[11]])\n",
    "plot.igraph(subgraph3.graph[[12]])\n",
    "plot.igraph(subgraph3.graph[[13]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# PART III: FOOD WEB data set\n",
    "\n",
    "# Kim N. Mouritsen, Robert Poulin, John P. McLaughlin and David W. Thieltges. 2011.\n",
    "# Food web including metazoan parasites for an intertidal ecosystem in New Zealand.\n",
    "# Ecology 92:2006.\n",
    " \n",
    "# Website: http://esapubs.org/archive/ecol/E092/173/\n",
    " \n",
    "# Otago Harbour: intertidal mudflat\n",
    "otago.links.data<-read.csv(\"~/Dropbox/Teaching/2016/BMS262/Module2/Exercises/All_Otago_Data_Files/Otago_Data_Links.csv\")\n",
    "otago.nodes.data<-read.csv(\"~/Dropbox/Teaching/2016/BMS262/Module2/Exercises/All_Otago_Data_Files/Otago_Data_Nodes.csv\")\n",
    " \n",
    "# Column names for data\n",
    "colnames(otago.links.data)\n",
    "colnames(otago.nodes.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Convert the data into a graph object using the first 2 columns of the dataset as an edgelist\n",
    "otago.graph<-graph.edgelist(as.matrix(otago.links.data[,1:2]))\n",
    "# Create graph object of just predator prey links\n",
    "otago.graph.p<-graph.edgelist(as.matrix(otago.links.data[1:1206,1:2]))\n",
    " \n",
    "# Get the web into matrix form\n",
    "otago.adjmatrix<-get.adjacency(otago.graph,sparse=F)\n",
    "otago.adjmatrix.p<-get.adjacency(otago.graph.p,sparse=F)\n",
    " \n",
    "# Get the basic network indices from the matrices with GenInd()\n",
    "ind.otago<-GenInd(otago.adjmatrix)\n",
    "ind.otago.p<-GenInd(otago.adjmatrix.p)\n",
    " \n",
    "# Now to plot these two webs to get a feel for what we are dealing with\n",
    " \n",
    "par(mar=c(.1,.1,.1,.1))\n",
    "plot.igraph(otago.graph,vertex.label=NA,vertex.size=3,edge.arrow.size=.25,layout=layout.circle)\n",
    "plot.igraph(otago.graph.p,vertex.label=NA,vertex.size=3,edge.arrow.size=.25,layout=layout.circle)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The NetIndices package also has a function to get some of the trophic properties of the food web\n",
    "# TrophInd() takes in an adjacency matrix and gives an output of the trophic level of each node,\n",
    "# as well as an index of the degree of omnivory for each node\n",
    " \n",
    "troph.otago<-TrophInd(otago.adjmatrix)\n",
    "troph.otago.p<-TrophInd(otago.adjmatrix.p)\n",
    " \n",
    "# An interesting aside, by adding parasites to the web it increases the trophic level of all species in\n",
    "# this web.\n",
    " \n",
    "plot(troph.otago[1:123,1]~troph.otago.p[,1],xlab=\"Level Without Parasites\",ylab=\"Level With Parasites\")\n",
    "abline(a=0,b=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# An interesting use for this trophic level function is to then use trophic level as a plotting parameter.\n",
    "# This way, you can plot the food web nodes according to trophic height. This adds greatly to a plot\n",
    "# of a food web, since you can gain more information about the trophic structure of the web by simply\n",
    "# glancing at the plot.\n",
    " \n",
    "# First we need to create a two-column matrix identifying the x and y values for each node.\n",
    "layout.matrix.1<-matrix(nrow=length(V(otago.graph)), ncol=2) # Rows equal to the number of vertices\n",
    "layout.matrix.1[,1]<-runif(length(V(otago.graph))) # randomly assign along x-axis\n",
    "layout.matrix.1[,2]<-troph.otago$TL # y-axis value based on trophic level\n",
    " \n",
    "layout.matrix.1p<-matrix(nrow=length(V(otago.graph.p)), ncol=2) # Rows equal to the number of vertices\n",
    "layout.matrix.1p[,1]<-runif(length(V(otago.graph.p)))\n",
    "layout.matrix.1p[,2]<-troph.otago.p$TL\n",
    " \n",
    "# Now use these matrices to define the layout instead of using the circle layout\n",
    "par(mar=c(.1,.1,.1,.1),mfrow=c(1,2))\n",
    "options(repr.plot.width=10, repr.plot.height=12)\n",
    " \n",
    "plot.igraph(otago.graph,\n",
    "vertex.label.cex=.35,\n",
    "vertex.size=3,\n",
    "edge.arrow.size=.25,\n",
    "layout=layout.matrix.1)\n",
    " \n",
    "plot.igraph(otago.graph.p,\n",
    "vertex.label.cex=.35,\n",
    "vertex.size=3,\n",
    "edge.arrow.size=.25,\n",
    "layout=layout.matrix.1p)\n",
    " \n",
    "# Note: Using runif() means that there is some chance that two nodes with the same trophic level\n",
    "# will be right on top of one another"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The inclusion of parasites also impacts community detection\n",
    "wtc.otago<-walktrap.community(otago.graph)\n",
    "wtc.otago.p<-walktrap.community(otago.graph.p)\n",
    " \n",
    "par(mar=c(0,0,0,0),mfrow=c(1,2))\n",
    "#par(mar=c(.1,.1,.1,.1),mfrow=c(1,2))\n",
    "\n",
    "options(repr.plot.width=10, repr.plot.height=12)\n",
    "\n",
    "plot.igraph(otago.graph,\n",
    "vertex.label.cex=.35,\n",
    "vertex.size=3,\n",
    "edge.arrow.size=.25,\n",
    "layout=layout.matrix.1,\n",
    "mark.groups=wtc.otago$membership,\n",
    "mark.col=\"green\")\n",
    " \n",
    "plot.igraph(otago.graph.p,\n",
    "vertex.label.cex=.35,\n",
    "vertex.size=3,\n",
    "edge.arrow.size=.25,\n",
    "layout=layout.matrix.1p,\n",
    "mark.groups=wtc.otago.p$membership,\n",
    "mark.col=\"green\")\n",
    " \n",
    "# It is clear that the increase in the connectivity of the web with parasites has led to\n",
    "# a larger densely connected community"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# It is clear that the increase in the connectivity of the web with parasites has led to\n",
    "# a larger densely connected community\n",
    " \n",
    "# The degree distribution of a food web can tell us a lot about the amount of specialization and\n",
    "# generalization in the web (in degree), as well as vulnerability (out degree)\n",
    " \n",
    "deg.otago<-degree(otago.graph)\n",
    "deg.otago.p<-degree(otago.graph.p)\n",
    " \n",
    "# Using the degree distribution gives a better way to visualize any differences\n",
    "# Looking at the in degree tells us about how general the diets of consumers are\n",
    "dd.otago.in<-degree.distribution(otago.graph,mode=\"in\",cumulative=T)\n",
    "dd.otago.in.p<-degree.distribution(otago.graph.p,mode=\"in\",cumulative=T)\n",
    " \n",
    "# Out degree is a measure of the vulnerability of organisms, telling us how many consumers\n",
    "# eat each species.\n",
    "dd.otago.out<-degree.distribution(otago.graph,mode=\"out\",cumulative=T)\n",
    "dd.otago.out.p<-degree.distribution(otago.graph.p,mode=\"out\",cumulative=T)\n",
    " \n",
    "par(mfrow=c(2,2))\n",
    "plot(dd.otago.in,xlim=c(0,80))\n",
    "plot(dd.otago.out,xlim=c(0,80))\n",
    "plot(dd.otago.in.p,xlim=c(0,80))\n",
    "plot(dd.otago.out.p,xlim=c(0,80))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# And finally the degree (\"all\") simply tells us about how well connected that species is\n",
    "# within the network\n",
    "dd.otago<-degree.distribution(otago.graph,mode=\"all\",cumulative=T)\n",
    "dd.otago.p<-degree.distribution(otago.graph.p,mode=\"all\",cumulative=T)\n",
    " \n",
    "power.fit<-power.law.fit(deg.otago)\n",
    "power.fit.p<-power.law.fit(deg.otago.p)\n",
    " \n",
    "par(mfrow=c(1,2))\n",
    "plot(dd.otago,log=\"xy\")\n",
    "lines(1:180,10*(1:180)^((-power.fit$alpha)+1))\n",
    " \n",
    "plot(dd.otago.p,log=\"xy\")\n",
    "lines(1:100,10*(1:100)^((-power.fit.p$alpha)+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# I can look at the diameter of the two versions of the web\n",
    "# For food webs the diameter is going to be the longest food chain\n",
    "# since energy only flows in one direction, the diameter will read from\n",
    "# basal species to top predator.\n",
    " \n",
    "get.diameter(otago.graph)\n",
    "get.diameter(otago.graph.p)\n",
    " \n",
    "# I think that here it is interesting to note that the diameter of the predator-prey only\n",
    "# food web (which we expect to be smaller) is not a subset of the diameter for the\n",
    "# larger parasites included network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# The next few properties are all related to the small world-ness of the network:\n",
    "transitivity(otago.graph)\n",
    "transitivity(otago.graph.p)\n",
    " \n",
    "# Betweenness is the number of shortest paths going through a specified node or edge\n",
    "otago.between<-betweenness(otago.graph)\n",
    "otago.between.p<-betweenness(otago.graph.p)\n",
    " \n",
    "plot(otago.between[1:123]~otago.between.p)\n",
    "abline(a=0,b=1)\n",
    " \n",
    "otago.edge.between<-edge.betweenness(otago.graph)\n",
    "otago.edge.between.p<-edge.betweenness(otago.graph.p)\n",
    " \n",
    "head(closeness(otago.graph))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Here are the adjacency matrices for each of the 13 subgraphs again\n",
    "s1<-matrix(c(0,1,0,0,0,1,0,0,0),nrow=3,ncol=3)\n",
    "s2<-matrix(c(0,1,1,0,0,1,0,0,0),nrow=3,ncol=3)\n",
    "s3<-matrix(c(0,1,0,0,0,1,1,0,0),nrow=3,ncol=3)\n",
    "s4<-matrix(c(0,0,1,0,0,1,0,0,0),nrow=3,ncol=3)\n",
    "s5<-matrix(c(0,1,1,0,0,0,0,0,0),nrow=3,ncol=3)\n",
    "d2<-matrix(c(0,1,1,1,0,1,0,0,0),nrow=3,ncol=3)\n",
    "d1<-matrix(c(0,1,1,0,0,1,0,1,0),nrow=3,ncol=3)\n",
    "d3<-matrix(c(0,0,1,1,0,0,1,0,0),nrow=3,ncol=3)\n",
    "d4<-matrix(c(0,0,0,1,0,1,0,1,0),nrow=3,ncol=3)\n",
    "d5<-matrix(c(0,1,1,0,0,1,1,0,0),nrow=3,ncol=3)\n",
    "d6<-matrix(c(0,1,1,1,0,1,1,1,0),nrow=3,ncol=3)\n",
    "d7<-matrix(c(0,1,1,1,0,1,1,0,0),nrow=3,ncol=3)\n",
    "d8<-matrix(c(0,1,1,1,0,0,1,0,0),nrow=3,ncol=3)\n",
    " \n",
    "# Turn them into a convenient list\n",
    "subgraph3.mat<-list(s1,s2,s3,s4,s5,d1,d2,d3,d4,d5,d6,d7,d8)\n",
    "# And then into a list of graph objects\n",
    "subgraph3.graph<-lapply(subgraph3.mat,graph.adjacency)\n",
    " \n",
    "# Count the number of the 13 different 3-node subgraphs in the two webs\n",
    "subgraph.freq.otago<-c()\n",
    "subgraph.freq.otago.p<-c()\n",
    "for(i in 1:13){\n",
    "subgraph.freq.otago[i]<-\n",
    "graph.count.subisomorphisms.vf2(otago.graph,subgraph3.graph[[i]],vertex.color1 = NULL,vertex.color2 = NULL,edge.color1 = NULL,edge.color2=NULL)\n",
    "subgraph.freq.otago.p[i]<-graph.count.subisomorphisms.vf2(otago.graph.p,subgraph3.graph[[i]],vertex.color1 = NULL,vertex.color2 = NULL,edge.color1 = NULL,edge.color2=NULL)\n",
    "}\n",
    " \n",
    "plot(subgraph.freq.otago,type=\"o\",lty=3, xlab=\"Subgraph\",ylab=\"Frequency\")\n",
    "points(subgraph.freq.otago.p,type=\"o\",lty=2)\n",
    " \n",
    "plot(subgraph.freq.otago~subgraph.freq.otago.p)\n",
    "abline(a=0,b=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}