Approximate MST Weight.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximate the MST Weight\n",
    "Source: http://www.cs.princeton.edu/courses/archive/spring04/cos598B/bib/Czumaj-Sohler-EstimMST.pdf\n",
    "\n",
    "In the paper:\n",
    "\n",
    "$P$ is a metric space (In the code, $(P,d)$ is a metric space)  we have that: $\\left|P\\right| = n$\n",
    "\n",
    "$G^{(i)}$ and is the $(1+\\epsilon)^i$-radius graph on $(P,d)$\n",
    "\n",
    "The distances are rescaled in $O(n)$ time so that $d(x,y) \\in [0,(1+\\epsilon)^r]$ where $(1+\\epsilon)^r$ is the smallest power of $1+\\epsilon$ above $\\frac{4n}{\\epsilon}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import rv_discrete \n",
    "from scipy.spatial import distance\n",
    "from math import *\n",
    "import itertools\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"degree_estimate.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Degree_Estimate(P,d,threshold,p,constant=1):\n",
    "    n = len(P)\n",
    "    clogn = constant*log(n)\n",
    "    l = 1\n",
    "    while True:\n",
    "        l = 2*l\n",
    "        N = clogn*l\n",
    "        S = np.random.choice(n,N,replace=True)\n",
    "        #number of vertices in S that are adjacent to p in G_{threshold}\n",
    "        Y = sum(1 for elt in S if (d(p,elt) < threshold))\n",
    "        if( Y > clogn or l > n):\n",
    "            return float(Y*n)/float(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src =\"clique_tree_traversal.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#This is Clique_Tree_Traversal from the paper modified to have the stopping criteria outlined in the paper.\n",
    "#threshold = (1+eps)^i\n",
    "def Modified_Clique_Tree_Traversal(P,d,p,X,threshold,epsilon,r):\n",
    "    R = [p] #representive vertices\n",
    "    E = [p] #explored vertices\n",
    "    U = range(0,p)+range(p+1,len(P)) #unexplored vertices, everything but p\n",
    "    teps = threshold*epsilon #this is rhl being a C programmer and saving a multiply\n",
    "    edges = [(d(a,b),a,b) for a in R for b in U]\n",
    "    #todo remove edges of weight larger than threshold*(1+epsilon)\n",
    "    edges.sort(reverse=True)\n",
    "    for i in xrange(len(edges) - 1, -1, -1):\n",
    "            if edges[i][0] >= threshold*(1+epsilon):\n",
    "                del edges[i]\n",
    "    while True:\n",
    "        if (len(edges) == 0 or len(E) > X or len(R) > (4*r/epsilon)):\n",
    "            return 0\n",
    "        #Since all the unexplored vertices are now explored\n",
    "        #We know that we have explored all of the vertices connected to p in the threshold\n",
    "        #perhaps we need to keep track of connected components with a d.s. datastructure?\n",
    "        if (len(U) == 0): \n",
    "            return 1\n",
    "        w = edges[-1][0] #the minimum weight edge between the representative and unexplored vertices\n",
    "        a = edges[-1][1] #the representative vertex\n",
    "        b = edges[-1][2] #the unexplored vertex\n",
    "        E.append( b) #mark b as explored\n",
    "        U.remove( b) #remove b from the unexplored\n",
    "        if( d(a,b) > teps): #if d(a,b) > teps\n",
    "            R.append(b) #append the now explored vertex b to the representative vertex list\n",
    "            edges.extend( [(d(b,h),b,h) for h in U]) #compute distance from b to the remaining U guys.\n",
    "            edges.sort(reverse=True) #merge the two sorted lists\n",
    "        #remove all edges to b in the edges list\n",
    "        for i in xrange(len(edges) - 1, -1, -1):\n",
    "            if (edges[i][2] == b) or (edges[i][0] >= (threshold*(1+epsilon))):\n",
    "                del edges[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"number_of_connected_components.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#threshold = (1+eps)^i\n",
    "def Number_Of_Connected_Components(P,d,threshold,T,epsilon,r):\n",
    "    n = len(P)\n",
    "    #Build the distribution P[X >= k] = 1/k$\n",
    "    probabilities = np.array([1.0/k - 1.0/(k+1) for k in xrange(1,n+1)])\n",
    "    distrib = rv_discrete(values=((xrange(1,n+1)),probabilities))\n",
    "    s = 0\n",
    "    chat=0\n",
    "    while (s <= T):\n",
    "        s = s+1\n",
    "        p = np.random.randint(n) #draw a vertex uniformly, independently at random.\n",
    "        X = distrib.rvs(size=1)\n",
    "        D = Degree_Estimate(P,d,threshold,p)\n",
    "        if (D < 2*X):\n",
    "            chat = chat + Modified_Clique_Tree_Traversal(P,d,p,X,threshold,epsilon,r)\n",
    "    return (float(n)/s)*chat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"metric_mst_approx.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Metric_MST_Approximation(P,orig_d,epsilon,constant=1):\n",
    "    n = len(P)\n",
    "    #the paper asks for this.. for now i've turned it off\n",
    "    #epsilon = epsilon/4.0\n",
    "    r = int(ceil(log(4*n/epsilon)/log1p(epsilon)))\n",
    "    vtx = np.random.randint(n) #draw a vertex uniformly, independently at random.\n",
    "    #by the triangle inequality d(x,y) <= 2*max_d_to_vtx\n",
    "    max_d_to_vtx = max([orig_d(P[x,:],P[vtx,:]) for x in xrange(0,n)])\n",
    "    #multiplying all the distances by this makes them in the range [0,(1+epsilon)^r]\n",
    "    rescale_factor = pow(1+epsilon,r)/(2*max_d_to_vtx)\n",
    "    #we turn the original distance into the rescaled distance\n",
    "    d = lambda x,y: orig_d(P[x,:],P[y,:])*rescale_factor    \n",
    "    c = 0\n",
    "    T = int(constant*n/pow(epsilon,6.0))\n",
    "    for i in xrange(0,r):\n",
    "        threshold = pow((1+epsilon),i)\n",
    "        print threshold\n",
    "        sys.stdout.flush()\n",
    "        cc = Number_Of_Connected_Components(P,d,threshold,T,epsilon,r)\n",
    "        print \"cc =\", cc\n",
    "        c = c + threshold*cc\n",
    "    return (n - pow(1+epsilon,r)) + epsilon*c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "P = np.random.random((10,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3670356698\n"
     ]
    }
   ],
   "source": [
    "print distance.euclidean(P[7,:],P[4,:]) #we just need to be able to do stuff like this"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "/** This is currently totally wrong, the estimates on the # of connected components are way off **/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#print Metric_MST_Approximation(P,distance.euclidean,.4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}