{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "density_circle.ipynb\n", "\n", "Discussion: Nonuniform sampling = density function; circular region\n", "\n", "Licensing: This code is distributed under the GNU LGPL license.\n", " \n", "Modified: 08 November 2016\n", "\n", "Author: John Burkardt, Lukas Bystricky" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using matplotlib backend: agg\n" ] } ], "source": [ "# Import necessary libraries and set plot option\n", "%matplotlib\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import math\n", "import scipy.spatial as spatial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# More about Density Functions #\n", "\n", "In the previous notebook, we looked at how density functions affect the \n", "computation of centers of mass and the CVT iteration, in a triangle and\n", "a square. We used uniform sampling, and a density function rho(x,y).\n", "It turns out that we can combine sampling and nonuniform density into\n", "a single operation, a nonuniform sampling method. This simply picks\n", "some points more often than others. If we can find an appropriate\n", "nonuniform sampling method, then the CVT iteration is simplified.\n", "\n", "We will look at:\n", "* using nonuniform sampling in a circle as another way to enforce a density.\n", "* displaying plots of the data during a CVT iteration.\n", "* adding the boundary or Voronoi diagram to such plots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sample the circle nonuniformly #\n", "\n", "There are actually two ways of computing a CVT with a density function.\n", "\n", "Way number 1:\n", "* Uniformly randomly sample the region;\n", "* Compute the density function at each sample point;\n", "* Modify calculations using weighted averaging\n", " \n", "Way number 2:\n", "* Modify the sampling technique for \"weighted\" random sampling.\n", "* The rest of the CVT algorithm behaves as though it was a uniform problem.\n", " \n", "Let's try to look at the circle this way. We know how to sample uniformly.\n", "Instead, let's do a nonuniform sampling, like this:\n", "```python\n", "def circle_sample_nonuniform ( m )\n", " import numpy as np\n", " s = np.zeros ( [ m, 2 ] )\n", " r = np.random.random ( m )\n", " t = 2.0 * np.pi * np.random.random ( m )\n", " s[:,0] = r * np.cos ( t )\n", " s[:,1] = r * np.sin ( t )\n", " return s\n", "```" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Function to sample the circle nonuniformly\n", "\n", "def circle_sample_nonuniform ( m ):\n", " import numpy as np\n", " s = np.zeros ( [ m, 2 ] )\n", " r = np.random.random ( m )\n", " t = 2.0 * np.pi * np.random.random ( m )\n", " s[:,0] = r * np.cos ( t )\n", " s[:,1] = r * np.sin ( t )\n", " return s\n", "#\n", "# Test the function. \n", "#\n", "m = 500\n", "s = circle_sample_nonuniform ( m )\n", "plt.plot ( s[:,0], s[:,1], 'bo' )\n", "plt.axis ( 'Equal')\n", "plt.axis ( [ -1.0, 1.0, -1.0, 1.0 ] )\n", "#\n", "# There seem to be various ways in MATPLOTLIB to draw a circle.\n", "# I can't get them to work. So I did it this way.\n", "#\n", "c = np.zeros ( [65, 2 ] )\n", "t = np.linspace ( 0, 2.0*np.pi, 65)\n", "c[:,0] = np.cos ( t )\n", "c[:,1] = np.sin ( t )\n", "plt.plot ( c[:,0], c[:,1], 'r' )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# CVT step for Circle, using nonuniform sampling #\n", "\n", "Once we have our nonuniform sampling method, the CVT step\n", "can be written as though we were solving a uniform problem\n", "with no weighting. \n", "\n", "You can more or less copy the CVT code we used for the\n", "square with uniform sampling, and then just change the name\n", "of the sampling function:\n", "```python\n", "def cvt_step_circle_nonuniform ( g, m ):\n", " import numpy as np\n", " s = circle_sample_nonuniform ( m )\n", " ***\n", " return c\n", "```" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0 0.0000 0.0000 -0.0672 -0.2790\n", " 1 0.2000 0.2000 0.4716 0.2041\n", " 2 -0.1000 0.3000 -0.3314 0.4111\n" ] } ], "source": [ "# Code for CVT step for Circle, using nonuniform sampling.\n", "#\n", "def cvt_step_circle_nonuniform ( g, m ):\n", " import numpy as np\n", " s = circle_sample_nonuniform ( m )\n", " n = g.shape[0]\n", " ni = np.zeros ( n )\n", " c = np.zeros ( [ n, 2 ] )\n", " for i in range ( 0, m ):\n", " k = -1\n", " d = np.Inf\n", " for j in range ( 0, n ):\n", " dj = np.linalg.norm ( s[i,:] - g[j,:] )\n", " if ( dj < d ):\n", " d = dj\n", " k = j\n", " ni[k] = ni[k] + 1\n", " c[k,:] = c[k,:] + s[i,:]\n", " \n", " for i in range ( 0, n ):\n", " if ( ni[i] == 0 ):\n", " c[i,:] = g[i,:]\n", " else:\n", " c[i,:] = c[i,:] / float ( ni[i] )\n", " \n", " return c\n", "\n", "#\n", "# Test the function.\n", "# We just expect C values to be near G values.\n", "#\n", "n = 3\n", "g = np.array ( [ \\\n", " [ 0.0 ,0.0 ], \\\n", " [ 0.2, 0.2 ], \\\n", " [ -0.1, 0.3 ] ] )\n", "\n", "m = 5000\n", "c = cvt_step_circle_nonuniform ( g, m )\n", "\n", "for i in range ( 0, n ):\n", " print ( '%2d %8.4f %8.4f %8.4f %8.4f' % ( i, g[i,0], g[i,1], c[i,0], c[i,1] ) )\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# CVT iteration for Circle with Nonuniform Density #\n", "\n", "As usual, the CVT step is the hard part of the coding.\n", "To do a CVT iteration, we just need to initialize the\n", "generators, and then carry out a loop that repeatedly\n", "computes centers of mass C and uses them to overwrite G.\n", "```python\n", "def cvt_circle_nonuniform ( n, m, step_num ):\n", " g = circle_sample_nonuniform ( m )\n", " for step in range ( 0, step_num ):\n", " ***\n", " return g\n", "```" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "What?\n" ] } ], "source": [ "# Code for CVT iteration for circle with nonuniform density\n", "#\n", "import matplotlib.pyplot as plt\n", "\n", "def cvt_circle_nonuniform ( n, m, step_num ):\n", " g = circle_sample_nonuniform ( n )\n", " plt.plot ( g[:,0], g[:,1], 'bo' )\n", " plt.axis ( 'Equal')\n", " plt.show ( )\n", " \n", " for step in range ( 0, step_num ):\n", " c = cvt_step_circle_nonuniform ( g, m )\n", " g = c.copy ( )\n", " \n", " plt.plot ( g[:,0], g[:,1], 'ro' )\n", " plt.axis ( 'Equal' )\n", " plt.show ( )\n", " return g\n", "\n", "#\n", "# Test the function.\n", "#\n", "n = 50\n", "m = 10000\n", "step_num = 10\n", "g = cvt_circle_nonuniform ( n, m, step_num )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Overlay the circular boundary on the plots\n", "\n", "Modify the cvt_circle_nonuniform code so that the two plots of\n", "the generators now also include a red circle, indicating the\n", "boundary of the region.\n", "```python\n", "def cvt_circle_nonuniform2 ( n, m, step_num ):\n", " t = np.linspace ( n, m, step_num ):\n", " xc = np.cos ( t )\n", " yc = np.sin ( t )\n", " ***\n", " return g\n", "```\n", "Although MATPLOTLIB apparently has ways to draw a circle, I\n", "could not get them to work for me, so I just computed 50\n", "points on a circle and plotted them." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Code for CVT iteration for circle with nonuniform density\n", "#\n", "import matplotlib.pyplot as plt\n", "\n", "def cvt_circle_nonuniform2 ( n, m, step_num ):\n", " \n", " t = np.linspace ( 0, 2.0 * np.pi, 50 )\n", " xc = np.cos ( t )\n", " yc = np.sin ( t )\n", "\n", " g = circle_sample_nonuniform ( n )\n", " plt.plot ( g[:,0], g[:,1], 'go' )\n", " plt.plot ( xc, yc, 'r-')\n", " plt.axis ( 'Equal')\n", " plt.show ( )\n", " \n", " for step in range ( 0, step_num ):\n", " c = cvt_step_circle_nonuniform ( g, m )\n", " g = c.copy ( )\n", " \n", " plt.plot ( g[:,0], g[:,1], 'bo' )\n", " plt.plot ( xc, yc, 'r-')\n", " plt.axis ( 'Equal' )\n", " plt.show ( )\n", " return g\n", "\n", "#\n", "# Test the function.\n", "#\n", "n = 50\n", "m = 10000\n", "step_num = 10\n", "g = cvt_circle_nonuniform2 ( n, m, step_num )\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Overlay the Voronoi diagram on the plots.\n", "\n", "Instead of showing the circular boundary, we could show the\n", "Voronoi diagram.\n", "```python\n", "def cvt_circle_nonuniform3 ( n, m, step_num ):\n", " g = circle_sample_nonuniform ( n )\n", " vor = spatial.Voronoi ( g )\n", " spatial.voronoi_plot_2d ( vor )\n", " ***\n", " return g\n", "```\n", "\n", "Unfortunately, the Voronoi diagram plotting routine does not\n", "make a very good judgment about the size of the plot, so that\n", "sometimes the interesting part of the plot is much too small." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Code for CVT iteration for circle with nonuniform density\n", "#\n", "import matplotlib.pyplot as plt\n", "\n", "def cvt_circle_nonuniform3 ( n, m, step_num ):\n", " g = circle_sample_nonuniform ( n )\n", " vor = spatial.Voronoi ( g )\n", " spatial.voronoi_plot_2d ( vor )\n", " plt.plot ( g[:,0], g[:,1], 'go' )\n", " plt.axis ( 'Equal')\n", " plt.show ( )\n", " \n", " for step in range ( 0, step_num ):\n", " c = cvt_step_circle_nonuniform ( g, m )\n", " g = c.copy ( )\n", " \n", " vor = spatial.Voronoi ( g )\n", " spatial.voronoi_plot_2d ( vor )\n", " plt.plot ( g[:,0], g[:,1], 'bo' )\n", " plt.axis ( 'Equal' )\n", " plt.show ( )\n", " return g\n", "\n", "#\n", "# Test the function.\n", "#\n", "n = 50\n", "m = 10000\n", "step_num = 10\n", "g = cvt_circle_nonuniform3 ( n, m, step_num )\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "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.8" } }, "nbformat": 4, "nbformat_minor": 0 }