Using the Polygons package for dynamic and repeated games

In this notebook, we introduce the Polygons class and show how it can be used to compute equilibrium sets of repeated games.

Set up

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:

In [1]:
#Pkg.clone("https://github.com/squipbar/Polygons.jl")
#Pkg.clone("https://github.com/cc7768/CHull2d.jl")

Note that you do not need to do this every time. Once you have it, re-downloading will yield errors.

Defining polygon objects

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:

$$ g_i z \le m_i \ \forall \ i = 1, \ldots, 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.

Let's define our first polygons:

In [2]:
using Polygons
In [3]:
Z = [ 1 -1; 1 1; -1 1; -1 -1 ]
G = [ 1 0; 0 1; -1 0; 0 -1 ]
m = [ 1, 1, 1, 1 ]
a = Polygon( Z, G, m ) ;
b = Polygon( pts=Z ) ;
c = Polygon( dirs=G, dists=m ) ;

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.

In [4]:
badZ = Z[ [1, 3, 4, 2], : ]
d = Polygon(pts=badZ)
Out[4]:
Polygons.Polygon(4x2 Array{Float64,2}:
  1.0  -1.0
  1.0   1.0
 -1.0   1.0
 -1.0  -1.0,4x2 Array{Float64,2}:
  1.0   0.0
  0.0   1.0
 -1.0   0.0
  0.0  -1.0,[1.0,1.0,1.0,1.0])

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:

In [5]:
a
Out[5]:
Polygons.Polygon(4x2 Array{Float64,2}:
  1.0  -1.0
  1.0   1.0
 -1.0   1.0
 -1.0  -1.0,4x2 Array{Float64,2}:
  1.0   0.0
  0.0   1.0
 -1.0   0.0
  0.0  -1.0,[1.0,1.0,1.0,1.0])

Plotting a polygon is easy. Just use the polyPlot function:

In [6]:
polyPlot(a)
Out[6]:
x -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 y

Look - it's a square. Awesome!

Vector addition, scalar multiplication

We can add a vector to a polygon as follows

In [7]:
f = a + [ 1, 1 ] ;

Scalar multiplication will expand the set along rays through the origin

In [8]:
g = 1.5 * a ;

Becuase polyPlot also accepts polygon arrays, we can easily compare multiple polygon objects in the same plot

In [9]:
polyPlot( [a, f, g] )
Out[9]:
x -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 y

Approximate set sums

Theory

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$:

$$ m^1_i = max Z_1 h_i' \qquad\qquad\qquad m^2_i = max Z_2 h_i' $$

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$:

$$ \begin{aligned} P &= \{ v_1 + v_2 | v_1 \in P_1, \ v_2 \in P_2 \} \\ &\subseteq \{ z | h_i.z \le m_i \ \forall i \} \end{aligned} $$

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$:

$$ \hat P \subseteq P$$

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.

Practice

So, let's try these methods:

In [10]:
h = Polygon( pts=[ 0 1 ; 1 0 ; -1 0 ; 0 -1 ] )
exactdirs = [ a.dirs ; h.dirs ]
# The union of the face directions.  Will give an exact result
exact = setSum( a, h, exactdirs )
# The exact set sum
apxdirs = [ .1 1; .2 .8 ; -1 1 ; -.7 -.3 ; .4 -.6; .6 -.4  ]
outer = setSum( a, h, apxdirs ) ;
inner = setSum( a, h, apxdirs, false ) ;
polyPlot( [ a, h, exact, inner, outer ])
Out[10]:
x -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 -20 -10 0 10 20 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 -20 -10 0 10 20 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 y

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.

In [11]:
N = 40
theta = 2 * pi * rand( N )
manydirs = zeros( N, 2 )
for i in 1:N
    manydirs[i,:] = [ cos(theta[i]), sin(theta[i]) ]
end
manyouter = setSum( a, h, manydirs ) ;
manyinner = setSum( a, h, manydirs, false ) ;
polyPlot( [ exact, manyinner, manyouter ])
Out[11]:
x -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 -9.0 -8.8 -8.6 -8.4 -8.2 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 -10 -5 0 5 10 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 -9.0 -8.8 -8.6 -8.4 -8.2 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 -10 -5 0 5 10 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 y

In this case, the inner approximation is the exact soluttion, as it picks out all the vertices perfectly

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.

Convex hull

The polygon package also provides a routine to generate a convex hull using the Graham Scan method

In [26]:
using Gadfly
using Colors
N = 40
srand(42)
manypts = randn( N, 2 )
chull( manypts )
ch = Polygon( pts = chull( manypts ) )
chplot = ch.pts[ [ 1:end; 1], : ]
plot( layer( x=manypts[:,1], y=manypts[:,2], Geom.point ), layer( x=chplot[:,1], y=chplot[:,2],  Geom.path,
          Theme(default_color=colorant"red") ) )
Out[26]:
x -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -20 -10 0 10 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 y

Cropping polygons

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.

In [27]:
chr = crop( ch, 1, .5, true )
chdown = crop( ch, 2, .5, false )
polyPlot( [chr, chdown, ch])
Out[27]:
x -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -20 -10 0 10 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 y

Unions

The union function forms the convex union of two polygons (or an array of polygons).

In [28]:
polyPlot( [ exact, ch, union( [ exact, ch ] ) ] )
Out[28]:
x -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -20 -10 0 10 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 -6 -3 0 3 6 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 y