{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"density_square.ipynb\n",
"\n",
"Discussion: CVT calculations in a square with a density function.\n",
"\n",
"Licensing: This code is distributed under the GNU LGPL license.\n",
" \n",
"Modified: 07 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": [
"# Density functions #\n",
"\n",
"If we want to think of a geometric region R in the plane as a physical\n",
"object, we could imagine cutting out the shape from a sheet of\n",
"cardboard. Physical objects can have a thickness, though, and\n",
"we can apply our Voronoi techniques to some interesting problems\n",
"if we can see how to include a concept of thickness. This will\n",
"allow us to assign more importance to some parts of the region,\n",
"corresponding to population of a country, or actual height of\n",
"a surface, or thickness of a material object. \n",
"\n",
"We can do all these things by specifying a density function, which \n",
"we will call rho(x,y). We will assume that rho(x,y) is positive \n",
"in R, although we can allow it to be zero at some relatively \n",
"small portion of R. When a density function is involved, then\n",
"instead of the geometric centroid of R, we will want to compute\n",
"the center of mass of R. If rho(x,y)=1 everywhere, these \n",
"quantities are equal.\n",
"\n",
"We will look at\n",
"* how sampling can estimate the center of mass of the unit triangle;\n",
"* how to do a CVT in the square, including a density;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The center of mass #\n",
"\n",
"Give a density rho(x,y), the center of mass of a region R\n",
"is the point c whose coordinates can be computed by:\n",
"* mass = integral ( x, y in R ) rho(x,y) 1 dx dy \n",
"* c_x = integral ( x, y in R ) rho(x,y) x dx dy / mass\n",
"* c_y = integral ( x, y in R ) rho(x,y) y dx dy / mass\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"# Sample the Unit Triangle #\n",
"\n",
"Let's assume our region R is the unit triangle whose vertices\n",
"are (0,0), (1,0), (0,1). \n",
"\n",
"Let's use sampling to estimate the integrals involved in the\n",
"center of mass. \n",
"\n",
"We looked at sampling any triangle in the notebook on triangles. \n",
"Use a copy of that code, but write it so it looks like a function.\n",
"Inside the code, you can take advantage of the fact that we\n",
"know we are working with the unit triangle, so a=(0,0),\n",
"b=(1,0) and c=(0,1), which should simplify your formulas.\n",
"\n",
"```python\n",
" def sample_unit_triangle ( m )\n",
" import numpy as np\n",
" x = np.zeros ( m )\n",
" y = np.zeros ( m )\n",
" ***\n",
" return x, y\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(0.3328231691671793, 0.15003861256765297, 0.4828617817348323)\n",
"(0.059412541338811679, 0.7415640800513531, 0.80097662139016479)\n",
"(0.31756906116274403, 0.50455454121935417, 0.82212360238209814)\n",
"(0.39398089218790633, 0.055730661989240274, 0.44971155417714659)\n",
"(0.14280066802703043, 0.30010741738312069, 0.44290808541015114)\n",
"(0.36558303235661738, 0.23346103472821803, 0.59904406708483537)\n",
"(0.2114035113264334, 0.3869664112931428, 0.59836992261957622)\n",
"(0.044025748232755471, 0.5013345484228503, 0.54536029665560581)\n",
"(0.58796551616373438, 0.13223273770710789, 0.72019825387084224)\n",
"(0.13871293504100515, 0.25635663351807514, 0.39506956855908026)\n"
]
}
],
"source": [
"# Code to sample the unit triangle\n",
"#\n",
"def sample_unit_triangle ( m ):\n",
" import numpy as np\n",
"\n",
" x = np.zeros ( m )\n",
" y = np.zeros ( m )\n",
"\n",
" for indx in range ( 0, m ):\n",
" r1 = np.random.random ( )\n",
" r2 = np.random.random ( )\n",
" r2 = np.sqrt ( r2 ) \n",
" i = 1.0 - r2\n",
" j = ( 1.0 - r1 ) * r2\n",
" k = r1 * r2\n",
" x[indx] = i * 0.0 + j * 1.0 + k * 0.0\n",
" y[indx] = i * 0.0 + j * 0.0 + k * 1.0\n",
"\n",
" return x, y\n",
"#\n",
"# Test the sampler.\n",
"#\n",
"# As long as x and y are nonnegative, and x+y is no greater than 1,\n",
"# we know (x,y) is in the unit triangle.\n",
"#\n",
"# To be sure this was a uniform sample, we would probably start\n",
"# by looking at a plot.\n",
"#\n",
"m = 10\n",
"x, y = sample_unit_triangle ( m )\n",
"for i in range ( 0, m ):\n",
" print ( x[i], y[i], x[i] + y[i] )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Density Functions for the Unit Triangle #\n",
"\n",
"Define three density functions:\n",
"\n",
"```python\n",
" def rho1(x,y):\n",
" value = 1.0\n",
" return value\n",
" \n",
" def rho2(x,y):\n",
" value = x\n",
" return value\n",
" \n",
" def rho3(x,y):\n",
" import numpy as np\n",
" value = np.exp(x+y)\n",
" return value\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Code for density functions for the unit triangle.\n",
"#\n",
"def rho1 ( x, y ):\n",
" value = 1.0\n",
" return value\n",
"\n",
"def rho2 ( x, y ):\n",
" value = x\n",
" return value\n",
"\n",
"def rho3 ( x, y ):\n",
" import numpy as np\n",
" value = np.exp ( x + y )\n",
" return value"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Center of Mass of Unit Triangle #\n",
"\n",
"Just like the centroid, we expect that the center of mass of a triangle \n",
"will always be a point inside the triangle, never outside! \n",
"\n",
"(This will also be true for any convex polgon.)\n",
"\n",
"For the uniform density rho1(x,y), we know that the center of mass\n",
"will be the same as the centroid, (1/3,1/3).\n",
"\n",
"Write a function that uses sampling to estimate the center of\n",
"mass of the unit triangle with density function rho.\n",
"If we let rho be an input quantity, we can use the same function\n",
"to check densities rho1, rho2 and rho3:\n",
"```python\n",
" def com_unit_triangle ( m, rho ):\n",
" ***\n",
" return cx, cy\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(0.32824193707783161, 0.33160961813424233)\n",
"(0.49917727466154121, 0.25318444600415191)\n",
"(0.35550779664983151, 0.35902281112936879)\n"
]
}
],
"source": [
"# Code to compute the center of mass of the unit triangle.\n",
"#\n",
"def com_unit_triangle ( m, rho ):\n",
" \n",
" import numpy as np\n",
"\n",
" mass = 0.0\n",
" cx = 0.0\n",
" cy = 0.0\n",
" x, y = sample_unit_triangle ( m )\n",
" mass = 0.0\n",
" cx = 0.0\n",
" cy = 0.0\n",
" for i in range ( 0, m ):\n",
" r = rho ( x[i], y[i] )\n",
" mass = mass + r\n",
" cx = cx + r * x[i]\n",
" cy = cy + r * y[i]\n",
" \n",
" cx = cx / mass\n",
" cy = cy / mass\n",
" return cx, cy\n",
"#\n",
"# Test each density function.\n",
"#\n",
"m = 1000\n",
"\n",
"x, y = com_unit_triangle ( m, rho1 )\n",
"print ( x, y )\n",
"\n",
"x, y = com_unit_triangle ( m, rho2 )\n",
"print ( x, y )\n",
"\n",
"x, y = com_unit_triangle ( m, rho3 )\n",
"print ( x, y )"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"# CVT for the square with a density function #\n",
"\n",
"Let's try to do a CVT calculation with a density function.\n",
"\n",
"Our region R will be the square [-1,+1]x[-1,+1].\n",
"\n",
"To estimate integrals, we still want to use uniform sampling. Because our\n",
"region has changed slightly, we will need to compute random numbers in [0,1]\n",
"and then multiply by 2 and subtract 1, so that they lie in [-1,+1]. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sample the Square\n",
"\n",
"It is easy to sample a square, but this time we have stretched\n",
"the square a little. So while the random number generator will\n",
"give us a value r between 0 and 1, we need to modify r\n",
"to get x and y values between -1 and +1. One formula for this is\n",
"\n",
"* x = 2 * r - 1\n",
"\n",
"Following the pattern we set for the unit triangle, let us wrap\n",
"up our sampling code into a function, which returns m random\n",
"values in our square. However, instead of returning x and y\n",
"as separate vectors, let us return a single Mx2 array:\n",
"```python\n",
" def sample_unit_square ( m )\n",
" import numpy as np\n",
" ***\n",
" return xy\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.73409643 -0.23928758]\n",
" [ 0.9018774 0.17383062]\n",
" [-0.06717249 -0.23507539]\n",
" [ 0.12112358 -0.07705584]\n",
" [-0.67087453 -0.78814242]]\n"
]
}
],
"source": [
"# Code to sample the square\n",
"#\n",
"def sample_square ( m ):\n",
" import numpy as np\n",
" xy = 2.0 * np.random.rand ( m, 2 ) - 1.0\n",
" return xy\n",
"#\n",
"# Test functions.\n",
"#\n",
"n = 5\n",
"g = sample_square ( n )\n",
"print ( g )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Voronoi diagram for square #\n",
"\n",
"We will want to make a Voronoi diagram from time to time\n",
"for a set of generators G. It only takes a few lines of code,\n",
"but they are hard to remember, so let us get them right and\n",
"make them a function:\n",
"```python\n",
" def display_voronoi_square ( g ):\n",
" ***compute and plot the Voronoi diagram of G\n",
" return\n",
"```\n",
"If we work at it a little, we can add a red box that\n",
"indicates where the boundary of our square is."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Code for Voronoi diagram for square\n",
"#\n",
"def display_voronoi_square ( g ):\n",
" \n",
" vor = spatial.Voronoi ( g )\n",
"#\n",
"# Because I have access to an OUTDATED copy of this code, I can't shut off\n",
"# the idiotic vertex plotting, which ruins the plot, because some of the\n",
"# vertices run off towards infinity!\n",
"#\n",
"# spatial.voronoi_plot_2d ( vor, show_vertices = False )\n",
" spatial.voronoi_plot_2d ( vor )\n",
"#\n",
"# Add a red box showing the extent of our square.\n",
"#\n",
" plt.plot ( [-1.0, 1.0, 1.0, -1.0, -1.0 ], [ -1.0, -1.0, 1.0, 1.0, -1.0 ], 'r' )\n",
" plt.axis ( 'Equal')\n",
"\n",
" return\n",
"#\n",
"# Test the function.\n",
"#\n",
"import numpy as np\n",
"g = np.array ( [ \\\n",
" [ -0.5, 0.5 ], \\\n",
" [ -0.5, -0.5 ], \\\n",
" [ 0.0, 0.0 ], \\\n",
" [ 0.5, 0.0 ] ] )\n",
"display_voronoi_square ( g )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Density for the square\n",
"\n",
"We want an interesting density thatt will vary strongly over the\n",
"square. The example we chose is\n",
"\n",
"* rho ( x, y ) = x^4 + y^4\n",
" \n",
"This is actually 0 at the origin, but that is OK. We just do not\n",
"ever want to deal with negative densities, or densities that are\n",
"zero over a substantial patch of area.\n",
"\n",
"Write a function that takes a variable holding M pairs of (x,y) values,\n",
"in an Mx2 array, and returns an M vector of density values:\n",
"```python\n",
" def density_square ( m, xy ):\n",
" ***\n",
" return value\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0 -0.5 0.5 0.125\n",
"1 -0.5 -0.5 0.125\n",
"2 0 0 0\n",
"3 0.5 0 0.0625\n"
]
}
],
"source": [
"# Code to evaluate density for the square\n",
"#\n",
"def density_square ( m, xy ):\n",
" d = np.zeros ( m )\n",
" for i in range ( 0, m ):\n",
" d[i] = xy[i,0] ** 4 + xy[i,1] ** 4\n",
" return d\n",
"#\n",
"# Test the function.\n",
"#\n",
"m = 4\n",
"xy = np.array ( [ \\\n",
" [ -0.5, 0.5 ], \\\n",
" [ -0.5, -0.5 ], \\\n",
" [ 0.0, 0.0 ], \\\n",
" [ 0.5, 0.0 ] ] )\n",
"d = density_square ( m, xy )\n",
"for i in range ( 0, m ):\n",
" print ( '%d %g %g %g' % ( i, xy[i,0], xy[i,1], d[i] ))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Centers of mass for Voronoi subregions in the square\n",
"\n",
"If we have a set of generators G in the square, we can compute a set of \n",
"M sample points in an Mx2 array S. \n",
"\n",
"Each sample point S(i) is closest to one of the generators G(j).\n",
"\n",
"We estimate the mass and the X and Y coordinates of the center of mass\n",
"of the subregion associated with G(j) by\n",
"* compute the density D of every point S(i) closest to G(j)\n",
"* to MASS[J] add the density \n",
"* to COM[J,0] add the density times the x coordinate of S(i)\n",
"* to COM[J,1] add the density times the y coordinate of S(i)\n",
"\n",
"If MASS[J] is not zero, normalize COM[J,0] and COM[J,1] by dividing by MASS[J].\n",
"Otherwise, set COM[J,0] and COM[J,1] to G[J,0] and G[J,1].\n",
"\n",
"This is your estimate for the centers of mass of each subregion."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 -0.8106 -0.8863 -0.6515 -0.9077\n",
" 1 -0.0870 0.9464 -0.1607 0.8680\n",
" 2 0.7271 -0.0192 0.8789 0.1207\n",
" 3 -0.4772 -0.0594 -0.8094 0.2644\n",
" 4 -0.0083 -0.6879 0.2203 -0.8642\n",
" 5 -0.7949 -0.6719 -0.8309 -0.5554\n",
" 6 0.3591 -0.1923 0.7414 -0.7300\n",
" 7 0.0448 -0.1959 0.0202 -0.3483\n",
" 8 0.2815 0.0432 0.4102 0.4254\n",
" 9 0.0619 -0.1233 -0.0766 0.1689\n"
]
}
],
"source": [
"# Code for centers of mass of Voronoi subregions of square with density\n",
"#\n",
"def center_of_mass_square ( g, m ):\n",
"\n",
" s = sample_square ( m )\n",
" d = density_square ( m, s )\n",
" mass = np.zeros ( n )\n",
" com = np.zeros ( [ n, 2 ] )\n",
"\n",
" for i in range ( 0, m ):\n",
" jmin = -1\n",
" dmin = np.Inf\n",
" for j in range ( 0, n ):\n",
" dj = np.linalg.norm ( s[i,:] - g[j,:] )\n",
" if ( dj < dmin ):\n",
" dmin = dj\n",
" jmin = j\n",
" mass[jmin] = mass[jmin] + d[i]\n",
" com[jmin,:] = com[jmin,:] + d[i] * s[i,:]\n",
" \n",
" for j in range ( 0, n ):\n",
" if ( mass[j] != 0.0 ):\n",
" com[j,:] = com[j,:] / mass[j]\n",
" else:\n",
" com[j,:] = g[j,:]\n",
"\n",
" return com\n",
"#\n",
"# Test function.\n",
"#\n",
"# All we can hope is that the C's are close to the G's.\n",
"#\n",
"n = 10\n",
"g = sample_square ( n )\n",
"m = 1000\n",
"c = center_of_mass_square ( g, m )\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] ) )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Compute a CVT for the square with a density #\n",
"\n",
"After all that work, the actual CVT calculation isn't too hard.\n",
"\n",
"We will\n",
"* select generators G at random using sample_square ()\n",
"* take step_num steps:\n",
"** compute the centers of mass\n",
"** replace the generators by the centers of mass\n",
"\n",
"As we go, we might try to display the Voronoi diagram each time,\n",
"except that my Voronoi plotter doesn't do a good job. So I'll\n",
"just plot the point locations."
]
},
{
"cell_type": "code",
"execution_count": 15,
"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 of square with density\n",
"#\n",
"def cvt_square_nonuniform ( n, m, step_num ):\n",
"#\n",
"# Initial G:\n",
"#\n",
" g = sample_square ( n )\n",
"# display_voronoi_square ( g )\n",
" plt.plot ( g[:,0], g[:,1], 'bo' )\n",
" plt.plot ( [-1, +1, +1, -1, -1], [-1, -1, +1, +1, -1], 'r-')\n",
" plt.axis ( 'Equal' )\n",
" plt.title ( 'Initial points')\n",
" plt.show ( )\n",
"\n",
" for step in range ( 0, step_num ):\n",
" c = center_of_mass_square ( g, m )\n",
" g = c.copy ( )\n",
" \n",
" return g\n",
"#\n",
"# Test functions.\n",
"#\n",
"n = 200\n",
"m = 10000\n",
"step_num = 10\n",
"g = cvt_square_nonuniform ( n, m, step_num )\n",
"plt.plot ( g[:,0], g[:,1], 'bo' )\n",
"plt.plot ( [-1, +1, +1, -1, -1], [-1, -1, +1, +1, -1], 'r-')\n",
"plt.axis ( 'Equal' )\n",
"plt.title ( 'Final results')\n",
"plt.show ( )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"\n"
]
}
],
"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
}