{
 "metadata": {
  "name": "",
  "signature": "sha256:b5ffa093f4dbdba8647b976ddbcb15fb6af6a2769b033d3a25f7757704c39f98"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# mc_compute_stationary_sympy: demonstration"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from mc_compute_stationary_sympy import mc_compute_stationary_sympy"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = [[.4, .6], [.2, .8]]\n",
      "mc_compute_stationary_sympy(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "([1.00000000000000],\n",
        " [array([0.250000000000000, 0.750000000000000], dtype=object)])"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "type(_[1][0][0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "sympy.core.numbers.Float"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "B = np.identity(2)\n",
      "mc_compute_stationary_sympy(B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "([1.00000000000000],\n",
        " [array([1.00000000000000, 0], dtype=object),\n",
        "  array([0, 1.00000000000000], dtype=object)])"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "C = np.identity(3)\n",
      "mc_compute_stationary_sympy(C)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "([1.00000000000000],\n",
        " [array([1.00000000000000, 0, 0], dtype=object),\n",
        "  array([0, 1.00000000000000, 0], dtype=object),\n",
        "  array([0, 0, 1.00000000000000], dtype=object)])"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## KMR Markov matrices"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us try matrices generated by the KMR model\n",
      "where action 1 is risk-dominant, \n",
      "so that the theory says that the unique stationary distribution is close to\n",
      "the distribution that concentrates on the state $N$\n",
      "(where the state is the number of players playing action 1)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from test_mc_compute_stationary import KMR_Markov_matrix_sequential, KMR_Markov_matrix_simultaneous"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "P = KMR_Markov_matrix_sequential(N=3, p=1./3, epsilon=1e-14)\n",
      "eigvals_P, eigvecs_P = mc_compute_stationary_sympy(P)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "eigvals_P, eigvecs_P"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "([1.00000000000000, 1.00000000000000],\n",
        " [array([-428914250225761., -2.14457125112880, 2.14457125112881,\n",
        "         428914250225762.], dtype=object),\n",
        "  array([4.99999999999991e-15, 7.49999999999982e-29, 1.49999999999996e-14,\n",
        "         0.999999999999980], dtype=object)])"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[val.evalf(50) for val in eigvals_P]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "[1.0000000000000017763568394002504646778106689453125,\n",
        " 1.0000000000000048849813083506887778639793395996094]"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[np.dot(vec, P) for vec in eigvecs_P]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "[array([-428914250225759., -2.14457125112881, 2.14457125112881,\n",
        "        428914250225760.], dtype=object),\n",
        " array([4.99999999999991e-15, 7.49999999999987e-29, 1.49999999999997e-14,\n",
        "        0.999999999999980], dtype=object)]"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Q = KMR_Markov_matrix_simultaneous(N=5, p=1./3, epsilon=1e-15)\n",
      "eigvals_Q, eigvecs_Q = mc_compute_stationary_sympy(Q)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "eigvals_Q, eigvecs_Q"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "([1.00000000000000, 1.00000000000000],\n",
        " [array([1.00000000000001, 2.50000000000002e-15, 2.50000000000003e-30,\n",
        "         -2.79568440577553e-44, -2.51812864961921e-29, -9.07251459847683e-15], dtype=object),\n",
        "  array([0, 0, 0, 3.08148791101960e-30, 2.77555756156288e-15,\n",
        "         0.999999999999997], dtype=object)])"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[val.evalf(50) for val in eigvals_Q]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 13,
       "text": [
        "[1.0000000000000024424906541753443889319896697998047,\n",
        " 1.0000000000000026645352591003756970167160034179688]"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[np.dot(vec, Q) for vec in eigvecs_Q]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "[array([1.00000000000001, 2.50000000000003e-15, 2.50000000000003e-30,\n",
        "        -2.67068440577554e-44, -2.51812864961921e-29, -9.07251459847683e-15], dtype=object),\n",
        " array([5.27109897161535e-77, 4.74778387287991e-61, 1.71056941445903e-45,\n",
        "        3.08148791101961e-30, 2.77555756156289e-15, 0.999999999999997], dtype=object)]"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "R = KMR_Markov_matrix_sequential(N=27, p=1./3, epsilon=1e-2)\n",
      "eigvals_R, eigvecs_R = mc_compute_stationary_sympy(R)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "eigvals_R, eigvecs_R"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "([], [])"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "mc_compute_stationary_sympy fails to return any eigenvalue or eigenvector for this $28 \\times 28$ matrix."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## SymPy\n",
      "\n",
      "mc_compute_stationary_sympy uses `Matrix.eigenvects()` from SymPy."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy.matrices import Matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Matrix(P).transpose().eigenvects()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "[(1.00000000000000, 1, [Matrix([\n",
        "   [   -0.999999999999998],\n",
        "   [-4.99999999999999e-15],\n",
        "   [              5.0e-15],\n",
        "   [                  1.0]])]), (1.00000000000000, 1, [Matrix([\n",
        "   [5.00000000000001e-15],\n",
        "   [7.49999999999997e-29],\n",
        "   [1.49999999999999e-14],\n",
        "   [                 1.0]])]), (0.666666666666665, 1, [Matrix([\n",
        "   [ 5.00000000000001e-15],\n",
        "   [-5.00000000000003e-15],\n",
        "   [                 -1.0],\n",
        "   [                  1.0]])]), (-3.99638918680497e-18, 1, [Matrix([\n",
        "   [-0.999999999999998],\n",
        "   [  2.99999999999999],\n",
        "   [              -3.0],\n",
        "   [               1.0]])])]"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Matrix(Q).transpose().eigenvects()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "[(0, 4, [Matrix([\n",
        "   [-1.0],\n",
        "   [ 1.0],\n",
        "   [   0],\n",
        "   [   0],\n",
        "   [   0],\n",
        "   [   0]]), Matrix([\n",
        "   [   0],\n",
        "   [   0],\n",
        "   [-1.0],\n",
        "   [ 1.0],\n",
        "   [   0],\n",
        "   [   0]]), Matrix([\n",
        "   [   0],\n",
        "   [   0],\n",
        "   [-1.0],\n",
        "   [   0],\n",
        "   [ 1.0],\n",
        "   [   0]]), Matrix([\n",
        "   [   0],\n",
        "   [   0],\n",
        "   [-1.0],\n",
        "   [   0],\n",
        "   [   0],\n",
        "   [ 1.0]])]), (1.00000000000000, 1, [Matrix([\n",
        "   [   -110223024625157.0],\n",
        "   [   -0.275557561562892],\n",
        "   [-2.75557561562893e-16],\n",
        "   [ 3.08148791101961e-30],\n",
        "   [ 2.77555756156289e-15],\n",
        "   [                  1.0]])]), (1.00000000000000, 1, [Matrix([\n",
        "   [                   0],\n",
        "   [                   0],\n",
        "   [                   0],\n",
        "   [3.08148791101961e-30],\n",
        "   [2.77555756156289e-15],\n",
        "   [                 1.0]])])]"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Matrix(R).transpose().eigenvects()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "[]"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    }
   ],
   "metadata": {}
  }
 ]
}