Diffusion+Limited+Aggregation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as plt\n",
    "import matplotlib.image as mpimg\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class DLA():\n",
    "\n",
    "    def __init__(self, size):\n",
    "        self.size = size\n",
    "        self.state = np.zeros((size, size), dtype = int)\n",
    "        self.state[size / 2, size / 2] = 1\n",
    "        \n",
    "    def printState(self):\n",
    "        plt.pyplot.imshow(1 - self.state, cmap='gray')\n",
    "\n",
    "    def getSeed(self):\n",
    "        '''\n",
    "            Returns a randomly sampled initial position\n",
    "            for a particle. \n",
    "            Returns (x, y) tuple\n",
    "        '''\n",
    "        \n",
    "        p = int(np.random.rand() / 0.25)    # select edge\n",
    "        \n",
    "        possible = [(0, np.random.randint(self.size - 1)),               # Each edge has size - 1 points\n",
    "                    (np.random.randint(self.size - 1), self.size - 1),   # Corners shouldn't get double prob\n",
    "                    (self.size - 1, np.random.randint(1, self.size)), \n",
    "                    (np.random.randint(1, self.size), 0)]\n",
    "        \n",
    "        return possible[p]\n",
    "    \n",
    "    def getAdjacentPoints(self, curr):\n",
    "        '''\n",
    "            Returns points adjacent to curr within image bounds\n",
    "            Assumption : Adjacent includes diagonal neighbors (max 8)\n",
    "                            A|A|A\n",
    "                            A|X|A\n",
    "                            A|A|A\n",
    "                 \n",
    "            Input args:\n",
    "                curr : (x, y) tuple\n",
    "            Returns:\n",
    "                List of adjacent points\n",
    "        '''\n",
    "        \n",
    "        x, y = curr\n",
    "        adjacentPoints = [(x - 1, y - 1), (x - 1, y), (x - 1, y + 1),\n",
    "                          (x, y - 1),                 (x, y + 1),\n",
    "                          (x + 1, y - 1), (x + 1, y), (x + 1, y + 1)]\n",
    "        \n",
    "        # Remove points outside the image\n",
    "        adjacentPoints = filter(lambda x : x[0] > -1 and x[0] < self.size and \\\n",
    "                                           x[1] > -1 and x[1] < self.size, adjacentPoints)\n",
    "        return adjacentPoints\n",
    "    \n",
    "    def checkIfTerminate(self, curr):\n",
    "        '''\n",
    "            Check if curr sticks in image.\n",
    "            Will happen if any adjacent block is 1\n",
    "     \n",
    "            Input args:\n",
    "                curr : (x, y) tuple\n",
    "            Returns True/False\n",
    "        '''\n",
    "        adjacentPoints = self.getAdjacentPoints(curr)\n",
    "        return any(map(lambda x : self.state[x] == 1, adjacentPoints))\n",
    "    \n",
    "    def getNextPosition(self, curr):\n",
    "        \n",
    "        '''\n",
    "            Get next point. Brownian motion is order 1 markov process\n",
    "            \n",
    "            Input args:\n",
    "                curr : (x, y) tuple\n",
    "            Returns:\n",
    "                (x, y) coordinate of next position\n",
    "        '''\n",
    "\n",
    "        adjacentPoints = self.getAdjacentPoints(curr) # List of adjacent points\n",
    "        s = np.random.randint(len(adjacentPoints))    # Get random point\n",
    "        return adjacentPoints[s]\n",
    "\n",
    "        \n",
    "    def addPoint(self, numPoints = 1):\n",
    "\n",
    "        count = 0\n",
    "        while count < numPoints:\n",
    "            \n",
    "            if (count + 1) % 10 == 0:\n",
    "                print 'added', count + 1\n",
    "            \n",
    "            curr = self.getSeed() # Get initial position\n",
    "            while not self.checkIfTerminate(curr):\n",
    "                curr = self.getNextPosition(curr)\n",
    "            self.state[curr] = 1\n",
    "\n",
    "            count += 1\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dla = DLA(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dla.addPoint(500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP4AAAD8CAYAAABXXhlaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEFtJREFUeJzt3W+oZPV9x/H3p7vZWA2Nq16Wza52tygJEkgNF6tYStGE\nWhuiDyQYQlmKsE/SxvyBRNsHoc8qhBgflMCiDUuRxHQjVSQkpBvzoE+2XqM0uqtxq0Z3WfUKmpQ8\naZd8+2DObafj3JnfmTln5pz5fl5w2TtzzznznaPf+/nNb35zriICM8vlt5ZdgJktnhvfLCE3vllC\nbnyzhNz4Zgm58c0ScuObJTRX40u6WdILkk5LuruposysXZp1AY+kHcDPgY8DZ4AngU9HxMnmyjOz\nNuycY99rgdMR8RKApO8AtwLbNv5ll10WBw4cmOMhzWySV155hbfeekvTtpun8fcBrw3dPgP8wehG\nkg4DhwGuuOIKNjY25nhIM5tkfX29aLvWJ/ci4khErEfE+traWtsPZ2YF5mn8s8DlQ7f3V/eZWcfN\n0/hPAldJOihpF3AH8FgzZZlZm2Z+jR8R5yX9JfBDYAfwDxHxXGOVmVlr5pncIyK+D3y/oVrMbEG8\ncs8sITe+WUJufLOE3PhmCbnxzRJy45sl5MY3S8iNb5aQG98sITe+WUJufLOE3PhmCbnxzRJy45sl\n5MY3S8iNb5aQG98sITe+WUJufLOE3PhmCbnxzRJy45sl5MY3S8iNb5aQG98sITe+WUJufLOE3Phm\nCbnxzRJy45sl5MY3S8iNb5aQG98sITe+WUJTG1/S5ZKekHRS0nOS7qruv0TSjyS9WP27u/1yzawJ\nJYl/HvhSRFwNXAd8VtLVwN3A8Yi4Cjhe3TazHpja+BFxLiJ+Wn3/n8ApYB9wK3C02uwocFtbRZpZ\ns2q9xpd0ALgGOAHsiYhz1Y9eB/Y0WpmtJElIWnYZ6RU3vqT3Ad8DPh8Rvxr+WUQEENvsd1jShqSN\nzc3NuYo1s2YUNb6k9zBo+oci4pHq7jck7a1+vhd4c9y+EXEkItYjYn1tba2Jmnslc8JtPffhr4hg\nkBO2TCWz+gIeBE5FxNeHfvQYcKj6/hDwaPPlmVkbdhZscwPw58DPJD1T3ffXwN8B35V0J/AL4FPt\nlNhvW+m2lfp10m7SPqOjiC6m6OhzH/6+i/VmMrXxI+Jfge3Gqjc1W46ZLUJJ4tscmkj6cXMEfUhM\np3t3ecmuWUJO/A7ZLun7kph9mHewASe+WUJO/JbVSe++J32X6y15h6TL9TfNiW+WkBvfLCEP9QvM\nMxScZd8mFv3U3W/V9fVlVFuc+GYJOfEnWHTSz/N22DKTbLuFRl1M1XG1drHOtjnxzRJy4k8wT5LV\n2bdkie60fZpW5zk3/QEka58T3ywhJ36B0VSqk1YlKT56nC58OKfkOfcx6Wd5Hqs4H+DEN0vIiV/D\noj9EU2c+oO0kmjQbPq2mLiVmSU19eHdiXk58s4Tc+GYJeahfYJahbFtD8DaGnSVD8XkmxfoyVO5b\nvfNw4psl5MSfoHQibVIKjk4UdWmia0tJ/SWjm0nH7ZMMye/EN0vIiT/BpGQfVpLiTV9jvsm/zjOp\n/pI5iz595LWphUh958Q3S8iJ34CS6+hN2naWRTnzLCMusegLiLSty4uKlsGJb5aQE78DSpe+1lW6\nf51Z/RJ9SNA+1NgmJ75ZQk78Dpi2PmDWlYGzpPYss/pdVmcE09fnOAsnvllCbnyzhDzU74AmJ9+a\nfptq1d4GK7m6UcniqD4+92FOfLOEnPg9M21BUNtLUetcbXeZk2XTHrvO4qhV5MQ3S6i48SXtkPS0\npMer2wclnZB0WtLDkna1V6aVkvS/XxHxrhFAkx/usf6qk/h3AaeGbt8L3BcRVwJvA3c2WZiZtaeo\n8SXtB/4MeKC6LeBG4Fi1yVHgtjYKtDKjKR8RjSb8uGNtd/ztRhxN8whmdqWJ/w3gy8BvqtuXAu9E\nxPnq9hlg37gdJR2WtCFpY3Nzc65izawZUxtf0ieANyPiqVkeICKORMR6RKyvra3NcggrMC5dR5N/\nlpHAdvME4443PNpoMum3q3f08Z3+5UrezrsB+KSkW4ALgN8B7gculrSzSv39wNn2yjSzJk1N/Ii4\nJyL2R8QB4A7gxxHxGeAJ4PZqs0PAo61VaWaNmud9/K8AX5R0msFr/gebKcmaNjr0rjPk324YP+54\nTSudJBw3odn2xGLf1Vq5FxE/AX5Sff8ScG3zJZlZ27xkN7F5P+/fhqY+BJTps/Wz8JJds4Sc+Nap\na8zPUkvfPyq8DE58s4Sc+DbzwpdFjw669LHfvnPimyXkxO+wtv/O26TjT9tneL+2krf07w1Mqtmj\ngfGc+GYJOfGXqPTyUG3NWtc5/jL+Wu52ST7t8mOT9i15vAyjBCe+WUJufLOEPNRfomlD5LYn90br\nGHf8OtfVb3LCcd7jlZ7bcftk4MQ3S8iJ3wHbvS21yIm0accv+WOT2x1rmeqc20yc+GYJOfE7qCSN\nFp3840x77JI5ijpJPMtzrbPYJxMnvllCTvyWlaTUPKld5zX4PEk5ev+k45W8E1Dn+LPMzG+3T/bX\n9luc+GYJOfErdZJ5kjqz4k2+Pm9qee88awqaSNk68wEl+9h4TnyzhNz4ZgmlG+rPsiS1zgKVacPd\ntq8PV+ez9eO2q/M2YVND+mnb1Bm+e4hfxolvllC6xK8zEVW6QGW7/cdt23QilXxOfrtt6+y7TF2q\nZVU48c0SSpf4bbxeLBkltPUaf5bX5E0tW/VbZ/3lxDdLKF3it3V12tLjz/v4095hmDSiqJPQ2314\npkt/Z89m58Q3Syhd4s9iWvq1nXR15gUmvY8/z3v0s1xlN/vlrbrMiW+WkBO/wLSLRixyNd40JSsP\n5/n476TH9Gv+/nDimyVU1PiSLpZ0TNLzkk5Jul7SJZJ+JOnF6t/dbRdrZs0oTfz7gR9ExIeAjwCn\ngLuB4xFxFXC8ur3SIoKIQNK7htRbX8uyVdO4Wraruy2Lfjyrb2rjS3o/8EfAgwAR8V8R8Q5wK3C0\n2uwocFtbRZpZs0oS/yCwCXxL0tOSHpB0EbAnIs5V27wO7GmryK5ZdrrD/yX8aNIP3ze6bVtJvN35\nmFTTaG3jnpe1p6TxdwIfBb4ZEdcAv2ZkWB+D/+pjO0HSYUkbkjY2NzfnrdfMGlDS+GeAMxFxorp9\njMEvgjck7QWo/n1z3M4RcSQi1iNifW1trYmabUjJ6GPZI5TheYdpI5UuzJdkMLXxI+J14DVJH6zu\nugk4CTwGHKruOwQ82kqFZta40gU8fwU8JGkX8BLwFwx+aXxX0p3AL4BPtVOi1Vn6WmehzbjjLypp\n63xseVELhLp0ftpW1PgR8QywPuZHNzVbjpktgpfsLlHppb2aSpkml/DWMWtiLvpDUCUXKlmV5che\nsmuWkBvfLCEP9Zdo2tCy6U//LXryquTPbk2yqJcms0wa9n3I78Q3S8iJ3wGTEnjcz9t8/CaTrKm3\n6Ea3aTptZzneot9qbJoT3ywhJ34HNfE6fvg4y/7bcyUpPstxZrmmX9PXAexb0m9x4psl5MRfMZNe\nr3cpleosmtlu3y2Triw8atxzb+I8deGc1uHEN0vIib/Cujzz3MTr6qaO36XzsihOfLOEnPgJdSHR\nFr1WoUQXalgUJ75ZQm58s4Q81F9hfZisaqu2Pjz3ZXLimyXkxF8RfUu4Nj7+O3zcaQt5VvE6enU4\n8c0ScuKviKYv2tEnJVcWnrRPyTarxolvlpAT31ZKkx+w6du8SR1OfLOEnPg9UPK6fd4LW9q7rfK8\niRPfLCEnfoeVpHidi1SuSlptmXQZrTYuGtr0cZfJiW+WkBvfLCEP9TugzgKSea6cuyrD1C2Leh6e\n3DOzleDE74B5/nR0SYqvWtKXaOLv7q3yeXPimyXkxO+ZVX293pa2r+bbV058s4SKGl/SFyQ9J+lZ\nSd+WdIGkg5JOSDot6WFJu9ou1mYTESudXlbf1MaXtA/4HLAeER8GdgB3APcC90XElcDbwJ1tFmpm\nzSkd6u8EflvSTuBC4BxwI3Cs+vlR4Lbmy8tN0rtmpbfu20rxcduYTTO18SPiLPA14FUGDf9L4Cng\nnYg4X212Btg3bn9JhyVtSNrY3Nxspmozm0vJUH83cCtwEPgAcBFwc+kDRMSRiFiPiPW1tbWZCzWz\n5pS8nfcx4OWI2ASQ9AhwA3CxpJ1V6u8HzrZXZk4lnwf323s2i5LX+K8C10m6UIP/u24CTgJPALdX\n2xwCHm2nRDNrWslr/BMMJvF+Cvys2ucI8BXgi5JOA5cCD7ZYZ0qjE3lOc2tK0cq9iPgq8NWRu18C\nrm28IjNrnZfsdtgsV4pdxY+QWvO8ZNcsISf+inHKWwknvllCbnyzhNz4Zgm58c0ScuObJeTGN0vI\njW+WkBvfLCE3vllCbnyzhNz4Zgm58c0ScuObJeTGN0vIjW+WkBvfLCE3vllCbnyzhNz4Zgm58c0S\ncuObJeTGN0vIjW+WkBvfLCE3vllCbnyzhNz4Zgm58c0ScuObJeTGN0vIjW+WkBvfLCE3vllCbnyz\nhBQRi3swaRP4NfDWwh50PpfRn1qhX/X2qVboT72/GxFr0zZaaOMDSNqIiPWFPuiM+lQr9KvePtUK\n/at3Gg/1zRJy45sltIzGP7KEx5xVn2qFftXbp1qhf/VOtPDX+Ga2fB7qmyW0sMaXdLOkFySdlnT3\noh63lKTLJT0h6aSk5yTdVd1/iaQfSXqx+nf3smvdImmHpKclPV7dPijpRHWOH5a0a9k1bpF0saRj\nkp6XdErS9V09t5K+UP0/8Kykb0u6oMvndhYLaXxJO4C/B/4UuBr4tKSrF/HYNZwHvhQRVwPXAZ+t\narwbOB4RVwHHq9tdcRdwauj2vcB9EXEl8DZw51KqGu9+4AcR8SHgIwzq7ty5lbQP+BywHhEfBnYA\nd9Dtc1tfRLT+BVwP/HDo9j3APYt47DlqfhT4OPACsLe6by/wwrJrq2rZz6BZbgQeB8RggcnOced8\nybW+H3iZak5p6P7OnVtgH/AacAmwszq3f9LVczvr16KG+lsnc8uZ6r5OknQAuAY4AeyJiHPVj14H\n9iyprFHfAL4M/Ka6fSnwTkScr2536RwfBDaBb1UvTR6QdBEdPLcRcRb4GvAqcA74JfAU3T23M/Hk\n3ghJ7wO+B3w+In41/LMY/Lpf+tsgkj4BvBkRTy27lkI7gY8C34yIaxgs2/5/w/oOndvdwK0Mfll9\nALgIuHmpRbVgUY1/Frh86Pb+6r5OkfQeBk3/UEQ8Ut39hqS91c/3Am8uq74hNwCflPQK8B0Gw/37\ngYsl7ay26dI5PgOciYgT1e1jDH4RdPHcfgx4OSI2I+K/gUcYnO+untuZLKrxnwSuqmZGdzGYLHls\nQY9dRJKAB4FTEfH1oR89Bhyqvj/E4LX/UkXEPRGxPyIOMDiXP46IzwBPALdXm3WiVoCIeB14TdIH\nq7tuAk7SwXPLYIh/naQLq/8ntmrt5Lmd2QInTW4Bfg78B/A3y57cGFPfHzIYav478Ez1dQuD187H\ngReBfwEuWXatI3X/MfB49f3vAf8GnAb+CXjvsusbqvP3gY3q/P4zsLur5xb4W+B54FngH4H3dvnc\nzvLllXtmCXlyzywhN75ZQm58s4Tc+GYJufHNEnLjmyXkxjdLyI1vltD/AM7K7FqrpMTSAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d186390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dla.printState()"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}