{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating and solving the game"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iteration1\n",
"Distance0.96599025767\n",
"Iteration2\n",
"Distance0.724492693252\n",
"Iteration5\n",
"Distance0.308667670965\n",
"Iteration10\n",
"Distance0.0803995842487\n",
"Iteration25\n",
"Distance0.00107441670934\n",
"34\n",
"TOC: Elapsed: 0:00:47.84\n"
]
}
],
"source": [
"import numpy as np\n",
"from quantecon.game_theory import *\n",
"from scipy.spatial import HalfspaceIntersection\n",
"\n",
"pd_payoff = np.array([[9.0, 1.0], [10.0, 3.0]])\n",
"\n",
"A = Player(pd_payoff)\n",
"B = Player(pd_payoff)\n",
"nfg = NormalFormGame((A, B))\n",
"\n",
"rpd = RepeatedGame(nfg, 0.75)\n",
"\n",
"result_32, dist_history, time_history, figure_history, (C, H, Z) = \\\n",
"rpd.outerapproximation(nH=32, tol=1e-4, linprog_method='simplex', verbose=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Algorithm animation"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# With 32 gradients\n",
"from matplotlib import pyplot as plt, animation, rcParams, rc\n",
"\n",
"rc('animation',html='html5')\n",
"rcParams['figure.figsize'] = 3, 3\n",
"\n",
"x = figure_history[0][:,0]\n",
"y = figure_history[0][:,1]\n",
"\n",
"def animate(i):\n",
" x = figure_history[i][:,0]\n",
" y = figure_history[i][:,1]\n",
" line.set_data(x,y)\n",
" return (line,)\n",
"\n",
"fig = plt.figure(dpi=80,figsize=(7,7))\n",
"ax = fig.add_subplot(111)\n",
"ax.set_xlim((-2,12.5))\n",
"ax.set_ylim((-2,12.5))\n",
"plt.scatter(x,y)\n",
"line, = ax.plot([],[],'o',lw=2)\n",
"animation.FuncAnimation(fig,animate,len(figure_history), interval=300)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"34\n",
"TOC: Elapsed: 0:00:18.02\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# With 16 gradients\n",
"result_16, dist_history, time_history, figure_history, (C, H, Z) = \\\n",
"rpd.outerapproximation(nH=16, tol=1e-4, linprog_method='simplex')\n",
"\n",
"x = figure_history[0][:,0]\n",
"y = figure_history[0][:,1]\n",
"\n",
"fig = plt.figure(dpi=80,figsize=(7,7))\n",
"ax = fig.add_subplot(111)\n",
"ax.set_xlim((-10,20))\n",
"ax.set_ylim((-10,20))\n",
"plt.scatter(x,y)\n",
"line, = ax.plot([],[],'o',lw=2)\n",
"animation.FuncAnimation(fig, animate, len(figure_history), interval=300)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"34\n",
"TOC: Elapsed: 0:01:14.02\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# With 31 gradients\n",
"result_31, dist_history, time_history, figure_history, (C, H, Z) = \\\n",
"rpd.outerapproximation(nH=31, tol=1e-4, linprog_method='simplex')\n",
"\n",
"x = figure_history[0][:,0]\n",
"y = figure_history[0][:,1]\n",
"\n",
"fig = plt.figure(dpi=80,figsize=(7,7))\n",
"ax = fig.add_subplot(111)\n",
"ax.set_xlim((-10,20))\n",
"ax.set_ylim((-10,20))\n",
"plt.scatter(x,y)\n",
"line, = ax.plot([],[],'o',lw=2)\n",
"animation.FuncAnimation(fig, animate, len(figure_history), interval=300)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"34\n",
"TOC: Elapsed: 0:22:37.88\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# With 128 gradients\n",
"result_128, dist_history, time_history, figure_history, (C, H, Z) = \\\n",
"rpd.outerapproximation(nH=128, tol=1e-4, linprog_method='simplex')\n",
"\n",
"x = figure_history[0][:,0]\n",
"y = figure_history[0][:,1]\n",
"\n",
"fig = plt.figure(dpi=80,figsize=(7,7))\n",
"ax = fig.add_subplot(111)\n",
"ax.set_xlim((-10,20))\n",
"ax.set_ylim((-10,20))\n",
"plt.scatter(x,y)\n",
"line, = ax.plot([],[],'o',lw=2)\n",
"animation.FuncAnimation(fig, animate, len(figure_history), interval=300)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"35\n",
"TOC: Elapsed: 0:18:26.39\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# With 127 gradients\n",
"result_127, dist_history, time_history, figure_history, (C, H, Z) = \\\n",
"rpd.outerapproximation(nH=127, tol=1e-4, linprog_method='simplex')\n",
"\n",
"x = figure_history[0][:,0]\n",
"y = figure_history[0][:,1]\n",
"\n",
"fig = plt.figure(dpi=80,figsize=(7,7))\n",
"ax = fig.add_subplot(111)\n",
"ax.set_xlim((-10,20))\n",
"ax.set_ylim((-10,20))\n",
"plt.scatter(x,y)\n",
"line, = ax.plot([],[],'o',lw=2)\n",
"animation.FuncAnimation(fig, animate, len(figure_history), interval=300)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAHMCAYAAADieBDmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAMTQAADE0B0s6tTgAAFfZJREFUeJzt3W+IXfd95/HP1zNDvKkgtbxNPDvjNAp2Cls3Sb0RTfZP\n0k1RUSjbB9GSEOwu5FHCNhDZBEoDCwvL0mygRpg8aEIKCakKIRF1oE1VVDZ/+s9E4CYbsVAnWK49\nk0m6sR0aETuZGf32gUZEkvVvZq71nbl6veAymnPvPfpyfKz3nHvvOVNjjAAAN94t3QMAwM1KhAGg\niQgDQBMRBoAmIgwATUQYAJqIMAA02VSEq+rWqnqkqh6vqm9U1Ymqumvjvi9X1emq+vrG7YGXZmQA\nmA6zW3jOJ5L8+RhjVNUHknwyya9u3PfAGOORSQ0HANNsU0fCY4wXxhhfHD+9zNajSV4z8akA4CZQ\n27lsZVV9JsmzY4wPVtWXk/yrJD9O8n+T/O4Y44nLPOfBJA+e/35mZmbhjjvu2PIMANBpeXn5J2OM\nl23luVuOcFV9OMl/SvJrY4wfVdWdY4ynq6qS/HaS/zrG+NfXWs/i4uJYWlra0gwA0K2qlscYi1t5\n7pY+HV1VH0ryziTvGGP8KEnGGE9vfB1jjI8leW1V3b6V9QPAzWDTEd54Ofk9SQ6MMX6wsWy2ql51\nwWMOJfneGOOZiU0KAFNmU5+OrqrFJL+f5IkkXzr3ynN+nOTtSf6sql6W5GyS7yf5zcmOCgDTZVMR\nHmMsJakr3P2m7Y8DADcPV8wCgCYiDABNRBgAmogwADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABN\nRBgAmogwADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmogwADQRYQBoIsIA0ESEAaCJCANA\nExEGgCYiDABNRBgAmogwADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmogwADQRYQBoIsIA\n0ESEAaCJCANAExEGgCYiDABNRBgAmogwADQRYQBoIsIA0GS2ewB2l7Pr6zl2/EROL69k38J8Dh08\nkFtmZrrHAtiVaozROsDi4uJYWlpqnYHr89ipUzl89GRWxm2ZzXrWMpP5ei5H7tufe++5p3s8gBZV\ntTzGWNzKc70czXU5u76ew0dPZnnszWrm8nxuzWrmsjz25vDRr+Xs+nr3iAC7jghzXY4dP5GVcVvW\nL3kHYz2zWRl7c+z4iabJAHYvEea6nF5eyWwuf7Q7l/WcXl65wRMB7H4izHXZtzCftVz+A1irmcm+\nhfkbPBHA7ifCXJdDBw9kvp7LTNYuWj6TtczXszl08EDTZAC7lwhzXW6ZmcmR+/ZnoZ7NXFbz8ryQ\nuaxmoZ7Jw/e/2WlKAFvgFCU2xXnCABfbzilKIgwA2+A8YQDYhUQYAJqIMAA0EWEAaCLCANBkUxGu\nqlur6pGqeryqvlFVJ6rqro37XllVx6vqW1V1qqre+tKMDADTYStHwp9I8gtjjDck+UKST24s/0iS\nR8cYdyd5b5I/rqq5yYwJANNnUxEeY7wwxvji+OnJxY8mec3Gn9+V5A82HncyyXeSvG1CcwLA1Nnu\ne8IfTPKFqro9ydwY47sX3Pdkkldvc/0AMLVmr/2Qy6uqDye5K8mvJfkXm3jeg0kePP/9K17xiq2O\nAAC72paOhKvqQ0nemeQdY4wfjTGeSbJWVXdc8LDXJHnq0ueOMR4aYyyev+3Zs2crIwDArrfpCG8c\nyb4nyYExxg8uuOtzSd6/8Zj9SRaSfGUSQwLANNrUy9FVtZjk95M8keRLVZUkPx5j/EqS30nymar6\nVpKfJLl/jLE64XkBYGpsKsJjjKUkdYX7vpfk1ycxFADcDFwxCwCaiDAANBFhAGgiwgDQRIQBoIkI\nA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgy\n2z0AANd2dn09x46fyOnllexbmM+hgwdyy8xM91hsU40xWgdYXFwcS0tLrTMA7GSPnTqVw0dPZmXc\nltmsZy0zma/ncuS+/bn3nnu6x7vpVdXyGGNxK8/1cjTADnZ2fT2Hj57M8tib1czl+dya1cxleezN\n4aNfy9n19e4R2QYRBtjBjh0/kZVxW9YvefdwPbNZGXtz7PiJpsmYBBEG2MFOL69kNpc/2p3Lek4v\nr9zgiZgkEQbYwfYtzGctl/8A1mpmsm9h/gZPxCSJMMAOdujggczXc5nJ2kXLZ7KW+Xo2hw4eaJqM\nSRBhgB3slpmZHLlvfxbq2cxlNS/PC5nLahbqmTx8/5udprTLOUUJYBdwnvDOtZ1TlEQYALbBecIA\nsAu5bOWEeckIgOvl5egJcmk5gJuPl6N3AJeWA2CzRHhCXFoOgM0S4QlxaTkANkuEJ8Sl5QDYLBGe\nEJeWA2CzRHhCXFoOgM1yitKEOU8Y4ObispUA0MR5wgCwC4kwADQRYQBoIsIA0ESEAaCJCANAExEG\ngCaz134IADcjFx966blYBwAv8tipUzl89GRWxm2ZzXrWMpP5ei5H7tufe++5p3u8HcXFOgCYmLPr\n6zl89GSWx96sZi7P59asZi7LY28OH/1azq5f/te2snkiDMBFjh0/kZVxW9YvecdyPbNZGXtz7PiJ\npsmmjwgDcJHTyyuZzeWPdueyntPLKzd4ouklwgBcZN/CfNZy+Q9grWYm+xbmb/BE00uEAbjIoYMH\nMl/PZSZrFy2fyVrm69kcOnigabLpI8IAXOSWmZkcuW9/FurZzGU1L88LmctqFuqZPHz/m52mNEFO\nUQLgspwnfH22c4qSCAPANjhPGAB2IREGgCYiDABNRBgAmogwADTZVISr6uGqerKqRlW98YLlX66q\n01X19Y3bA5MfFQCmy2Z/n/Dnk3w0yV9f5r4HxhiPbH8kALg5bCrCY4yvJklVvTTTAMBNZJLvCX+0\nqr5ZVZ+tqtde6UFV9WBVLZ2/nTlzZoIjAMDuMakI/9YY43VJXp/kr5L86ZUeOMZ4aIyxeP62Z8+e\nCY0AALvLRCI8xnh64+sYY3wsyWur6vZJrBsAptW2I1xVs1X1qgu+P5Tke2OMZ7a7bgCYZpv6YFZV\nfTzJbyS5I8lfVNUPk7whyZ9V1cuSnE3y/SS/OelBAWDabPbT0e+7wl1vmsAsAHBTccUsAGgiwgDQ\nRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAA\nNBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIM\nAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkI\nA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgi\nwgDQRIQBoIkIA0ATEQaAJpuKcFU9XFVPVtWoqjdesPyVVXW8qr5VVaeq6q2THxUApstmj4Q/n+Tf\nJ/nHS5Z/JMmjY4y7k7w3yR9X1dwE5gOAqTW7mQePMb6aJFV16V3vSnLXxmNOVtV3krwtyV9OYEYA\nmErbfk+4qm5PMjfG+O4Fi59M8urtrhsAptkN/2BWVT1YVUvnb2fOnLnRIwDAjrDtCI8xnkmyVlV3\nXLD4NUmeusLjHxpjLJ6/7dmzZ7sjAMCuNKkj4c8leX+SVNX+JAtJvjKhdQPAVNrsKUofr6qlJItJ\n/qKqvr1x1+8k+bdV9a0kn0py/xhjdaKTAsCU2eyno993heXfS/LrE5kIAG4SrpgFAE1EGACaiDAA\nNBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIM\nAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkI\nA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgi\nwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACa\niDAANBFhAGgiwgDQRIQBoMnsJFdWVU8m+XGS5zcW/d4Y47OT/DsAYFpMNMIb3j3G+PpLsF4AmCpe\njgaAJi9FhD9TVd+sqj+sqp+79M6qerCqls7fzpw58xKMAAA736Qj/NYxxi8luTfJ95N8+tIHjDEe\nGmMsnr/t2bNnwiMAwO4w0feExxhPbXxdraojSR6f5PoBYJpM7Ei4qn6mqn72gkXvSfL3k1o/AEyb\nSR4JvyrJsaqaSVJJnkjyXya4fgCYKhOL8BjjiSS/PKn1AcC0c4oSADQRYQBoIsIA0ESEAaCJCANA\nExEGgCYiDABNRBgAmogwADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmogwADQRYQBoIsIA\n0ESEAaCJCANAExEGgCYiDABNRBgAmogwADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmogw\nADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmsx2DwDAznR2fT3Hjp/I6eWV7FuYz6GDB3LL\nzEz3WFOlxhitAywuLo6lpaXWGQC42GOnTuXw0ZNZGbdlNutZy0zm67kcuW9/7r3nnu7xdpSqWh5j\nLG7luV6OBuAiZ9fXc/joySyPvVnNXJ7PrVnNXJbH3hw++rWcXV/vHnFqiDAAFzl2/ERWxm1Zv+Qd\ny/XMZmXszbHjJ5ommz4iDMBFTi+vZDaXP9qdy3pOL6/c4ImmlwgDcJF9C/NZy+U/gLWamexbmL/B\nE00vEQbgIocOHsh8PZeZrF20fCZrma9nc+jggabJpo8IA3CRW2ZmcuS+/VmoZzOX1bw8L2Quq1mo\nZ/Lw/W92mtIEOUUJgMtynvD12c4pSiIMANvgPGEA2IVEGACaiDAANBFhAGgiwgDQRIQBoIkIA0AT\nEQaAJiIMAE1EGACaiDAANJlohKvq7qr626p6vKpOVtUvTnL9ADBNJn0k/PEknxhjvC7J/0ryqQmv\nHwCmxsQiXFWvTPKmJH+0sehYkjur6q5J/R0AME0meSR8Z5KVMcZakoxzvyPxqSSvnuDfAQBT44Z/\nMKuqHqyqpfO3M2fO3OgRAGBHmGSEn04yX1WzSVJVlXNHwU9d+KAxxkNjjMXztz179kxwBADYPSYW\n4THGPyV5LMn9G4sOJVkaY3x7Un8HAEyT2Qmv731JPlVVH07yz0neO+H1A8DUmGiExxj/kOQtk1wn\nAEwrV8wCgCYiDABNRBgAmogwADQRYQBoIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmogwADQRYQBo\nIsIA0ESEAaCJCANAExEGgCYiDABNRBgAmsx2D8DucXZ9PceOn8jp5ZXsW5jPoYMHcsvMTPdYALtW\njTFaB1hcXBxLS0utM3Btj506lcNHT2Zl3JbZrGctM5mv53Lkvv259557uscDaFNVy2OMxa0818vR\nXNPZ9fUcPnoyy2NvVjOX53NrVjOX5bE3h49+LWfX17tHBNiVRJhrOnb8RFbGbVm/5N2L9cxmZezN\nseMnmiYD2N1EmGs6vbyS2Vz+aHcu6zm9vHKDJwKYDiLMNe1bmM9aLv8BrNXMZN/C/A2eCGA6iDDX\ndOjggczXc5nJ2kXLZ7KW+Xo2hw4eaJoMYHcTYa7plpmZHLlvfxbq2cxlNS/PC5nLahbqmTx8/5ud\npgSwRU5R4ro5TxjgxbZzipIIA8A2OE8YAHYhEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkI\nA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgi\nwgDQRIQBoIkIA0ATEQaAJiIMAE1EGACaiDAANBFhAGgiwgDQRIQBoIkIA0CTiUS4qv57Vf2/qvr6\nxu3oJNYLANNsdoLrOjrGODzB9QHAVPNyNAA0meSR8Luq6u1Jvp/kf4wxvnS5B1XVg0kevGDR2apa\nmeAc02hPkjPdQ+xgts/V2T7XZhtdne1zdXds9Yk1xrj2g6r+LsndV7j7l5OsJnlmjLFaVf8uyZ8k\n2T/G+MfrWPfSGGNxEzPfdGyjq7N9rs72uTbb6Opsn6vbzva5riPhMcZbrneFY4y/qaq/T/KmJNeM\nMADcrCb16ejFC/58d5I3JvnmJNYNANNqUu8J/8+q+jdJ1pKsJ/ntMcbj1/nchyY0wzSzja7O9rk6\n2+fabKOrs32ubsvb57reEwYAJs8pSgDQRIQBoIkIA0CTHRFh156+vKq6u6r+tqoer6qTVfWL3TPt\nJFX1ZFX9wwX7zbu7Z+pUVQ9vbJNRVW+8YPkrq+p4VX2rqk5V1Vs75+x0lW305ao6fcG+9EDnnF2q\n6taqemTj35xvVNWJqrpr476bfj+6xvbZ0j40yStmbZdrT7/Yx5N8Yozxqar6z0k+lWR/70g7zrvH\nGF/vHmKH+HySjyb560uWfyTJo2OMg1W1P8mfVNW+McbqDZ+w35W2UZI8MMZ45AbPsxN9IsmfjzFG\nVX0gySeT/GrsR+ddafskW9iHdsSRMC9WVa/MuQue/NHGomNJ7jz/Uxdcaozx1THG0mXueleSP9h4\nzMkk30nyths5205xlW1EkjHGC2OML46fnjbzaJLXbPz5pt+PrrF9tmQnRfhdVfV/qup/V9V/7B5m\nB7gzycoYYy1JNv6jP5Xk1a1T7TyfqapvVtUfVtXPdQ+z01TV7UnmxhjfvWDxk7EfXc5HN/alz1bV\na7uH2SE+mOQL9qMr+mCSL1zw/ab3oRsS4ar6u6r6/hVud+bcT1c/P8Z4fZL/luSzVfXzN2I2drW3\njjF+Kcm9OfeLQz7dPA+712+NMV6X5PVJ/irJnzbP066qPpzkriS/2z3LTnSZ7bOlfeiGRHiM8ZYx\nxr+8wu3pMcZ3z7+vMMb4myTnrz19M3s6yXxVzSZJVVXO/dT5VOtUO8gY46mNr6tJjiT5D70T7Txj\njGeSrFXVhb/l5TWxH11kjPH0xtcxxvhYktduHP3dlKrqQ0nemeQdY4wf2Y8udun2Sba+D+2Il6Nd\ne/rFxhj/lOSxJPdvLDqUZGmM8e2+qXaOqvqZqvrZCxa9J+d+eOPFPpfk/Umy8YGahSRfaZ1oB6mq\n2ap61QXfH0ryvY3w3HQ2ft3se5IcGGP84IK77Ee5/PbZzj60Iy5bWVWfTnLhtad/b4zx+d6p+lXV\nL+TcJ6JvT/LPSd47xripfzg5b+P9lmNJZpJUkieSfHCM8WTnXJ2q6uNJfiPnfrfpM0l+OMa4a+Mf\nh88k2ZfkJ0k+cKXf9z3tLreNkrwh52LysiRnc+6tjQfHGN/omrPLxgHR0zn3/9MPNxb/eIzxK/aj\nK2+fJG/PFvehHRFhALgZ7YiXowHgZiTCANBEhAGgiQgDQBMRBoAmIgwATUQYAJqIMAA0+f911WsM\nq6OG4gAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# A moderately good guess of the value set\n",
"guess_C = np.array([ 20.00017115, 25.45596399, 21.16595899,\n",
" 13.08174746, -2.08048976, 2.43087133])\n",
"\n",
"guess_Z = np.array([[20.00012836, 18.00009409, 12.34435949,\n",
" 1.24983985, 1.4998567, 20.00012836],\n",
" [1.24983985, 18.00003305, 19.12512836,\n",
" 20.00012836, 1.50015027, 1.24983985]])\n",
"\n",
"guess_H = np.array([[ 1.00000000e+00, 0.00000000e+00],\n",
" [ 7.07106781e-01, 7.07106781e-01],\n",
" [ 1.95090322e-01, 9.80785280e-01],\n",
" [ -7.07106781e-01, 7.07106781e-01],\n",
" [ -9.80785280e-01, -1.95090322e-01],\n",
" [ 1.95090322e-01, -9.80785280e-01]])\n",
"\n",
"guess_vertices = \\\n",
"HalfspaceIntersection(np.column_stack((guess_H, -guess_C)), np.mean(guess_Z, axis=1)).intersections\n",
"\n",
"x = guess_vertices[:][:, 0]\n",
"y = guess_vertices[:][:, 1]\n",
"plt.scatter(x,y)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"38\n",
"TOC: Elapsed: 0:00:2.51\n"
]
},
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"result_guess, dist_history, time_history, figure_history, (C, H, Z) = \\\n",
"rpd.outerapproximation(nH=6, tol=1e-4, linprog_method='simplex', init_guess=(guess_C, guess_H, guess_Z))\n",
"\n",
"x = figure_history[0][:,0]\n",
"y = figure_history[0][:,1]\n",
"\n",
"fig = plt.figure(dpi=80,figsize=(7,7))\n",
"ax = fig.add_subplot(111)\n",
"ax.set_xlim((-5,25))\n",
"ax.set_ylim((-5,25))\n",
"plt.scatter(x,y)\n",
"line, = ax.plot([],[],'o',lw=2)\n",
"animation.FuncAnimation(fig, animate, len(figure_history), interval=300)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Area comparison"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Number of gradients | \n",
" Area | \n",
" Area relative to min area | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 16 | \n",
" 26.3979 | \n",
" 102.6517 | \n",
"
\n",
" \n",
" 1 | \n",
" 31 | \n",
" 27.6691 | \n",
" 107.5949 | \n",
"
\n",
" \n",
" 2 | \n",
" 32 | \n",
" 26.197 | \n",
" 101.8704 | \n",
"
\n",
" \n",
" 3 | \n",
" 127 | \n",
" 26.1029 | \n",
" 101.5045 | \n",
"
\n",
" \n",
" 4 | \n",
" 128 | \n",
" 25.716 | \n",
" 100.0 | \n",
"
\n",
" \n",
" 5 | \n",
" 6 (not from unit circle) | \n",
" 34.282 | \n",
" 133.31 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Number of gradients Area Area relative to min area\n",
"0 16 26.3979 102.6517\n",
"1 31 27.6691 107.5949\n",
"2 32 26.197 101.8704\n",
"3 127 26.1029 101.5045\n",
"4 128 25.716 100.0\n",
"5 6 (not from unit circle) 34.282 133.31"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import pandas as pd\n",
"from scipy.spatial import ConvexHull\n",
"\n",
"columns = [\"Number of gradients\", \"Area\", \"Area relative to min area\"]\n",
"nb_gradients = np.array([\"16\", \"31\", \"32\", \"127\", \"128\", \"6 (not from unit circle)\"])\n",
"results = np.array([result_16, result_31, result_32, result_127, result_128, result_guess])\n",
"convex_hulls = np.array([ConvexHull(x) for x in results])\n",
"areas = np.array([round(x.area, 4) for x in convex_hulls])\n",
"relative_areas = np.array([round(area/min(areas)*100, 4) for area in areas])\n",
"\n",
"df = pd.DataFrame(data=np.transpose((nb_gradients, areas, relative_areas)), columns=columns)\n",
"df"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"34\n",
"TOC: Elapsed: 0:00:14.17\n",
"34\n",
"TOC: Elapsed: 0:00:43.36\n",
"34\n",
"TOC: Elapsed: 0:00:45.22\n",
"34\n",
"TOC: Elapsed: 0:01:12.22\n",
"36\n",
"TOC: Elapsed: 0:01:56.73\n",
"34\n",
"TOC: Elapsed: 0:03:1.44\n",
"34\n",
"TOC: Elapsed: 0:03:6.61\n",
"36\n",
"TOC: Elapsed: 0:03:50.57\n",
"34\n",
"TOC: Elapsed: 0:05:15.82\n",
"36\n",
"TOC: Elapsed: 0:07:19.41\n",
"35\n",
"TOC: Elapsed: 0:13:6.09\n",
"34\n",
"TOC: Elapsed: 0:26:38.57\n"
]
}
],
"source": [
"nHs = [16, 30, 32, 40, 50, 60, 64, 70, 80, 90, 110, 128]\n",
"output = []\n",
"\n",
"for nH in nHs:\n",
" output.append(rpd.outerapproximation(nH=nH, tol=1e-4, linprog_method='simplex'))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Number of gradients | \n",
" Area | \n",
" Area relative to min area | \n",
"
\n",
" \n",
" \n",
" \n",
" 0 | \n",
" 16.0 | \n",
" 26.3979 | \n",
" 102.7316 | \n",
"
\n",
" \n",
" 1 | \n",
" 30.0 | \n",
" 26.8183 | \n",
" 104.3676 | \n",
"
\n",
" \n",
" 2 | \n",
" 32.0 | \n",
" 26.1970 | \n",
" 101.9497 | \n",
"
\n",
" \n",
" 3 | \n",
" 40.0 | \n",
" 25.9072 | \n",
" 100.8219 | \n",
"
\n",
" \n",
" 4 | \n",
" 50.0 | \n",
" 26.3712 | \n",
" 102.6276 | \n",
"
\n",
" \n",
" 5 | \n",
" 60.0 | \n",
" 25.6960 | \n",
" 100.0000 | \n",
"
\n",
" \n",
" 6 | \n",
" 64.0 | \n",
" 25.7264 | \n",
" 100.1183 | \n",
"
\n",
" \n",
" 7 | \n",
" 70.0 | \n",
" 26.1435 | \n",
" 101.7415 | \n",
"
\n",
" \n",
" 8 | \n",
" 80.0 | \n",
" 25.8103 | \n",
" 100.4448 | \n",
"
\n",
" \n",
" 9 | \n",
" 90.0 | \n",
" 26.0449 | \n",
" 101.3578 | \n",
"
\n",
" \n",
" 10 | \n",
" 110.0 | \n",
" 25.9620 | \n",
" 101.0352 | \n",
"
\n",
" \n",
" 11 | \n",
" 128.0 | \n",
" 25.7160 | \n",
" 100.0778 | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Number of gradients Area Area relative to min area\n",
"0 16.0 26.3979 102.7316\n",
"1 30.0 26.8183 104.3676\n",
"2 32.0 26.1970 101.9497\n",
"3 40.0 25.9072 100.8219\n",
"4 50.0 26.3712 102.6276\n",
"5 60.0 25.6960 100.0000\n",
"6 64.0 25.7264 100.1183\n",
"7 70.0 26.1435 101.7415\n",
"8 80.0 25.8103 100.4448\n",
"9 90.0 26.0449 101.3578\n",
"10 110.0 25.9620 101.0352\n",
"11 128.0 25.7160 100.0778"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"convex_hulls = np.array([ConvexHull(x[0]) for x in output])\n",
"areas = np.array([round(x.area, 4) for x in convex_hulls])\n",
"relative_areas = np.array([round(area/min(areas)*100, 4) for area in areas])\n",
"\n",
"df = pd.DataFrame(data=np.transpose((nHs, areas, relative_areas)), columns=columns)\n",
"df"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAM0AAADFCAYAAAD373YEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACuhJREFUeJzt3VuMXVUdx/Hvz2nVaWMsUGwoVcuDqSaNUp0QvDVKwXoh\nQIwxkJAgMfpCFH2o0icSn9D6gE8miBeMWi+1VENii3diok2mDLHVUi8I2GlhSnBUcCKl/nyYM9rL\nXM46ew97nzO/TzKZc/bss86/yfl1r73O2mvLNhHRvRc1XUBEv0loIgolNBGFEpqIQglNRKGEJqJQ\nQhNRKKGJKJTQRBRattAOkr4CXA1M2N7Y2XY+8B1gPfAo8EHbf1uordWrV3v9+vUVyo1YHAcOHHjK\n9oXd7KuFptFI2gw8A3z9tNB8Dnja9h2SbgPOs/3phd5sZGTEo6Oj3dQV8YKSdMD2SDf7Ltg9s/0A\n8PRZm68F7uk8vge4rqjCiD7W6znNGtvHO4+fANbMtaOkj0oalTR64sSJHt8uoj0qDwR4un83Zx/P\n9l22R2yPXHhhV13GiFZbcCBgDk9Kusj2cUkXARNVitgzNs6OfUc4NjnF2lXDbNu6ges2XVylyYhF\n0+uR5ofATZ3HNwE/6LWAPWPjbN99kPHJKQyMT06xffdB9oyN99pkxKJaMDSSdgK/BjZIOirpw8Ad\nwFWS/ghc2Xnekx37jjB18tQZ26ZOnmLHviO9NhmxqBbsntm+YY4/bamjgGOTU0XbI5rW+IyAtauG\ni7ZHNK3x0GzbuoHh5UNnbBtePsS2rRsaqihifr2OntVmZpQso2fRLxoPDUwHJyGJftF49yyi3yQ0\nEYVa0T2L9sksjbklNHGOmVkaM186z8zSABIc0j2LWWSWxvwSmjhHZmnML6GJc2SWxvwSmjhHZmnM\nLwMBcY7M0phfpdBIuhX4CCDgS7bvrKWqaFwbZmm0ddi759BI2sh0YC4DngP2SrrP9p/qKi6WrjYP\ne1c5p3kdsN/2v2w/D/wSeH89ZcVS1+Zh7yqhOQS8XdIFklYA7wVeefZOWY0metHmYe+eQ2P7MPBZ\n4H5gL/AQcGqW/bIaTRRr87B3pSFn21+2/Sbbm4G/AX+op6xY6to87F119OwVtickvYrp85nL6ykr\nlro2D3tX/Z7m+5IuAE4Ct9ierKGmCKAdw96zqRQa22+vq5CIfpFpNBGFEpqIQglNRKGEJqJQQhNR\nKKGJKJTQRBRKaCIKJTQRhRKaiEIJTUShhCaiUKXQSPqkpN9JOiRpp6SX1lVYRFv1HBpJFwMfB0Zs\nbwSGgOvrKiyirap2z5YBw5KWASuAY9VLimi3KmsEjAOfBx4HjgN/t33/2ftlYY0YNFW6Z+cB1wKX\nAGuBlZJuPHu/LKwRg6ZK9+xK4C+2T9g+CewG3lJPWRHtVSU0jwOXS1ohScAW4HA9ZUW0V5Vzmv3A\nLuBB4GCnrbtqqiuitaourHE7cHtNtUT0hcwIiCiU0EQUyk2doi81ee+ahCb6TtP3rkn3LPpO0/eu\nSWii7zR975qEJvpO0/euSWii7zR975oMBETfafreNQlN9KUm712T7llEoYQmolBCE1GoypWbGyQ9\ndNrPPyR9os7iItqo54EA20eASwEkDQHjwL011RXRWnV1z7YAf7b9WE3tRbRWXaG5Htg52x+yGk0M\nmsqhkfRi4Brge7P9PavRxKCp40jzHuBB20/W0FZE69URmhuYo2sWMYiqLoC+EriK6TXPIpaEqqvR\nPAtcUFMtEX0hMwIiCiU0EYUSmohCCU1EoYQmolBCE1EooYkolNBEFEpoIgolNBGFEpqIQglNRKGE\nJqJQ1UsDVknaJelhSYclvbmuwiLaquqytF8A9tr+QOey5xU11BTRaj2HRtLLgc3AhwBsPwc8V09Z\nEe1VpXt2CXAC+KqkMUl3d67kPENWo4lBUyU0y4A3Al+0vQl4Frjt7J2yGk0MmiqhOQoctb2/83wX\n0yGKGGg9h8b2E8BfJc3cfmoL8PtaqoposaqjZx8DvtkZOXsEuLl6SRHtVnU1moeAkZpqiegLmREQ\nUSihiSiU0EQUSmgiCiU0EYUSmohCCU1EoYQmolBCE1EooYkolNBEFEpoIgolNBGFKs1ylvQo8E/g\nFPC87cx4joFX9XoagHfafqqGdiL6QrpnEYWqhsbATyQdkPTR2XbIajQxaKqG5m22LwXeA9wiafPZ\nO2Q1mhg0lUJje7zzewK4F7isjqIi2qzn0EhaKellM4+BdwGH6iosoq2qjJ6tAe6VNNPOt2zvraWq\niBbrOTS2HwHeUGMtEX0hQ84RhRKaiEIJTUShhCaiUEITUSihiSiU0EQUSmgiCiU0EYUSmohCCU1E\noYQmolDl0EgakjQm6b46CopouzqONLcCh2toJ6IvVAqNpHXA+4C76yknov2qHmnuBD4F/GeuHbKw\nRgyaKpc7Xw1M2D4w335ZWCMGTZUjzVuBazqrbH4buELSN2qpKqLFeg6N7e2219leD1wP/Mz2jbVV\nFtFS+Z4molAdazlj+xfAL+poK6LtcqSJKJTQRBRKaCIKJTQRhRKaiEIJTUShhCaiUEITUSihiSiU\n0EQUSmgiCiU0EYUSmohCPc9ylvRS4AHgJZ12dtm+vZe29oyNs2PfEY5NTrF21TDbtm7guk0X91pa\nxP8sxmeryqUB/wausP2MpOXAryT9yPZvShrZMzbO9t0HmTp5CoDxySm27z4IkOBEJYv12apy5aZt\nP9N5urzz49J2duw78r9/1Iypk6fYse9Ir6VFAIv32aq6hNOQpIeACeDHtvfPss+8q9Ecm5yate25\ntkd0a7E+W5VCY/uU7UuBdcBlkjbOss+8q9GsXTU8a9tzbY/o1mJ9tmoZPbM9CfwceHfpa7dt3cDw\n8qEztg0vH2Lb1g11lBZL2GJ9tqqMnl0InLQ9KWkYuAr4bGk7MydkGT2Lui3WZ0t28bn79Aul1wP3\nAENMH7G+a/sz871mZGTEo6OjPb1fxGKSdMD2SDf79nyksf1bYFOvr4/oV5kREFEooYko1PM5TU9v\nJp0AHnvB3vD/VgNPNfC+VfRbzf1WL5xZ86ttd7VC/wsamqZIGu32JK8t+q3mfqsXeq853bOIQglN\nRKGlEpq7mi6gB/1Wc7/VCz3WvCTOaSLqtFSONBG1SWgiCi2J0HSu+xmTdF/TtSxE0ipJuyQ9LOmw\npDc3XdNCJH1S0u8kHZK0s3MpfKtI+oqkCUmHTtt2vqQfS/pj5/d53bS1JEID3AocbrqILn0B2Gv7\ntcAbaHndki4GPg6M2N7I9ATe65utalZf49xLV24Dfmr7NcBPO88XNPChkbQOeB9wd9O1LETSy4HN\nwJcBbD/XuVap7ZYBw5KWASuAYw3Xcw7bDwBPn7X5WqZn6tP5fV03bQ18aIA7gU8B/2m6kC5cApwA\nvtrpTt4taWXTRc3H9jjweeBx4Djwd9v3N1tV19bYPt55/ASwppsXDXRoJF0NTNg+0HQtXVoGvBH4\nou1NwLN02WVoSuc84FqmA78WWCnpxmarKufp7166+v5loEMDvBW4RtKjwLeBKyR9o9mS5nUUOHra\nAiW7mA5Rm10J/MX2Cdsngd3AWxquqVtPSroIoPN7opsXDXRobG+3vc72eqZPTn9mu7X/C9p+Avir\npJmL2LcAv2+wpG48DlwuaYUkMV1zqwcvTvND4KbO45uAH3TzoiqLBcbi+BjwTUkvBh4Bbm64nnnZ\n3i9pF/Ag8DwwRgun1EjaCbwDWC3pKHA7cAfwXUkfZvqSlQ921Vam0USUGejuWcRiSGgiCiU0EYUS\nmohCCU1EoYQmolBCE1Hov1+O12CIrBd/AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = output[5][0][:, 0]\n",
"y = output[5][0][:, 1]\n",
"plt.scatter(x,y)\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}