{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using the `Polygons` package for dynamic and repeated games\n",
"\n",
"In this notebook, we introduce the Polygons class and show how it can be used to compute equilibrium sets of repeated games."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Set up\n",
"\n",
"Currently the package is not available through the Julia package distribution system. To download and use the package (and the convex hull routines one which it depends), remove the comment markers (#) and run the following in an interactive Julia session:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Pkg.clone(\"https://github.com/squipbar/Polygons.jl\")\n",
"#Pkg.clone(\"https://github.com/cc7768/CHull2d.jl\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that you do not need to do this every time. Once you have it, re-downloading will yield errors.\n",
"\n",
"\n",
"## Defining polygon objects\n",
"\n",
"The `Polygon` class stores a polygon object is two ways. First, as a collection of vertices, `pts`. This is not a terribly helpful representation, mathematically speaking, and is useful mostly for visualizing the object. The second representation is the normal-distance representation of the sides of the polygon. The is a mathematially useful object. If $G$ is an $n\\times2$ matrix of vector normals $g_i$ of the polygon faces, and $m$ is a length-$n$ vector of distances, then the point $z\\in\\mathbb R^2$ is inside the polygon if:\n",
"\n",
"$$ g_i z \\le m_i \\ \\forall \\ i = 1, \\ldots, n$$\n",
"\n",
"Ths representation is stored inside the `Polygon` class via the terms `dirs` (for $G$) and `dists` (for $m$). Except where explicitly noted below, polygons are assumed to be convex.\n",
"\n",
"Let's define our first polygons:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using Polygons"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Z = [ 1 -1; 1 1; -1 1; -1 -1 ]\n",
"G = [ 1 0; 0 1; -1 0; 0 -1 ]\n",
"m = [ 1, 1, 1, 1 ]\n",
"a = Polygon( Z, G, m ) ;\n",
"b = Polygon( pts=Z ) ;\n",
"c = Polygon( dirs=G, dists=m ) ;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because the class constructors standardize the polygon representation (points and normals are ordered anti-clockwise starting at lowest point, normals are unit-length), then we can throw badly-formed representations at the constructors and still preserve equality."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Polygons.Polygon(4x2 Array{Float64,2}:\n",
" 1.0 -1.0\n",
" 1.0 1.0\n",
" -1.0 1.0\n",
" -1.0 -1.0,4x2 Array{Float64,2}:\n",
" 1.0 0.0\n",
" 0.0 1.0\n",
" -1.0 0.0\n",
" 0.0 -1.0,[1.0,1.0,1.0,1.0])"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"badZ = Z[ [1, 3, 4, 2], : ]\n",
"d = Polygon(pts=badZ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So `d` and `a` are the same thing (sadly the `==` operator doesn't work as well as hoped here). But we can eyeball the two to check that everything is ok:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Polygons.Polygon(4x2 Array{Float64,2}:\n",
" 1.0 -1.0\n",
" 1.0 1.0\n",
" -1.0 1.0\n",
" -1.0 -1.0,4x2 Array{Float64,2}:\n",
" 1.0 0.0\n",
" 0.0 1.0\n",
" -1.0 0.0\n",
" 0.0 -1.0,[1.0,1.0,1.0,1.0])"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting a polygon is easy. Just use the `polyPlot` function:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"polyPlot(a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Look - it's a square. Awesome!\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"## Vector addition, scalar multiplication\n",
"\n",
"We can add a vector to a polygon as follows"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f = a + [ 1, 1 ] ;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Scalar multiplication will expand the set along rays through the origin"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"g = 1.5 * a ;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Becuase `polyPlot` also accepts polygon arrays, we can easily compare multiple `polygon` objects in the same plot"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"polyPlot( [a, f, g] )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Approximate set sums\n",
"\n",
"### Theory\n",
"\n",
"The `polygon` package also allows us to compute approximate set sums of polygons. Denote by $Z_1$ and $Z_2$ the vertices of two polygons $P_1$ and $P_2$. If $H$ is a $m \\times 2$ matrix of search directions, the extremal points of the two polygons in the search direction $h_i$ are given by the maximal projection of the vertices onto $h_i$:\n",
"\n",
"$$ m^1_i = max Z_1 h_i' \\qquad\\qquad\\qquad m^2_i = max Z_2 h_i' $$\n",
"\n",
"Let $m^j=(m_i^j)$ be the vector of these projections and $m = m^1 + m^2$. Then we can form a Judd-Yeltekin-Conklin outer aproximation to the set sum as the set with faces with normals given by the rows of $H$ and distances $m$:\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" P &= \\{ v_1 + v_2 | v_1 \\in P_1, \\ v_2 \\in P_2 \\} \\\\\n",
" &\\subseteq \\{ z | h_i.z \\le m_i \\ \\forall i \\}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"We can also connect the vertices which maximize the projected distance vectors to form an inner approximation to the setsum. Let $\\hat z_i^j$ be the vertex of polygon $j$ which achieves the maximum distance in the direction $h_i$. Then if $\\hat P$ is the polygon defined by the vertices $\\hat z_i = z_i^1 + z_i^2$:\n",
"\n",
"$$ \\hat P \\subseteq P$$\n",
"\n",
"Note that we could also just take the convex hull of all the possible point sums of the vertices $Z_1$ and $Z_2$ and recover a poterntially larger inner approximation. But this is, in general, a fairly slow algorithm. Projecting onto thevectors in $H$ is faster, as it involves less sorting.\n",
"\n",
"### Practice\n",
"\n",
"So, let's try these methods:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h = Polygon( pts=[ 0 1 ; 1 0 ; -1 0 ; 0 -1 ] )\n",
"exactdirs = [ a.dirs ; h.dirs ]\n",
"# The union of the face directions. Will give an exact result\n",
"exact = setSum( a, h, exactdirs )\n",
"# The exact set sum\n",
"apxdirs = [ .1 1; .2 .8 ; -1 1 ; -.7 -.3 ; .4 -.6; .6 -.4 ]\n",
"outer = setSum( a, h, apxdirs ) ;\n",
"inner = setSum( a, h, apxdirs, false ) ;\n",
"polyPlot( [ a, h, exact, inner, outer ])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The polygons `a` and `h` are in red and blue in the centre of the diagram. The exact set sum is in black, and is achieved when the search directions are the same as the union of nomals of the two sets. The cyan and green polygons defin the inner and outer approximations to this. Here, they are quite inaccuracte as there are few search directions and the differ in general from the true normals. If we increase the number of search directions then we will get much better approximations."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 40\n",
"theta = 2 * pi * rand( N )\n",
"manydirs = zeros( N, 2 )\n",
"for i in 1:N\n",
" manydirs[i,:] = [ cos(theta[i]), sin(theta[i]) ]\n",
"end\n",
"manyouter = setSum( a, h, manydirs ) ;\n",
"manyinner = setSum( a, h, manydirs, false ) ;\n",
"polyPlot( [ exact, manyinner, manyouter ])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, the inner approximation *is* the exact soluttion, as it picks out all the vertices perfectly\n",
"\n",
"Not discussed here, but useful functions to know about: `setsum` can apply to arrays, `weightedSum` can perform a weighted set sum, and `deeDoop` removes duplicate and near-duplicate points from a polygon.\n",
"\n",
"## Convex hull\n",
"\n",
"The `polygon` package also provides a routine to generate a convex hull using the Graham Scan method"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Gadfly\n",
"using Colors\n",
"N = 40\n",
"srand(42)\n",
"manypts = randn( N, 2 )\n",
"chull( manypts )\n",
"ch = Polygon( pts = chull( manypts ) )\n",
"chplot = ch.pts[ [ 1:end; 1], : ]\n",
"plot( layer( x=manypts[:,1], y=manypts[:,2], Geom.point ), layer( x=chplot[:,1], y=chplot[:,2], Geom.path,\n",
" Theme(default_color=colorant\"red\") ) ) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cropping polygons\n",
"\n",
"Cropping polygons is a fundamental operation when calculating incentive-compatilbility with a fixed deviating payoff (when the deviating payoff is not fixed, we need to use linear programming or other methods to find the deviating payoff). The `crop` function does this - chopping the polyon in either the x or y direction as required."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"