<commit_before>{
 "metadata": {
  "name": "",
  "signature": "sha256:91aa98cc703e74fa4d29d40a63c0c9f69520984f8e247bb2b8154aa41c6501cd"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "<center>Structural Analysis and Visualization of Networks</center>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "<center>Home Assignment #2: Network models</center>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "<center>Student: *<NAME>*</center>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "<hr />\n",
      "General Information"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Due Date:** 18.02.2015 23:59 <br \\>\n",
      "**Late submission policy:** -0.2 points per day <br \\>\n",
      "\n",
      "\n",
      "Please send your reports to <mailto:<EMAIL>> and <mailto:<EMAIL>> with message subject of the following structure:<br \\> **[HSE Networks 2015] *Nazarov* *Ivan* HA*1***\n",
      "\n",
      "Support your computations with figures and comments. <br \\>\n",
      "If you are using IPython Notebook you may use this file as a starting point of your report.<br \\>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<hr \\>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Problem statement"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Task 1\n",
      "\n",
      "Consider Barabasi and Albert dynamical grow model. Two main ingredients of this model are *network growing* and *prefferential attachment*. Implement two restricted B&A-based models:  \n",
      "\n",
      " **Model A**  \n",
      "  Lack of prefferential attachment, that is at each time-step form edges uniformly at random while network keeps growing.\n",
      "\n",
      " **Model B**  \n",
      "  Lack of growing, that is fix total number of nodes, on each time-step randomly choose one and form edges with prefferential attachment."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Task 2\n",
      "Consider the following  \"Vertex copying model\" of growing network.\n",
      "\n",
      "At every time step a random vertex from already existing vertices is selected and duplicated together with all edges, such that  every edge of the  vertex\n",
      "* is copied with probability $q$\n",
      "* is rewired to any other randomly selected vertex with probability $1-q$\n",
      "\n",
      "\n",
      "Starting state is defined by some small number of randomly connected vertices.\n",
      "\n",
      "The model can generate both directed and undirected networks."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Objectives\n",
      "1. Generate networks according to the models above ($N > 1000$ nodes)\n",
      "2. Compute CDF/PDF, describe the distribution and compute\\describe its properties.\n",
      "3. Illustate the following dependencies: \n",
      "    * average path length to the number of nodes\n",
      "    * average clustering coefficient to the number of nodes\n",
      "    * average node degee to the nodes \"age\"\n",
      "4. Is scale-free property conserved in these models?\n",
      "\n",
      "Analyse results with respect to various parameter settings"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<hr/><hr/>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Preamble\n",
      "We shall include all the necessary modules here, before any analysis."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## As usual, attach the numpy and the netwrokx modules\n",
      "import networkx as nx\n",
      "import numpy as np\n",
      "import numpy.random as rnd\n",
      "import numpy.matlib as ml\n",
      "from scipy.stats import chisquare\n",
      "\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "%matplotlib inline\n",
      "\n",
      "import scipy.sparse as sp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### The original model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First and foremost, but solely for reference purposes, let's state the Barabasi-Albert dynamical growth model in its original form (**Barabasi, Albert (1999)**).\n",
      "\n",
      "The evolution of the network is represented by a sequence of graphs ${(G_n)}_{n\\geq0}$ with $G_n = \\big(V_n, E_n\\big)$.\n",
      "The process starts with a graph $G_0$ of $m_0$ nodes. The initial number of vertices is usually a small number.  At every moment $n$ the network $G_n$ undergoes evolution, which can be divided in two phases  \n",
      " * a new vertex $w'\\notin V_n$ is introduced into the network: $V_{n+1} = V_n \\cup \\{w'\\}$;\n",
      " * the vertex forms $m$ ($\\leq m_0$) links to pre-existing nodes chosen at random without replacement according to some distribution $D_n$ over the nodes $V_n$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the Barabasi-Albert growth model $D_n$ is the ``preferential connection'' distribution: for every $v\\in V_n$  \n",
      "\n",
      "$$\\mathbb{P}\\big(w'\\leadsto v\\big) = \\frac{\\delta_n(v)}{\\sum_{u\\in V_n}\\delta_n(u) }$$\n",
      "where $\\delta_n(v)$ is the degree of vertex $v$ in graph $G_n$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## The original attachment distribution of Barabasi-Albert\n",
      "def ba_valency( G, nbunch = None ) :\n",
      "## Get the distribution of the valency of the existing nodes base on the observed\n",
      "##  connectivity and degrees.\n",
      "    deg = G.degree( nbunch )\n",
      "## For some reason the call to $order( G )$ takes exceptonally long time to return\n",
      "##  But still, we need the numebr of edges here.\n",
      "    deg_sum = float( sum( deg.values( ) ) )\n",
      "## The vector of valencies is in the same order as the thr nbunch (or nx.nodes(.))\n",
      "    nbunch = nx.nodes( G ) if nbunch is None else nbunch\n",
      "    return np.fromiter( ( deg[ n ] / deg_sum for n in nbunch ),\n",
      "        count = len( nbunch ), dtype = np.float )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Analysis of Barabasi-Albert model\n",
      "#### Motivation\n",
      "\n",
      "In their exposition Barabasi-Albert never state that the initial network $G_0$ is connected, although it should be noted that, unless $m\\neq m_0$, the first ever attachment step brings about a connected star graph. Nevertheless, the model is stated as follows : (see p.551 c.1 L 30-35 of **<NAME> (1999)**) \n",
      "> To incorporate the growing character of the network, starting with a small number ($m_0$) of vertices, at every time step we add a new vertex with $m$ ($\\leq \udbff\udc14m_0$) edges that link the new vertex to $m$ different vertices already present in the system.\n",
      "\n",
      "Furthermore : (see p.551 c.1 L 40-42 of **<NAME> (1999)**)\n",
      "> After $t$ time steps, the model leads to a random network with $t\udbff\udc04 + m_0$ vertices and $mt$ edges.\n",
      "\n",
      "Note that had the initial network $G_0$ been connected, at least minimally, it would have had at least $mt + (m_0-1)$ edges by the end of step $t$. The $m_0-1$ term comes form the fact that the minimal spanning tree is the minimal connected graph."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, the case $m < m_0$ does not seem unnatural as it models the situation when the edge fromation is dominated by vertex discovery rather than preferential attachment.\n",
      "\n",
      "In light of this observation the previous definition of attachment probability dstribution suffers from the following drawback: isolated vertices remain almost surely isolated.\n",
      "\n",
      "In my opinion, this does not seem to be what authors had had in mind when they devised their model and attempted to check if both evolutionary phases were indeed necessary for the development of the power law scaling (see p.551 c.2 L 8-9 of **Barabasi, Albert (1999)** ) : \n",
      "> In model B, we start with $N$ vertices and no edges. At each time step, we randomly select a vertex and connect it with probability $\\Pi\udbff\udc0a(k_i) =\udbff\udc06 \\frac{k_i}{\udbff\udc0b\\sum_j k_j}$ to vertex $i$ in the system.\n",
      "\n",
      "This definition of attachment probabilities would have no effect on the Barabasi-Albert model and its type-B restriction even in the case when $m < m_0$, if the evolution began with a connected graph $G_0$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "References:\n",
      "* Barab\u00e1si, <NAME>., & <NAME>. (1999). \"_Emergence of scaling in random networks._\" Science, 286(5439), 509-512.\n",
      "Chicago\t"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Bayesian preferential attachment\n",
      "Nevertheless it is easy to patch this \"vagueness\" in the model.\n",
      "Let's imagine for a moment that for some observed graph $G = (V,E)$ we are faced with the problem of estimating the probabilities ${(\\theta_v)}_{v\\in W}$ of getting attached to a node in some $W\\subseteq V$. Then the vertex degrees of the nodes in $W$ can be regarded as the evidence of attachement, whence the likelihood function should be given by:  \n",
      "\n",
      "$$L\\Big(G; {(\\theta_v)}_{v\\in W}\\Big) = \\prod_{v\\in W} \\theta_v^{\\delta(v)}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Brief theory\n",
      "For this task of parameter estimation, let's employ Bayesian approach instead of the common MLE, since the latter tends to overfit the model because it fails to incorporate the uncertainty of the parameters.\n",
      "\n",
      "Cosider a prior $\\pi(\\cdot)$ on the space of possible probability vectors ${(\\theta_v)}_{v\\in W}$ -- a $(n-1)$-dimensional simplex $S^\\circ_{n-1}\\subseteq {[0,1]}^n$ where $n=\\big\\rvert W\\big\\lvert$. Using Bayes' theorem, the posterior distribution of ${(\\theta_v)}_{v\\in W}$ given the observed data is  \n",
      "\n",
      "$$p\\Big(\\Big.{(\\theta_v)}_{v\\in W}\\,\\big.\\big\\rvert\\,G\\Big.\\Big)\n",
      "\\propto p\\Big(\\Big.G\\,\\big.\\big\\rvert\\,{(\\theta_v)}_{v\\in W}\\Big.\\Big) \\pi\\Big({(\\theta_v)}_{v\\in W}\\Big)\n",
      "= L\\Big(G; {(\\theta_v)}_{v\\in W}\\Big) \\pi\\Big({(\\theta_v)}_{v\\in W}\\Big)$$\n",
      "\n",
      "In general it is convenient to pick a prior $\\pi(\\cdot)$ from an \"**eigen-family**\" of this prior-to-posterior tranfromation. Such prior families are know nas the **conjugate prior families**. For this particular kind of likelihood such a family is the **Dirichlet distribution** given by  \n",
      "\n",
      "$$\\text{Dir}\\big({(\\alpha_v)}_{v\\in W}\\big) = \\frac{\\Gamma\\big( \\sum_{v\\in W} \\alpha_v \\big)}{ \\prod_{v\\in W} \\Gamma\\big(\\alpha_v \\big)} \\prod_{v\\in W} \\theta_v^{\\alpha_v-1}$$\n",
      "\n",
      "It is easy to demonstrate that the posterior is $\\text{Dir}\\big({(\\alpha_v)}_{v\\in W} + {(\\delta(v))}_{v\\in W}\\big)$. This implies that the posterior probability of a random vertex being connected to a node $v\\in W$ is given by  \n",
      "$$p\\Big(\\Big. v\\,\\big.\\big\\rvert\\,G\\Big.\\Big)\n",
      "= \\int p\\Big(\\Big. v\\,\\big.\\big\\rvert\\,{(\\theta_v)}_{v\\in W}\\Big.\\Big) p\\Big(\\Big. {(\\theta_v)}_{v\\in W}\\, \\big.\\big\\rvert\\,G\\Big.\\Big) d\\mathbf{\\theta} = \\ldots $$  \n",
      "\n",
      "which in our case resolves into  \n",
      "$$\\ldots = \\int \\theta_v \\text{Dir}\\big({(\\alpha_v)}_{v\\in W} + {(\\delta(v))}_{v\\in W}\\big) d\\mathbf{\\theta}\n",
      "= \\frac{\\alpha_v + \\delta(v)}{\\sum_{u\\in W} \\alpha_u + \\delta(u)}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Final probability vector\n",
      "Since there is no reason to assume any specific prior distribution of the values of the attachment probability vectors, let's use an uninformative prior, which assigna equal probability of every vector in the simplex. In the case of a Dirichlet family, this sets the hyperparameter ${(\\alpha_v)}_{v\\in W}$ of the prior to a vector of ones: $\\alpha_v=1$ for each $v\\in W$.\n",
      "\n",
      "Therefore the attachment probability distribution in $G_n$ is given by  \n",
      "\n",
      "$$p\\Big(\\Big. w'\\leadsto v\\,\\big.\\big\\rvert\\,G_n\\Big.\\Big) = \\frac{\\delta_n(v)+1}{\\big\\lvert V_n\\big\\rvert + \\sum_{u\\in S}\\delta_n(u)}$$  \n",
      "\n",
      "Listed below is the procedure, which computes the required distribution for a given graph $G$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## Bayesian valency distribution. Funny thing is that this correction\n",
      "##  turns out to is\n",
      "def ba_bayes_valency( G, nbunch = None ) :\n",
      "## Get the degree of nodes\n",
      "    deg = G.degree( nbunch )\n",
      "## Get the adjusted sum of vertex degrees\n",
      "    nbunch = nx.nodes( G ) if nbunch is None else nbunch\n",
      "    deg_sum = float( sum( deg.values( ) ) + len( nbunch ) )\n",
      "## Compute the posterior distribution of vertex valency\n",
      "    return np.fromiter( ( ( deg[ n ] + 1.0 ) / deg_sum for n in nbunch ),\n",
      "        count = len( nbunch ), dtype = np.float )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bayesian smoothing reduces the zero-degree penalty on the isolated vertices, thus enabling vertex discovery during the early stages of evolution. Notable thing is that for an uninformative prior the posterior distribution given $G$ is equivalent to the original attachment distribution for the same graph with an extra loop at every vertex."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### The long run distribution\n",
      "Using the heuristic \"mean field\" approach it is possible to derive an approximation to the digree distribution un the long run.\n",
      "\n",
      "Suppose vertices are born at intervals $\\frac{1}{n}$ with $m$ edges. Then at time $t\\in [0,1]$ there are $nt$ vertices and the sum of all node degrees is $2mnt$. Across the span of an infinitesimal interval $[t, t+dt]$ extra $n dt$ vertices are born, each with $m$ edges.\n",
      "\n",
      "Suppose a vertex $v$ is born at time $s$ with $m$ edges and at time $t>s$ this node $v$ has degree $k_s(t)$.  Then during $[t, t+dt]$ the vertex $v$ gets an incoming edge with probability $n dt \\frac{k_s(t)+1}{2mnt + nt} = \\frac{k_s(t)+1}{(2m+1)t} dt$. Ov average thus $v$ gets $\\frac{m}{2m+1} \\frac{k_s(t)+1}{t}dt$ edges, which is by how much its degree increses :  \n",
      "\n",
      "$$ \\frac{k_s(t+dt)-k_s(t)}{dt} = \\frac{m}{2m+1} \\frac{k_s(t)+1}{t}$$\n",
      "\n",
      "Note that the degree depends on the time thevertex was born and not on the vertex itself."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the continuous time limit as $n\\to \\infty$ the following differential equation approximates the dynamics of degree of the vertex $v$ :  \n",
      "\n",
      "$$\\dot{k}_s(t) = \\frac{m}{2m+1} \\frac{k_s(t)+1}{t}$$\n",
      "\n",
      "with initial value at time of birth being $m$ : $k_s(t=s) = m$.\n",
      "\n",
      "The solution to this ordinary differential equation for these initial conditions is given by  \n",
      "$$\\delta_s(t) = (m+1) \\Big(\\frac{t}{s}\\Big)^\\frac{m}{2m+1}$$\n",
      "Put $\\alpha = \\frac{m}{2m+1}$ and note that $\\frac{1}{\\alpha} = 2 + \\frac{1}{m}$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the continuous time limit the probability that at time $t$ a vertex born at $s$ has degree less than $k$ is given by \n",
      "$$F(k) = \\mathbb{P}(\\delta_s(t)\\leq k) = \\mathbb{P}\\Bigg(t \\Big(\\frac{m+1}{k+1}\\Big)^\\frac{1}{\\alpha} \\leq s\\Bigg) = \\frac{t - t \\Big(\\frac{m+1}{k+1}\\Big)^\\frac{1}{\\alpha}}{t} = 1 - \\Big(\\frac{m+1}{k+1}\\Big)^\\frac{1}{\\alpha}$$\n",
      "\n",
      "The distribution function is given by \n",
      "$$\\frac{d}{d k} F(k) = \\frac{1}{\\alpha} \\Big( k+1 \\Big)^{-\\frac{1}{\\alpha}-1} \\Big( m+1 \\Big)^\\frac{1}{\\alpha} = \\Big( k+1 \\Big)^{-3-\\frac{1}{m}} \\Big(2+\\frac{1}{m}\\Big) \\Big( m+1 \\Big)^{2+\\frac{1}{m}}$$\n",
      "\n",
      "Thus means that the limiting degree distribution of the Barabasi-Albert model with Bayesian valency behaves approximately like the power law with similar exponent as in the original model.\n",
      "\n",
      "For a more formally rigorous treatment of the original model see _**Bollobas, Riordan et al., (2001)**_."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "References:\n",
      "* Bollob\u00e1<NAME>., <NAME>., <NAME>., & Tusn\u00e1<NAME>. (2001). \"_The degree sequence of a scale\u2010free random graph process._\" Random Structures & Algorithms, 18(3), 279-290."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Model evolution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### _Complete model_\n",
      "A has been mentioned above the basic growth step in the complete Barabasi-Albert consisted of two phases:\n",
      " 1. growth;\n",
      " 2. attachment according to valency distribution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def barabasi_albert( n, m = 10, callback = None, avclu = False, avlen = False ) :\n",
      "## Make sure the callback is a proper function\n",
      "\tcallback_fn = callback if callable( callback ) else lambda G, _, __: None\n",
      "## Barabasi-Albert evolution always begins with the isolated graph\n",
      "\tG = nx.empty_graph( m )\n",
      "## To speed up the computations keep track of the node degree locally,\n",
      "##  but instead of using a dictionary, store them in a linear array.\n",
      "\tdegree = np.zeros( n + m, dtype = np.int )\n",
      "\tdegree_half_sum = 0\n",
      "## Local adjacency matrix : should be made sparse\n",
      "\tadjacency = np.zeros( ( n + m, n + m ) )\n",
      "## Shortest path matrix for the modified Floyd Warshall\n",
      "\tpi = 1 - ml.identity( n + m, np.float )\n",
      "\tpi[ pi == 1 ] = np.inf\n",
      "\tnpmin = np.minimum\n",
      "## Local tracking of the number of triangles for the clustering coefficient\n",
      "\ttriangles = np.zeros( n + m, dtype = np.int )\n",
      "## Note that the graph itself is used just to keep track of the edges.\n",
      "\tresult = list( )\n",
      "## Create an array of average path lengths\n",
      "\tavg_path_ln = np.full( n + m, np.inf, dtype = np.float )\n",
      "## Setup the average clustering metric storage\n",
      "\tavg_tri_clu = np.full( n + m, np.inf, dtype = np.float )\n",
      "## Age-degree: the dependence of the average node degree on the age of a vertex\n",
      "\tavg_age_deg = np.zeros( n + 1, dtype = np.float )\n",
      "## Identify vertices by consequitive integers. \n",
      "\tt = G.number_of_nodes( )\n",
      "## In the original Barabasi-Albert model the first new vertex\n",
      "##  is attached to all m vertices, thus making the graph connected.\n",
      "\tU = range( m )\n",
      "\twhile t < n + m :\n",
      "## Add pending edges\n",
      "\t\tG.add_edges_from( ( t, u ) for u in U )    \n",
      "## maintain a local copy of the adjacency matrix\n",
      "\t\tadjacency[ t, U ] = adjacency[ U, t ] = 1\n",
      "## Record the vertex degree\n",
      "\t\tdegree[ U ] += 1\n",
      "\t\tdegree[ t ] += len( U )\n",
      "\t\tdegree_half_sum += len( U )\n",
      "## Create a new vertex using stack allocation -- on top of others.\n",
      "\t\tt = G.number_of_nodes( )\n",
      "## Run the Floyd-Warshall algorithm\n",
      "\t\tif avlen :\n",
      "## Add the shortest paths from the affected vertices to the new node\n",
      "\t\t\tpi[ U, t-1 ] = pi[ t-1, U ] = 1\n",
      "## Compute the shortest path from the exisitng nodes to the new one\n",
      "\t\t\tfor u in U :\n",
      "\t\t\t\tpi[ :t, t-1 ] = npmin( pi[ :t, t-1 ], pi[ :t, u ] + pi[ u, t-1 ] )\n",
      "## Restore the distance symmetry due to bidirectionality\n",
      "\t\t\tpi[ t-1, :t ] = pi[ :t, t-1 ].T\n",
      "## Compare the shortest paths through the new vertex with the current paths\n",
      "##\tAt this point the matrix pi has the following property:\n",
      "##\t * for nodes $i,j\\in V$ the value $\\pi_{ij}$ is the shortest path\n",
      "##\t   through nodes other than $\\omega$;\n",
      "##\t * for any $x\\in V$ the element $\\pi_{x\\omega}$ (and $\\pi_{\\omega x}) is\n",
      "##\t   the length of the shortest path from $i$ to the new node $\\omega$ with\n",
      "##\t   intermediate vertices exclusively in $V$.\n",
      "\t\t\tpi[ :t, :t ] = npmin( pi[ :t, :t ], pi[ :t, t-1 ] + pi[ t-1, :t ] )\n",
      "## Numerically confirm the correctness shown elsewhere in this notebook\n",
      "##\t\t\tassert( np.all( nx.floyd_warshall_numpy( G ) == pi[ :t, :t ] ) )\n",
      "## Compute the average path length\n",
      "##             avg_path_ln[ t-1 ] = nx.average_shortest_path_length( G )\n",
      "\t\t\tavg_path_ln[ t-1 ] = np.mean( pi[ :t, :t ], dtype = np.float64 ) * t / ( t - 1.0 )\n",
      "## Compute the average clustering coefficient\n",
      "\t\tif avclu :\n",
      "\t\t\ttri_num = 0\n",
      "\t\t\tfor u in U :\n",
      "## Count the number of edges from u to U, since they are responsible for the new triangles containing t\n",
      "\t\t\t\tX = np.sum( adjacency[ u, U ] )\n",
      "## Update the number of the current endpoint\n",
      "\t\t\t\ttriangles[ u ] += X\n",
      "\t\t\t\ttri_num += X\n",
      "## Every new triangle could have been created only via the edges incident on t\n",
      "##  The loop above counts the number of trinagles containing t twice (once per each edge of a triangle)\n",
      "\t\t\ttriangles[ t-1 ] += tri_num // 2\n",
      "## Compute the clustering coefficinet of each node\n",
      "\t\t\tmax_clust = degree[ :t ] * ( degree[ :t ] - 1 )\n",
      "\t\t\tclust = np.where( max_clust != 0, ( 2.0 * triangles[ :t ] ) / max_clust, 0 )\n",
      "## And the average\n",
      "\t\t\tavg_tri_clu[ t-1 ] = np.mean( clust, dtype = np.float64 )\n",
      "# \t\t\tcc = nx.average_clustering( G )\n",
      "# \t\t\tavg_tri_clu[ t-1 ] = np.abs( ( np.mean( clust, dtype = np.float64 ) - cc ) / ( cc + 1e-12 ) )\n",
      "## Compute the average node degree of created vertices of a given age.\n",
      "## * A[0:t-m0] -- current mean degrees of vertices of age 0:t-m0 (right ends not included)\n",
      "## * D[0:m0] -- the current degrees of the oldest (initial) m0 vertices\n",
      "## * D[m0:t] -- the degrees of nodes created at iterations 0:t-m0 respectively\n",
      "\t\tavg_age_deg[ :t-m ] += ( degree[ t-1 :m-1 :-1 ] - avg_age_deg[ :t-m ]\n",
      "\t\t\t) / np.arange( t-m, 0, -1, dtype = np.float )\n",
      "## Now incorporate the degree of the initial $m$ nodes, which are the oldest ones.\n",
      "\t\tavg_age_deg[ t-m ] = np.mean( degree[ :m ], dtype = np.float64 )\n",
      "## The evolution step is complete, call the required function\n",
      "\t\tresult.append( callback_fn( G, degree[ :t ], pi[ :t, :t ] ) )\n",
      "## The initial nodes contribute to the average node degree of the next generation of vertices\n",
      "## The nodes which the new vertex is attached to are picked\n",
      "##  according to the preferential attachment distribution.\n",
      "\t\tprobabilities = degree[ :t ] / ( 2.0 * degree_half_sum )\n",
      "\t\tU = rnd.choice( t, size = m, p = probabilities, replace = False )\n",
      "#     assert( np.all( nx.floyd_warshall_numpy( G ) == pi ) )\n",
      "\treturn G, result, avg_age_deg, avg_tri_clu, avg_path_ln\n",
      "## https://networkx.github.io/documentation/latest/reference/classes.graph.html\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## %lprun -f barabasi_albert barabasi_albert( 1000, m = 2, avlen = True, avclu = True )\n",
      "# G1, R1, A1, C1, L1 = barabasi_albert( 1000, m = 10, avlen = True, avclu = True )\n",
      "# plt.loglog( C1 )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### _Model A_\n",
      "In **model A**, the first restricted verison of the evolution, the network grows with time but the newly added nodes are attached to the existing ones uniformly at random.\n",
      "The code below implements the basic steps of both the Bayesian and the first restricted model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## This evolution step is similar to the original BA step, except it\n",
      "##  relaxes the requirement that m_0 be equal to m, and it implemets\n",
      "##  the egalitarian attachment and bayesian valency.\n",
      "def barabasi_albert_model_a( n, m0 = 10, m = 5, bayes = False, callback = None, avclu = False, avlen = False ) :\n",
      "\tcallback_fn = callback if callable( callback ) else lambda G, _, __: None\n",
      "\tresult = list( )\n",
      "\tG = nx.empty_graph( m0 )\n",
      "## Local degree and shortest path tracking\n",
      "\tdegree = np.zeros( n + m0, dtype = np.int )\n",
      "\tdegree_half_sum = 0\n",
      "\tadjacency = np.zeros( ( n + m0, n + m0 ) )\n",
      "\tpi = 1 - ml.identity( n + m0, np.float )\n",
      "\tpi[ pi == 1 ] = np.inf\n",
      "\tnpmin = np.minimum\n",
      "\ttriangles = np.zeros( n + m0, dtype = np.int )\n",
      "## Average path length and clustering\n",
      "\tavg_path_ln = np.full( n + m0, np.inf, dtype = np.float )\n",
      "\tavg_tri_clu = np.full( n + m0, np.inf, dtype = np.float )\n",
      "\tavg_age_deg = np.zeros( n + 1, dtype = np.float )\n",
      "## Begin the simulation\n",
      "\tt = G.number_of_nodes( )\n",
      "\twhile t < n + m0 :\n",
      "## Select either the bayesian valency distribution, or uniform attachment\n",
      "\t\ttheta = None if not bayes else ( degree[ :t ] + 1.0 ) / ( 2.0 * degree_half_sum + t )\n",
      "\t\tU = rnd.choice( t, size = m, p = theta, replace = False )\n",
      "\t\tG.add_edges_from( ( t, u ) for u in U )\n",
      "\t\tadjacency[ t, U ] = adjacency[ U, t ] = 1\n",
      "## Accumulate degree information\n",
      "\t\tdegree[ U ] += 1\n",
      "\t\tdegree[ t ] += len( U )\n",
      "\t\tdegree_half_sum += len( U )\n",
      "## Preallocate the next vertex\n",
      "\t\tt = G.number_of_nodes( )\n",
      "## Compute the shortest paths\n",
      "\t\tif avlen :\n",
      "\t\t\tpi[ U, t-1 ] = pi[ t-1, U ] = 1\n",
      "\t\t\tfor u in U :\n",
      "\t\t\t\tpi[ :t, t-1 ] = npmin( pi[ :t, t-1 ], pi[ :t, u ] + pi[ u, t-1 ] )\n",
      "\t\t\tpi[ t-1, :t ] = pi[ :t, t-1 ].T\n",
      "\t\t\tpi[ :t, :t ] = npmin( pi[ :t, :t ], pi[ :t, t-1 ] + pi[ t-1, :t ] )\n",
      "\t\t\tavg_path_ln[ t-1 ] = np.mean( pi[ :t, :t ], dtype = np.float64 ) * t / ( t - 1.0 )\n",
      "## Get the average clustering coefficient\n",
      "\t\tif avclu :\n",
      "\t\t\ttri_num = 0\n",
      "\t\t\tfor u in U :\n",
      "\t\t\t\tX = np.sum( adjacency[ u, U ] )\n",
      "\t\t\t\ttriangles[ u ] += X\n",
      "\t\t\t\ttri_num += X\n",
      "\t\t\ttriangles[ t-1 ] += tri_num // 2\n",
      "\t\t\tmax_clust = degree[ :t ] * ( degree[ :t ] - 1 )\n",
      "\t\t\tclust = np.where( max_clust != 0, ( 2.0 * triangles[ :t ] ) / max_clust, 0 )\n",
      "\t\t\tavg_tri_clu[ t-1 ] = np.mean( clust, dtype = np.float64 )\n",
      "## Finalize this step's statistics\n",
      "\t\tavg_age_deg[ :t-m0 ] += ( degree[ t-1 :m0-1 :-1 ] - avg_age_deg[ :t-m0 ]\n",
      "\t\t\t) / np.arange( t-m0, 0, -1, dtype = np.float )\n",
      "\t\tavg_age_deg[ t-m0 ] = np.mean( degree[ :m0 ], dtype = np.float64 )\n",
      "\t\tresult.append( callback_fn( G, degree[ :t ], pi[ :t, :t ] ) )\n",
      "\treturn G, result, avg_age_deg, avg_tri_clu, avg_path_ln\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## %lprun -f barabasi_albert_model_a barabasi_albert_model_a( 300, m = 1, avlen = True, avclu = False )\n",
      "# G1, R1, A1, C1, L1 = barabasi_albert_model_a( 1000, m0 = 10, m=10, avlen = True, avclu = True )\n",
      "# plt.loglog( C1 )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### _Model B_\n",
      "**Model B** retains the preferential attachment phase, but foregoes the growth step. Instead, a random existing node is picked and is linked with other nodes, preferentially chosing vertices with high degree. Initally this model begins its evolution from a set of isolated nodes.\n",
      "\n",
      "For this model to work as intended the only option is to use the Bayesian valency distribution, since the originally proposed probabilities handle poorly the isolated nodes, as it has been alread stated. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###### A couple of nuances\n",
      "Let $w'\\in V$ be a randomly picked vertex. Since multiple connections and loops are forbidden in the model, the preferential attachment of a vertex $w'$ must use conditional node valency distribution. If a subset $S\\subseteq V$, $w'\\in S$, is **prohibited**, then for any vertex $v\\in V\\setminus S$  \n",
      "\n",
      "$$p\\Big(\\Big. w'\\leadsto v\\,\\big.\\big\\rvert\\, w'\\not\\!\\leadsto S,\\,G \\Big.\\Big)\n",
      "= \\frac{p\\Big(\\Big. w'\\leadsto v\\,\\big.\\big\\rvert G \\Big.\\Big)}{p\\Big(\\Big. w'\\not\\!\\leadsto S\\,\\big.\\big\\rvert G \\Big.\\Big)}\n",
      "= \\frac{\\delta_n(v)+1}{\\big\\lvert V\\setminus S\\big\\rvert + \\sum_{u\\in V\\setminus S}\\delta_n(u)} $$  \n",
      "\n",
      "In particular, the prohibited set for a vertex $w'$ is itself and the set of its immediate neighbours:  \n",
      "\n",
      "$$S = \\big\\{w'\\big\\}\\cup N_G(w')$$\n",
      "\n",
      "Another subtle point is the connection strategy when $\\big\\lvert N_G(w') \\big\\rvert\\geq n-1-m$. If $\\big\\lvert N_G(w') \\big\\rvert < n-1$, the model forcefully connects $w'$ to all of the remaining vertices. If however $w'$ is saturated, then the attachment step for this node is skipped as this vertex is just unable to form new connections."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## Another relaxation of the complete Barabasi-Albert model, which eliminates\n",
      "##  the growth phase. This restriction uses the Bayesian valency distribution.\n",
      "##  when performing preferential attachment.\n",
      "def barabasi_albert_model_b( n, m = 5, callback = None, avclu = False, avlen = False ) :\n",
      "\tcallback_fn = callback if callable( callback ) else lambda G, _, __: None\n",
      "\tG = nx.empty_graph( n )\n",
      "## Local node degree cache\n",
      "\tdegree = np.zeros( n, dtype = np.int64 )\n",
      "\tdegree_sum = 0\n",
      "\tadjacency = np.zeros( ( n, n ), dtype = np.int )\n",
      "## Shortest path matrix for the modified Floyd Warshall\n",
      "\tpi = 1 - ml.identity( n, np.float )\n",
      "\tpi[ pi == 1 ] = np.inf\n",
      "\tnpmin = np.minimum\n",
      "\ttriangles = np.zeros( n, dtype = np.int )\n",
      "## Average path length and clustering\n",
      "\tavg_tri_clu = np.full( n, np.inf, dtype = np.float )\n",
      "\tavg_path_ln = np.full( n, np.inf, dtype = np.float )\n",
      "\tavg_age_deg = np.zeros( n, dtype = np.float )\n",
      "## Begin the simulation\n",
      "\tt = 0\n",
      "\tresult = list( )\n",
      "\twhile t < n :\n",
      "## \"Pick an existing one at random\" (cf. p.551 c.3 L 9 of [Barabasi, Albert; 1999])\n",
      "##  ... and get its neighbours.\n",
      "\t\tv = rnd.randint( n, size = 1 )[ 0 ]\n",
      "\t\tS = nx.neighbors( G, v )\n",
      "## If it has more neighbours than n-1-m, connect it to the remaining nodes\n",
      "##  and continue. If the node is staturated, ignore it.\n",
      "\t\tm_prime = min( m, n - 1 - len( S ) )\n",
      "\t\tif m_prime > 0 :\n",
      "## Get the correction terms due to conditioning (see the formula).\n",
      "\t\t\tc_term = degree[ v ] + 1.0 + np.sum( degree[ S ] ) + len( S )\n",
      "## Compute the conditional valency of each vertex in the current configuration\n",
      "\t\t\ttheta = ( degree + 1.0 ) / ( 0.0 + degree_sum + n - c_term )\n",
      "## Zero out the prohibited vertices, so that they'll never be picked.\n",
      "\t\t\ttheta[ v ] = 0 ; theta[ S ] = 0\n",
      "## Pick $m'$ other nodes to preferentially link the chosen one to.\n",
      "\t\t\tU = rnd.choice( n, size = m_prime, p = theta, replace = False )\n",
      "\t\t\tG.add_edges_from( ( v, u ) for u in U )\n",
      "\t\t\tadjacency[ v, U ] = adjacency[ U, v ] = 1\n",
      "\t\t\tdegree[ U ] += 1\n",
      "\t\t\tdegree[ v ] += m_prime\n",
      "\t\t\tdegree_sum += m_prime * 2 ;\n",
      "## Compute the shortest paths\n",
      "\t\tif avlen :\n",
      "\t\t\tif m_prime > 0 :\n",
      "\t\t\t\tif False :\n",
      "## Incrementally update the shortest path matrix by every new bidirectional edge\n",
      "\t\t\t\t\tfor u in U :\n",
      "## Compute the path through this edge\n",
      "\t\t\t\t\t\tpath = pi[ :, u ] + 1.0 + pi[ v, : ]\n",
      "## Choose the optimal edge direction and whether the new path is the shortest\n",
      "\t\t\t\t\t\tpi = npmin( pi, npmin( path, path.T ) )\n",
      "\t\t\t\telse :\n",
      "\t\t\t\t\tpi[ U, v ] = pi[ v, U ] = 1\n",
      "\t\t\t\t\tif False :\n",
      "\t\t\t\t\t\tfor u in U :\n",
      "\t\t\t\t\t\t\tpi = npmin( pi, pi[ :, u ] + pi[ u, : ] )\n",
      "\t\t\t\t\telse :\n",
      "\t\t\t\t\t\tfor u in U :\n",
      "\t\t\t\t\t\t\tpi[ :, v ] = npmin( pi[ :, v ], pi[ :, u ] + pi[ u, v ] )\n",
      "\t\t\t\t\t\tpi[ v, : ] = pi[ :, v ].T\n",
      "\t\t\t\t\tpi = npmin( pi, pi[ :, v ] + pi[ v, : ] )\n",
      "# \t\t\t\tassert( np.all( nx.floyd_warshall_numpy( G ) == pi ) )\n",
      "\t\t\tavg_path_ln[ t ] = np.mean( pi ) * n / ( n - 1.0 ) if m_prime > 0 else avg_path_ln[ t-1 ]\n",
      "## Get the average clustering coefficient\n",
      "\t\tif avclu :\n",
      "\t\t\ttri_num = 0\n",
      "\t\t\tfor u in G[ v ] :\n",
      "\t\t\t\tX = np.sum( adjacency[ u, U ] )\n",
      "\t\t\t\ttriangles[ u ] += X\n",
      "\t\t\t\tif u in U :\n",
      "\t\t\t\t\ttri_num += X\n",
      "\t\t\t\telse :\n",
      "\t\t\t\t\ttriangles[ v ] += X\n",
      "## it is necessary to inform the vertices adjacent to a preexisting neighbour\n",
      "## that the edge (v,u) induced a new triangle through a newly added edge.\n",
      "\t\t\t\t\ttriangles[ U ] += adjacency[ U, u ]\n",
      "\t\t\ttriangles[ v ] += tri_num // 2\n",
      "\t\t\tmax_clust = degree * ( degree - 1 )\n",
      "\t\t\tclust = np.where( max_clust != 0, ( 2.0 * triangles ) / max_clust, 0 )\n",
      "\t\t\tavg_tri_clu[ t ] = np.mean( clust, dtype = np.float64 ) if m_prime > 0 else avg_tri_clu[ t-1 ]\n",
      "## In this model all nodes are the same age, which means that the degree-age\n",
      "##  relationship is just the average degree with respect to the generation.\n",
      "\t\tavg_age_deg[ t ] = np.mean( degree ) if m_prime > 0 else avg_age_deg[ t-1 ]\n",
      "\t\tresult.append( callback_fn( G, degree, pi ) )\n",
      "\t\tt += 1\n",
      "\treturn G, result, avg_age_deg, avg_tri_clu, avg_path_ln\n"
     ],
     "language": "python",<commit_msg>Justified a lean approach to shortest path computation upon edge addition<commit_after>{
 "metadata": {
  "name": "",
  "signature": "sha256:9663df4d4781a990931333362dde3f5e9d295343b2c431827629a0c2db07fd1c"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "<center>Structural Analysis and Visualization of Networks</center>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "<center>Home Assignment #2: Network models</center>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "<center>Student: *<NAME>*</center>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "<hr />\n",
      "General Information"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Due Date:** 18.02.2015 23:59 <br \\>\n",
      "**Late submission policy:** -0.2 points per day <br \\>\n",
      "\n",
      "\n",
      "Please send your reports to <mailto:<EMAIL>> and <mailto:<EMAIL>> with message subject of the following structure:<br \\> **[HSE Networks 2015] *Nazarov* *Ivan* HA*1***\n",
      "\n",
      "Support your computations with figures and comments. <br \\>\n",
      "If you are using IPython Notebook you may use this file as a starting point of your report.<br \\>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<hr \\>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Problem statement"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Task 1\n",
      "\n",
      "Consider Barabasi and Albert dynamical grow model. Two main ingredients of this model are *network growing* and *prefferential attachment*. Implement two restricted B&A-based models:  \n",
      "\n",
      " **Model A**  \n",
      "  Lack of prefferential attachment, that is at each time-step form edges uniformly at random while network keeps growing.\n",
      "\n",
      " **Model B**  \n",
      "  Lack of growing, that is fix total number of nodes, on each time-step randomly choose one and form edges with prefferential attachment."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Task 2\n",
      "Consider the following  \"Vertex copying model\" of growing network.\n",
      "\n",
      "At every time step a random vertex from already existing vertices is selected and duplicated together with all edges, such that  every edge of the  vertex\n",
      "* is copied with probability $q$\n",
      "* is rewired to any other randomly selected vertex with probability $1-q$\n",
      "\n",
      "\n",
      "Starting state is defined by some small number of randomly connected vertices.\n",
      "\n",
      "The model can generate both directed and undirected networks."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Objectives\n",
      "1. Generate networks according to the models above ($N > 1000$ nodes)\n",
      "2. Compute CDF/PDF, describe the distribution and compute\\describe its properties.\n",
      "3. Illustate the following dependencies: \n",
      "    * average path length to the number of nodes\n",
      "    * average clustering coefficient to the number of nodes\n",
      "    * average node degee to the nodes \"age\"\n",
      "4. Is scale-free property conserved in these models?\n",
      "\n",
      "Analyse results with respect to various parameter settings"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<hr/><hr/>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Preamble\n",
      "We shall include all the necessary modules here, before any analysis."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## As usual, attach the numpy and the netwrokx modules\n",
      "import networkx as nx\n",
      "import numpy as np\n",
      "import numpy.random as rnd\n",
      "import numpy.matlib as ml\n",
      "from scipy.stats import chisquare\n",
      "\n",
      "import matplotlib.pyplot as plt\n",
      "\n",
      "%matplotlib inline\n",
      "\n",
      "import scipy.sparse as sp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### The original model"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First and foremost, but solely for reference purposes, let's state the Barabasi-Albert dynamical growth model in its original form (**Barabasi, Albert (1999)**).\n",
      "\n",
      "The evolution of the network is represented by a sequence of graphs ${(G_n)}_{n\\geq0}$ with $G_n = \\big(V_n, E_n\\big)$.\n",
      "The process starts with a graph $G_0$ of $m_0$ nodes. The initial number of vertices is usually a small number.  At every moment $n$ the network $G_n$ undergoes evolution, which can be divided in two phases  \n",
      " * a new vertex $w'\\notin V_n$ is introduced into the network: $V_{n+1} = V_n \\cup \\{w'\\}$;\n",
      " * the vertex forms $m$ ($\\leq m_0$) links to pre-existing nodes chosen at random without replacement according to some distribution $D_n$ over the nodes $V_n$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the Barabasi-Albert growth model $D_n$ is the ``preferential connection'' distribution: for every $v\\in V_n$  \n",
      "\n",
      "$$\\mathbb{P}\\big(w'\\leadsto v\\big) = \\frac{\\delta_n(v)}{\\sum_{u\\in V_n}\\delta_n(u) }$$\n",
      "where $\\delta_n(v)$ is the degree of vertex $v$ in graph $G_n$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## The original attachment distribution of Barabasi-Albert\n",
      "def ba_valency( G, nbunch = None ) :\n",
      "## Get the distribution of the valency of the existing nodes base on the observed\n",
      "##  connectivity and degrees.\n",
      "    deg = G.degree( nbunch )\n",
      "## For some reason the call to $order( G )$ takes exceptonally long time to return\n",
      "##  But still, we need the numebr of edges here.\n",
      "    deg_sum = float( sum( deg.values( ) ) )\n",
      "## The vector of valencies is in the same order as the thr nbunch (or nx.nodes(.))\n",
      "    nbunch = nx.nodes( G ) if nbunch is None else nbunch\n",
      "    return np.fromiter( ( deg[ n ] / deg_sum for n in nbunch ),\n",
      "        count = len( nbunch ), dtype = np.float )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Analysis of Barabasi-Albert model\n",
      "#### Motivation\n",
      "\n",
      "In their exposition Barabasi-Albert never state that the initial network $G_0$ is connected, although it should be noted that, unless $m\\neq m_0$, the first ever attachment step brings about a connected star graph. Nevertheless, the model is stated as follows : (see p.551 c.1 L 30-35 of **<NAME> (1999)**) \n",
      "> To incorporate the growing character of the network, starting with a small number ($m_0$) of vertices, at every time step we add a new vertex with $m$ ($\\leq \udbff\udc14m_0$) edges that link the new vertex to $m$ different vertices already present in the system.\n",
      "\n",
      "Furthermore : (see p.551 c.1 L 40-42 of **Barabasi, Albert (1999)**)\n",
      "> After $t$ time steps, the model leads to a random network with $t\udbff\udc04 + m_0$ vertices and $mt$ edges.\n",
      "\n",
      "Note that had the initial network $G_0$ been connected, at least minimally, it would have had at least $mt + (m_0-1)$ edges by the end of step $t$. The $m_0-1$ term comes form the fact that the minimal spanning tree is the minimal connected graph."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, the case $m < m_0$ does not seem unnatural as it models the situation when the edge fromation is dominated by vertex discovery rather than preferential attachment.\n",
      "\n",
      "In light of this observation the previous definition of attachment probability dstribution suffers from the following drawback: isolated vertices remain almost surely isolated.\n",
      "\n",
      "In my opinion, this does not seem to be what authors had had in mind when they devised their model and attempted to check if both evolutionary phases were indeed necessary for the development of the power law scaling (see p.551 c.2 L 8-9 of **Barabasi, Albert (1999)** ) : \n",
      "> In model B, we start with $N$ vertices and no edges. At each time step, we randomly select a vertex and connect it with probability $\\Pi\udbff\udc0a(k_i) =\udbff\udc06 \\frac{k_i}{\udbff\udc0b\\sum_j k_j}$ to vertex $i$ in the system.\n",
      "\n",
      "This definition of attachment probabilities would have no effect on the Barabasi-Albert model and its type-B restriction even in the case when $m < m_0$, if the evolution began with a connected graph $G_0$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "References:\n",
      "* Barab\u00e1si, <NAME>., & <NAME>. (1999). \"_Emergence of scaling in random networks._\" Science, 286(5439), 509-512.\n",
      "Chicago\t"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Bayesian preferential attachment\n",
      "Nevertheless it is easy to patch this \"vagueness\" in the model.\n",
      "Let's imagine for a moment that for some observed graph $G = (V,E)$ we are faced with the problem of estimating the probabilities ${(\\theta_v)}_{v\\in W}$ of getting attached to a node in some $W\\subseteq V$. Then the vertex degrees of the nodes in $W$ can be regarded as the evidence of attachement, whence the likelihood function should be given by:  \n",
      "\n",
      "$$L\\Big(G; {(\\theta_v)}_{v\\in W}\\Big) = \\prod_{v\\in W} \\theta_v^{\\delta(v)}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Brief theory\n",
      "For this task of parameter estimation, let's employ Bayesian approach instead of the common MLE, since the latter tends to overfit the model because it fails to incorporate the uncertainty of the parameters.\n",
      "\n",
      "Cosider a prior $\\pi(\\cdot)$ on the space of possible probability vectors ${(\\theta_v)}_{v\\in W}$ -- a $(n-1)$-dimensional simplex $S^\\circ_{n-1}\\subseteq {[0,1]}^n$ where $n=\\big\\rvert W\\big\\lvert$. Using Bayes' theorem, the posterior distribution of ${(\\theta_v)}_{v\\in W}$ given the observed data is  \n",
      "\n",
      "$$p\\Big(\\Big.{(\\theta_v)}_{v\\in W}\\,\\big.\\big\\rvert\\,G\\Big.\\Big)\n",
      "\\propto p\\Big(\\Big.G\\,\\big.\\big\\rvert\\,{(\\theta_v)}_{v\\in W}\\Big.\\Big) \\pi\\Big({(\\theta_v)}_{v\\in W}\\Big)\n",
      "= L\\Big(G; {(\\theta_v)}_{v\\in W}\\Big) \\pi\\Big({(\\theta_v)}_{v\\in W}\\Big)$$\n",
      "\n",
      "In general it is convenient to pick a prior $\\pi(\\cdot)$ from an \"**eigen-family**\" of this prior-to-posterior tranfromation. Such prior families are know nas the **conjugate prior families**. For this particular kind of likelihood such a family is the **Dirichlet distribution** given by  \n",
      "\n",
      "$$\\text{Dir}\\big({(\\alpha_v)}_{v\\in W}\\big) = \\frac{\\Gamma\\big( \\sum_{v\\in W} \\alpha_v \\big)}{ \\prod_{v\\in W} \\Gamma\\big(\\alpha_v \\big)} \\prod_{v\\in W} \\theta_v^{\\alpha_v-1}$$\n",
      "\n",
      "It is easy to demonstrate that the posterior is $\\text{Dir}\\big({(\\alpha_v)}_{v\\in W} + {(\\delta(v))}_{v\\in W}\\big)$. This implies that the posterior probability of a random vertex being connected to a node $v\\in W$ is given by  \n",
      "$$p\\Big(\\Big. v\\,\\big.\\big\\rvert\\,G\\Big.\\Big)\n",
      "= \\int p\\Big(\\Big. v\\,\\big.\\big\\rvert\\,{(\\theta_v)}_{v\\in W}\\Big.\\Big) p\\Big(\\Big. {(\\theta_v)}_{v\\in W}\\, \\big.\\big\\rvert\\,G\\Big.\\Big) d\\mathbf{\\theta} = \\ldots $$  \n",
      "\n",
      "which in our case resolves into  \n",
      "$$\\ldots = \\int \\theta_v \\text{Dir}\\big({(\\alpha_v)}_{v\\in W} + {(\\delta(v))}_{v\\in W}\\big) d\\mathbf{\\theta}\n",
      "= \\frac{\\alpha_v + \\delta(v)}{\\sum_{u\\in W} \\alpha_u + \\delta(u)}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Final probability vector\n",
      "Since there is no reason to assume any specific prior distribution of the values of the attachment probability vectors, let's use an uninformative prior, which assigna equal probability of every vector in the simplex. In the case of a Dirichlet family, this sets the hyperparameter ${(\\alpha_v)}_{v\\in W}$ of the prior to a vector of ones: $\\alpha_v=1$ for each $v\\in W$.\n",
      "\n",
      "Therefore the attachment probability distribution in $G_n$ is given by  \n",
      "\n",
      "$$p\\Big(\\Big. w'\\leadsto v\\,\\big.\\big\\rvert\\,G_n\\Big.\\Big) = \\frac{\\delta_n(v)+1}{\\big\\lvert V_n\\big\\rvert + \\sum_{u\\in S}\\delta_n(u)}$$  \n",
      "\n",
      "Listed below is the procedure, which computes the required distribution for a given graph $G$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## Bayesian valency distribution. Funny thing is that this correction\n",
      "##  turns out to is\n",
      "def ba_bayes_valency( G, nbunch = None ) :\n",
      "## Get the degree of nodes\n",
      "    deg = G.degree( nbunch )\n",
      "## Get the adjusted sum of vertex degrees\n",
      "    nbunch = nx.nodes( G ) if nbunch is None else nbunch\n",
      "    deg_sum = float( sum( deg.values( ) ) + len( nbunch ) )\n",
      "## Compute the posterior distribution of vertex valency\n",
      "    return np.fromiter( ( ( deg[ n ] + 1.0 ) / deg_sum for n in nbunch ),\n",
      "        count = len( nbunch ), dtype = np.float )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bayesian smoothing reduces the zero-degree penalty on the isolated vertices, thus enabling vertex discovery during the early stages of evolution. Notable thing is that for an uninformative prior the posterior distribution given $G$ is equivalent to the original attachment distribution for the same graph with an extra loop at every vertex."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### The long run distribution\n",
      "Using the heuristic \"mean field\" approach it is possible to derive an approximation to the digree distribution un the long run.\n",
      "\n",
      "Suppose vertices are born at intervals $\\frac{1}{n}$ with $m$ edges. Then at time $t\\in [0,1]$ there are $nt$ vertices and the sum of all node degrees is $2mnt$. Across the span of an infinitesimal interval $[t, t+dt]$ extra $n dt$ vertices are born, each with $m$ edges.\n",
      "\n",
      "Suppose a vertex $v$ is born at time $s$ with $m$ edges and at time $t>s$ this node $v$ has degree $k_s(t)$.  Then during $[t, t+dt]$ the vertex $v$ gets an incoming edge with probability $n dt \\frac{k_s(t)+1}{2mnt + nt} = \\frac{k_s(t)+1}{(2m+1)t} dt$. Ov average thus $v$ gets $\\frac{m}{2m+1} \\frac{k_s(t)+1}{t}dt$ edges, which is by how much its degree increses :  \n",
      "\n",
      "$$ \\frac{k_s(t+dt)-k_s(t)}{dt} = \\frac{m}{2m+1} \\frac{k_s(t)+1}{t}$$\n",
      "\n",
      "Note that the degree depends on the time thevertex was born and not on the vertex itself."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the continuous time limit as $n\\to \\infty$ the following differential equation approximates the dynamics of degree of the vertex $v$ :  \n",
      "\n",
      "$$\\dot{k}_s(t) = \\frac{m}{2m+1} \\frac{k_s(t)+1}{t}$$\n",
      "\n",
      "with initial value at time of birth being $m$ : $k_s(t=s) = m$.\n",
      "\n",
      "The solution to this ordinary differential equation for these initial conditions is given by  \n",
      "$$\\delta_s(t) = (m+1) \\Big(\\frac{t}{s}\\Big)^\\frac{m}{2m+1}$$\n",
      "Put $\\alpha = \\frac{m}{2m+1}$ and note that $\\frac{1}{\\alpha} = 2 + \\frac{1}{m}$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the continuous time limit the probability that at time $t$ a vertex born at $s$ has degree less than $k$ is given by \n",
      "$$F(k) = \\mathbb{P}(\\delta_s(t)\\leq k) = \\mathbb{P}\\Bigg(t \\Big(\\frac{m+1}{k+1}\\Big)^\\frac{1}{\\alpha} \\leq s\\Bigg) = \\frac{t - t \\Big(\\frac{m+1}{k+1}\\Big)^\\frac{1}{\\alpha}}{t} = 1 - \\Big(\\frac{m+1}{k+1}\\Big)^\\frac{1}{\\alpha}$$\n",
      "\n",
      "The distribution function is given by \n",
      "$$\\frac{d}{d k} F(k) = \\frac{1}{\\alpha} \\Big( k+1 \\Big)^{-\\frac{1}{\\alpha}-1} \\Big( m+1 \\Big)^\\frac{1}{\\alpha} = \\Big( k+1 \\Big)^{-3-\\frac{1}{m}} \\Big(2+\\frac{1}{m}\\Big) \\Big( m+1 \\Big)^{2+\\frac{1}{m}}$$\n",
      "\n",
      "Thus means that the limiting degree distribution of the Barabasi-Albert model with Bayesian valency behaves approximately like the power law with similar exponent as in the original model.\n",
      "\n",
      "For a more formally rigorous treatment of the original model see _**<NAME>., (2001)**_."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "References:\n",
      "* Bollob\u0<NAME>., <NAME>., <NAME>., & Tusn\u0<NAME>. (2001). \"_The degree sequence of a scale\u2010free random graph process._\" Random Structures & Algorithms, 18(3), 279-290."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Model evolution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### _Complete model_\n",
      "A has been mentioned above the basic growth step in the complete Barabasi-Albert consisted of two phases:\n",
      " 1. growth;\n",
      " 2. attachment according to valency distribution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def barabasi_albert( n, m = 10, callback = None, avclu = False, avlen = False ) :\n",
      "## Make sure the callback is a proper function\n",
      "\tcallback_fn = callback if callable( callback ) else lambda G, _, __: None\n",
      "## Barabasi-Albert evolution always begins with the isolated graph\n",
      "\tG = nx.empty_graph( m )\n",
      "## To speed up the computations keep track of the node degree locally,\n",
      "##  but instead of using a dictionary, store them in a linear array.\n",
      "\tdegree = np.zeros( n + m, dtype = np.int )\n",
      "\tdegree_half_sum = 0\n",
      "## Local adjacency matrix : should be made sparse\n",
      "\tadjacency = np.zeros( ( n + m, n + m ) )\n",
      "## Shortest path matrix for the modified Floyd Warshall\n",
      "\tpi = 1 - ml.identity( n + m, np.float )\n",
      "\tpi[ pi == 1 ] = np.inf\n",
      "\tnpmin = np.minimum\n",
      "## Local tracking of the number of triangles for the clustering coefficient\n",
      "\ttriangles = np.zeros( n + m, dtype = np.int )\n",
      "## Note that the graph itself is used just to keep track of the edges.\n",
      "\tresult = list( )\n",
      "## Create an array of average path lengths\n",
      "\tavg_path_ln = np.full( n + m, np.inf, dtype = np.float )\n",
      "## Setup the average clustering metric storage\n",
      "\tavg_tri_clu = np.full( n + m, np.inf, dtype = np.float )\n",
      "## Age-degree: the dependence of the average node degree on the age of a vertex\n",
      "\tavg_age_deg = np.zeros( n + 1, dtype = np.float )\n",
      "## Identify vertices by consequitive integers. \n",
      "\tt = G.number_of_nodes( )\n",
      "## In the original Barabasi-Albert model the first new vertex\n",
      "##  is attached to all m vertices, thus making the graph connected.\n",
      "\tU = range( m )\n",
      "\twhile t < n + m :\n",
      "## Add pending edges\n",
      "\t\tG.add_edges_from( ( t, u ) for u in U )    \n",
      "## maintain a local copy of the adjacency matrix\n",
      "\t\tadjacency[ t, U ] = adjacency[ U, t ] = 1\n",
      "## Record the vertex degree\n",
      "\t\tdegree[ U ] += 1\n",
      "\t\tdegree[ t ] += len( U )\n",
      "\t\tdegree_half_sum += len( U )\n",
      "## Create a new vertex using stack allocation -- on top of others.\n",
      "\t\tt = G.number_of_nodes( )\n",
      "## Run the Floyd-Warshall algorithm\n",
      "\t\tif avlen :\n",
      "## Add the shortest paths from the affected vertices to the new node\n",
      "\t\t\tpi[ U, t-1 ] = pi[ t-1, U ] = 1\n",
      "## Compute the shortest path from the exisitng nodes to the new one\n",
      "\t\t\tfor u in U :\n",
      "\t\t\t\tpi[ :t, t-1 ] = npmin( pi[ :t, t-1 ], pi[ :t, u ] + pi[ u, t-1 ] )\n",
      "## Restore the distance symmetry due to bidirectionality\n",
      "\t\t\tpi[ t-1, :t ] = pi[ :t, t-1 ].T\n",
      "## Compare the shortest paths through the new vertex with the current paths\n",
      "##\tAt this point the matrix pi has the following property:\n",
      "##\t * for nodes $i,j\\in V$ the value $\\pi_{ij}$ is the shortest path\n",
      "##\t   through nodes other than $\\omega$;\n",
      "##\t * for any $x\\in V$ the element $\\pi_{x\\omega}$ (and $\\pi_{\\omega x}) is\n",
      "##\t   the length of the shortest path from $i$ to the new node $\\omega$ with\n",
      "##\t   intermediate vertices exclusively in $V$.\n",
      "\t\t\tpi[ :t, :t ] = npmin( pi[ :t, :t ], pi[ :t, t-1 ] + pi[ t-1, :t ] )\n",
      "## Numerically confirm the correctness shown elsewhere in this notebook\n",
      "##\t\t\tassert( np.all( nx.floyd_warshall_numpy( G ) == pi[ :t, :t ] ) )\n",
      "## Compute the average path length\n",
      "##             avg_path_ln[ t-1 ] = nx.average_shortest_path_length( G )\n",
      "\t\t\tavg_path_ln[ t-1 ] = np.mean( pi[ :t, :t ], dtype = np.float64 ) * t / ( t - 1.0 )\n",
      "## Compute the average clustering coefficient\n",
      "\t\tif avclu :\n",
      "\t\t\ttri_num = 0\n",
      "\t\t\tfor u in U :\n",
      "## Count the number of edges from u to U, since they are responsible for the new triangles containing t\n",
      "\t\t\t\tX = np.sum( adjacency[ u, U ] )\n",
      "## Update the number of the current endpoint\n",
      "\t\t\t\ttriangles[ u ] += X\n",
      "\t\t\t\ttri_num += X\n",
      "## Every new triangle could have been created only via the edges incident on t\n",
      "##  The loop above counts the number of trinagles containing t twice (once per each edge of a triangle)\n",
      "\t\t\ttriangles[ t-1 ] += tri_num // 2\n",
      "## Compute the clustering coefficinet of each node\n",
      "\t\t\tmax_clust = degree[ :t ] * ( degree[ :t ] - 1 )\n",
      "\t\t\tclust = np.where( max_clust != 0, ( 2.0 * triangles[ :t ] ) / max_clust, 0 )\n",
      "## And the average\n",
      "\t\t\tavg_tri_clu[ t-1 ] = np.mean( clust, dtype = np.float64 )\n",
      "# \t\t\tcc = nx.average_clustering( G )\n",
      "# \t\t\tavg_tri_clu[ t-1 ] = np.abs( ( np.mean( clust, dtype = np.float64 ) - cc ) / ( cc + 1e-12 ) )\n",
      "## Compute the average node degree of created vertices of a given age.\n",
      "## * A[0:t-m0] -- current mean degrees of vertices of age 0:t-m0 (right ends not included)\n",
      "## * D[0:m0] -- the current degrees of the oldest (initial) m0 vertices\n",
      "## * D[m0:t] -- the degrees of nodes created at iterations 0:t-m0 respectively\n",
      "\t\tavg_age_deg[ :t-m ] += ( degree[ t-1 :m-1 :-1 ] - avg_age_deg[ :t-m ]\n",
      "\t\t\t) / np.arange( t-m, 0, -1, dtype = np.float )\n",
      "## Now incorporate the degree of the initial $m$ nodes, which are the oldest ones.\n",
      "\t\tavg_age_deg[ t-m ] = np.mean( degree[ :m ], dtype = np.float64 )\n",
      "## The evolution step is complete, call the required function\n",
      "\t\tresult.append( callback_fn( G, degree[ :t ], pi[ :t, :t ] ) )\n",
      "## The initial nodes contribute to the average node degree of the next generation of vertices\n",
      "## The nodes which the new vertex is attached to are picked\n",
      "##  according to the preferential attachment distribution.\n",
      "\t\tprobabilities = degree[ :t ] / ( 2.0 * degree_half_sum )\n",
      "\t\tU = rnd.choice( t, size = m, p = probabilities, replace = False )\n",
      "#     assert( np.all( nx.floyd_warshall_numpy( G ) == pi ) )\n",
      "\treturn G, result, avg_age_deg, avg_tri_clu, avg_path_ln\n",
      "## https://networkx.github.io/documentation/latest/reference/classes.graph.html\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## %lprun -f barabasi_albert barabasi_albert( 1000, m = 2, avlen = True, avclu = True )\n",
      "# G1, R1, A1, C1, L1 = barabasi_albert( 1000, m = 10, avlen = True, avclu = True )\n",
      "# plt.loglog( C1 )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### _Model A_\n",
      "In **model A**, the first restricted verison of the evolution, the network grows with time but the newly added nodes are attached to the existing ones uniformly at random.\n",
      "The code below implements the basic steps of both the Bayesian and the first restricted model."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## This evolution step is similar to the original BA step, except it\n",
      "##  relaxes the requirement that m_0 be equal to m, and it implemets\n",
      "##  the egalitarian attachment and bayesian valency.\n",
      "def barabasi_albert_model_a( n, m0 = 10, m = 5, bayes = False, callback = None, avclu = False, avlen = False ) :\n",
      "\tcallback_fn = callback if callable( callback ) else lambda G, _, __: None\n",
      "\tresult = list( )\n",
      "\tG = nx.empty_graph( m0 )\n",
      "## Local degree and shortest path tracking\n",
      "\tdegree = np.zeros( n + m0, dtype = np.int )\n",
      "\tdegree_half_sum = 0\n",
      "\tadjacency = np.zeros( ( n + m0, n + m0 ) )\n",
      "\tpi = 1 - ml.identity( n + m0, np.float )\n",
      "\tpi[ pi == 1 ] = np.inf\n",
      "\tnpmin = np.minimum\n",
      "\ttriangles = np.zeros( n + m0, dtype = np.int )\n",
      "## Average path length and clustering\n",
      "\tavg_path_ln = np.full( n + m0, np.inf, dtype = np.float )\n",
      "\tavg_tri_clu = np.full( n + m0, np.inf, dtype = np.float )\n",
      "\tavg_age_deg = np.zeros( n + 1, dtype = np.float )\n",
      "## Begin the simulation\n",
      "\tt = G.number_of_nodes( )\n",
      "\twhile t < n + m0 :\n",
      "## Select either the bayesian valency distribution, or uniform attachment\n",
      "\t\ttheta = None if not bayes else ( degree[ :t ] + 1.0 ) / ( 2.0 * degree_half_sum + t )\n",
      "\t\tU = rnd.choice( t, size = m, p = theta, replace = False )\n",
      "\t\tG.add_edges_from( ( t, u ) for u in U )\n",
      "\t\tadjacency[ t, U ] = adjacency[ U, t ] = 1\n",
      "## Accumulate degree information\n",
      "\t\tdegree[ U ] += 1\n",
      "\t\tdegree[ t ] += len( U )\n",
      "\t\tdegree_half_sum += len( U )\n",
      "## Preallocate the next vertex\n",
      "\t\tt = G.number_of_nodes( )\n",
      "## Compute the shortest paths\n",
      "\t\tif avlen :\n",
      "\t\t\tpi[ U, t-1 ] = pi[ t-1, U ] = 1\n",
      "\t\t\tfor u in U :\n",
      "\t\t\t\tpi[ :t, t-1 ] = npmin( pi[ :t, t-1 ], pi[ :t, u ] + pi[ u, t-1 ] )\n",
      "\t\t\tpi[ t-1, :t ] = pi[ :t, t-1 ].T\n",
      "\t\t\tpi[ :t, :t ] = npmin( pi[ :t, :t ], pi[ :t, t-1 ] + pi[ t-1, :t ] )\n",
      "\t\t\tavg_path_ln[ t-1 ] = np.mean( pi[ :t, :t ], dtype = np.float64 ) * t / ( t - 1.0 )\n",
      "## Get the average clustering coefficient\n",
      "\t\tif avclu :\n",
      "\t\t\ttri_num = 0\n",
      "\t\t\tfor u in U :\n",
      "\t\t\t\tX = np.sum( adjacency[ u, U ] )\n",
      "\t\t\t\ttriangles[ u ] += X\n",
      "\t\t\t\ttri_num += X\n",
      "\t\t\ttriangles[ t-1 ] += tri_num // 2\n",
      "\t\t\tmax_clust = degree[ :t ] * ( degree[ :t ] - 1 )\n",
      "\t\t\tclust = np.where( max_clust != 0, ( 2.0 * triangles[ :t ] ) / max_clust, 0 )\n",
      "\t\t\tavg_tri_clu[ t-1 ] = np.mean( clust, dtype = np.float64 )\n",
      "## Finalize this step's statistics\n",
      "\t\tavg_age_deg[ :t-m0 ] += ( degree[ t-1 :m0-1 :-1 ] - avg_age_deg[ :t-m0 ]\n",
      "\t\t\t) / np.arange( t-m0, 0, -1, dtype = np.float )\n",
      "\t\tavg_age_deg[ t-m0 ] = np.mean( degree[ :m0 ], dtype = np.float64 )\n",
      "\t\tresult.append( callback_fn( G, degree[ :t ], pi[ :t, :t ] ) )\n",
      "\treturn G, result, avg_age_deg, avg_tri_clu, avg_path_ln\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "## %lprun -f barabasi_albert_model_a barabasi_albert_model_a( 300, m = 1, avlen = True, avclu = False )\n",
      "# G1, R1, A1, C1, L1 = barabasi_albert_model_a( 1000, m0 = 10, m=10, avlen = True, avclu = True )\n",
      "# plt.loglog( C1 )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### _Model B_\n",
      "**Model B** retains the preferential attachment phase, but foregoes the growth step. Instead, a random existing node is picked and is linked with other nodes, preferentially chosing vertices with high degree. Initally this model begins its evolution from a set of isolated nodes.\n",
      "\n",
      "For this model to work as intended the only option is to use the Bayesian valency distribution, since the originally proposed probabilities handle poorly the isolated nodes, as it has been alread stated. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "###### A couple of nuances\n",
      "Let $w'\\in V$ be a randomly picked vertex. Since multiple connections and loops are forbidden in the model, the preferential attachment of a vertex $w'$ must use conditional node valency distribution. If a subset $S\\subseteq V$, $w'\\in S$, is **prohibited**, then for any vertex $v\\in V\\setminus S$  \n",
      "\n",
      "$$p\\Big(\\Big. w'\\leadsto v\\,\\big.\\big\\rvert\\, w'\\not\\!\\leadsto S,\\,G \\Big.\\Big)\n",
      "= \\frac{p\\Big(\\Big. w'\\leadsto v\\,\\big.\\big\\rvert G \\Big.\\Big)}{p\\Big(\\Big. w'\\not\\!\\leadsto S\\,\\big.\\big\\rvert G \\Big.\\Big)}\n",
      "= \\frac{\\delta_n(v)+1}{\\big\\lvert V\\setminus S\\big\\rvert + \\sum_{u\\in V\\setminus S}\\delta_n(u)} $$  \n",
      "\n",
      "In particular, the prohibited set for a vertex $w'$ is itself and the set of its immediate neighbours:  \n",
      "\n",
      "$$S = \\big\\{w'\\big\\}\\cup N_G(w')$$\n",
      "\n",
      "Another subtle point is the connection strategy when $\\big\\lvert N_G(w') \\big\\rvert\\geq n-1-m$. If $\\big\\lvert N_G(w') \\big\\rvert < n-1$, the model forcefully connects $w'$ to all of the remaining vertices. If however $w'$ is saturated, then the attachment step for this node is skipped as this vertex is just unable to form new connections."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "## Another relaxation of the complete Barabasi-Albert model, which eliminates\n",
      "##  the growth phase. This restriction uses the Bayesian valency distribution.\n",
      "##  when performing preferential attachment.\n",
      "def barabasi_albert_model_b( n, m = 5, callback = None, avclu = False, avlen = False ) :\n",
      "\tcallback_fn = callback if callable( callback ) else lambda G, _, __: None\n",
      "\tG = nx.empty_graph( n )\n",
      "## Local node degree cache\n",
      "\tdegree = np.zeros( n, dtype = np.int64 )\n",
      "\tdegree_sum = 0\n",
      "\tadjacency = np.zeros( ( n, n ), dtype = np.int )\n",
      "## Shortest path matrix for the modified Floyd Warshall\n",
      "\tpi = 1 - ml.identity( n, np.float )\n",
      "\tpi[ pi == 1 ] = np.inf\n",
      "\tnpmin = np.minimum\n",
      "\ttriangles = np.zeros( n, dtype = np.int )\n",
      "## Average path length and clustering\n",
      "\tavg_tri_clu = np.full( n, np.inf, dtype = np.float )\n",
      "\tavg_path_ln = np.full( n, np.inf, dtype = np.float )\n",
      "\tavg_age_deg = np.zeros( n, dtype = np.float )\n",
      "## Begin the simulation\n",
      "\tt = 0\n",
      "\tresult = list( )\n",
      "\twhile t < n :\n",
      "## \"Pick an existing one at random\" (cf. p.551 c.3 L 9 of [Barabasi, Albert; 1999])\n",
      "##  ... and get its neighbours.\n",
      "\t\tv = rnd.randint( n, size = 1 )[ 0 ]\n",
      "\t\tS = nx.neighbors( G, v )\n",
      "## If it has more neighbours than n-1-m, connect it to the remaining nodes\n",
      "##  and continue. If the node is staturated, ignore it.\n",
      "\t\tm_prime = min( m, n - 1 - len( S ) )\n",
      "\t\tif m_prime > 0 :\n",
      "## Get the correction terms due to conditioning (see the formula).\n",
      "\t\t\tc_term = degree[ v ] + 1.0 + np.sum( degree[ S ] ) + len( S )\n",
      "## Compute the conditional valency of each vertex in the current configuration\n",
      "\t\t\ttheta = ( degree + 1.0 ) / ( 0.0 + degree_sum + n - c_term )\n",
      "## Zero out the prohibited vertices, so that they'll never be picked.\n",
      "\t\t\ttheta[ v ] = 0 ; theta[ S ] = 0\n",
      "## Pick $m'$ other nodes to preferentially link the chosen one to.\n",
      "\t\t\tU = rnd.choice( n, size = m_prime, p = theta, replace = False )\n",
      "\t\t\tG.add_edges_from( ( v, u ) for u in U )\n",
      "\t\t\tadjacency[ v, U ] = adjacency[ U, v ] = 1\n",
      "\t\t\tdegree[ U ] += 1\n",
      "\t\t\tdegree[ v ] += m_prime\n",
      "\t\t\tdegree_sum += m_prime * 2 ;\n",
      "## Compute the shortest paths\n",
      "\t\tif avlen :\n",
      "\t\t\tif m_prime > 0 :\n",
      "\t\t\t\tif False :\n",
      "## Incrementally update the shortest path matrix by every new bidirectional edge\n",
      "\t\t\t\t\tfor u in U :\n",
      "## Compute the path through this edge\n",
      "\t\t\t\t\t\tpath = pi[ :, u ] + 1.0 + pi[ v, : ]\n",
      "## Choose the optimal edge direction and whether the new path is the shortest\n",
      "\t\t\t\t\t\tpi = npmin( pi, npmin( path, path.T ) )\n",
      "\t\t\t\telse :\n",
      "\t\t\t\t\tpi[ U, v ] = pi[ v, U ] = 1\n",
      "## Construct the shortest paths thorugh each affected vertex U wit\n",
      "##  intermediate nodes other than v. Before the addition of edges\n",
      "##  the shortest path from any node to v could not have included\n",
      "##  the vertex u in U as the last intermediate node, since there was no\n",
      "##  edge between u and v prior to this generation. If it included it\n",
      "##  as anything farther from v than one hop, then the shortest path has\n",
      "##  to be updated. Also any prior path s~u through v with a new edge u-v\n",
      "##  would create a cycle and be longer that the existing s~v (shorter that s~v~u~v).\n",
      "## Therefore this step wold coorrectly compute the length of the shortest\n",
      "##  paths to v\n",
      "\t\t\t\t\tfor u in U :\n",
      "\t\t\t\t\t\tpi[ :, v ] = npmin( pi[ :, v ], pi[ :, u ] + pi[ u, v ] )\n",
      "\t\t\t\t\tpi[ v, : ] = pi[ :, v ].T\n",
      "## See if paths x~v~y are shorter than x~y without v\n",
      "\t\t\t\t\tpi = npmin( pi, pi[ :, v ] + pi[ v, : ] )\n",
      "# \t\t\t\tassert( np.all( nx.floyd_warshall_numpy( G ) == pi ) )\n",
      "\t\t\tavg_path_ln[ t ] = np.mean( pi ) * n / ( n - 1.0 ) if m_prime > 0 else avg_path_ln[ t-1 ]\n",
      "## Get the average clustering coefficient\n",
      "\t\tif avclu :\n",
      "\t\t\ttri_num = 0\n",
      "\t\t\tfor u in G[ v ] :\n",
      "\t\t\t\tX = np.sum( adjacency[ u, U ] )\n",
      "\t\t\t\ttriangles[ u ] += X\n",
      "\t\t\t\tif u in U :\n",
      "\t\t\t\t\ttri_num += X\n",
      "\t\t\t\telse :\n",
      "\t\t\t\t\ttriangles[ v ] += X\n",
      "## it is necessary to inform the vertices adjacent to a preexisting neighbour\n",
      "## that the edge (v,u) induced a new triangle through a newly added edge.\n",
      "\t\t\t\t\ttriangles[ U ] += adjacency[ U, u ]\n",
      "\t\t\ttriangles[ v ] += tri_num // 2\n",
      "\t\t\tmax_clust = degree * ( degree - 1 )\n",
      "\t\t\tclust = np.where( max_clust != 0, ( 2.0 * triangles ) / max_clust, 0 )\n",
      "\t\t\tavg_tri_clu[ t ] = np.mean( clust, dtype = np.float64 ) if m_prime > 0 else avg_tri_clu[ t-1 ]\n",
      "## In this model all nodes are the same age, which means that the degree-age\n",
      "##  relationship is just the average degree with respect to the generation.\n",
      "\t\tavg_age_deg[ t ] = np.mean( degree ) if m_prime > 0 else avg_age_deg[ t-1 ]\n",
      "\t\tresult.append( callback_fn( G, degree, pi ) )\n",
      "\t\tt += 1\n",
      "\treturn G, result, avg_age_deg, avg_tri_clu, avg_path_ln\n"
     ],
     "language": "python",
{'fraction_non_alphanumeric': 0.1840829804798807, 'fraction_numerical': 0.012516311150169103, 'mean_word_length': 2.32593773526416, 'pattern_counts': {'":': 308, '<': 71, '<?xml version=': 0, '>': 79, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2773144', 'n_tokens_mistral': 27359, 'n_tokens_neox': 26136, 'n_words': 8450}