{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simplex Rounding Simulations Generator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Libraries, Helpfer Functions, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-22T14:34:54.435368Z",
     "start_time": "2018-07-22T14:34:53.807257Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`polytope` failed to import `cvxopt.glpk`.\n",
      "will use `scipy.optimize.linprog`\n"
     ]
    }
   ],
   "source": [
    "%reset -f\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "\n",
    "import os\n",
    "import pickle\n",
    "import numpy as np\n",
    "from itertools import product\n",
    "\n",
    "from polytope2 import Polytope2\n",
    "from ellipsoid import Ellipsoid\n",
    "from polytope2 import make_standard_simplex\n",
    "from randomwalks import DikinWalk\n",
    "\n",
    "def print_func(lbl, obj):\n",
    "    print(60*'-', '{:^60}'.format(lbl), 60*'-', sep='\\n')\n",
    "    print(obj)\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check One Run - Standard Simplex\n",
    "Take `n_steps = 1000` steps with a given value of r to see if we get roughly 30% of the non-lazy steps being rejected due to the metropolis filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-22T14:34:57.345254Z",
     "start_time": "2018-07-22T14:34:54.436253Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rejection frequency in 1000 steps:\n",
      "    metropolis:     0.3640\n",
      "          none:     0.6360\n"
     ]
    }
   ],
   "source": [
    "# parameters \n",
    "dim = 50\n",
    "n_steps = 1000\n",
    "r = .35\n",
    "\n",
    "# body and walk\n",
    "P = make_standard_simplex(dim)\n",
    "x0 = np.random.dirichlet(np.ones(dim + 1))[:dim]\n",
    "walk = DikinWalk(body=P, r=r, x0=x0, proposal='gaussian', lazy_steps=False)\n",
    "\n",
    "# steps and step types\n",
    "X = np.full((n_steps, dim), np.nan)\n",
    "rejections = np.full(n_steps, 'not avail', dtype='<U10')\n",
    "for i in range(n_steps):\n",
    "    walk.step()\n",
    "    X[i, :] = walk.x.copy()\n",
    "    rejections[i] = walk.reject_type \n",
    "\n",
    "# check rejection rate\n",
    "counts = np.unique(rejections, return_counts=True)\n",
    "print('Rejection frequency in {} steps:'.format(n_steps))\n",
    "for label, count in zip(counts[0], counts[1]):\n",
    "    print('{:>14}'.format(label) + ':', \n",
    "          '{:>10.4f}'.format(count/n_steps))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check One Run - Distorted Simplex\n",
    "Take `n_steps = 1000` steps with a given value of r to see if we get roughly 30% of the non-lazy steps being rejected due to the metropolis filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-22T14:35:00.186253Z",
     "start_time": "2018-07-22T14:34:57.346254Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------------------------------\n",
      "                   new condition number:                    \n",
      "------------------------------------------------------------\n",
      "132650.99998698864\n",
      "\n",
      "Rejection frequency in 1000 steps:\n",
      "    metropolis:     0.3460\n",
      "          none:     0.6540\n"
     ]
    }
   ],
   "source": [
    "# parameters (same as before but now specify kappa)\n",
    "kappa = (dim + 1)**3\n",
    "a = np.ones(dim + 1)\n",
    "a0 = a.sum()\n",
    "c = (1/(a0**2 * (a0 + 1)))\n",
    "sigma_full = np.outer(a=-a*c, b=a)\n",
    "sigma_full[np.diag_indices(dim + 1)] = a*(a0 - a)*c\n",
    "sigma = sigma_full[0:dim, 0:dim]\n",
    "vals, vecs = np.linalg.eigh(sigma)\n",
    "U = vecs.copy()\n",
    "D = np.diag([np.sqrt(a0/kappa)] + [1]*(dim -1))\n",
    "L = np.matmul(U, np.matmul(D, U.T))\n",
    "sigma_new = np.matmul(L, np.matmul(sigma, L.T))\n",
    "new_vals = np.linalg.eigvalsh(sigma_new)\n",
    "print_func('new condition number:', new_vals.max()/new_vals.min())\n",
    "\n",
    "# make new body and walk on new body\n",
    "P_new = P.copy()\n",
    "P_new.A = np.matmul(P.A, np.linalg.inv(L))\n",
    "x0 = np.matmul(L, np.random.dirichlet(a)[:dim])\n",
    "walk = DikinWalk(body=P_new, r=r, x0=x0, proposal='gaussian', lazy_steps=False)\n",
    "\n",
    "# steps and step types\n",
    "X = np.full((n_steps, dim), np.nan)\n",
    "rejections = np.full(n_steps, 'not avail', dtype='<U10')\n",
    "for i in range(n_steps):\n",
    "    walk.step()\n",
    "    X[i, :] = walk.x.copy()\n",
    "    rejections[i] = walk.reject_type \n",
    "\n",
    "# check rejection rate\n",
    "counts = np.unique(rejections, return_counts=True)\n",
    "print('Rejection frequency in {} steps:'.format(n_steps))\n",
    "for label, count in zip(counts[0], counts[1]):\n",
    "    print('{:>14}'.format(label) + ':', \n",
    "          '{:>10.4f}'.format(count/n_steps))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulations\n",
    "Let's first generate the parameters over which we'll run simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-22T14:35:00.254254Z",
     "start_time": "2018-07-22T14:35:00.187254Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'d': 50, 'k': 4},\n",
       " {'d': 50, 'k': 3},\n",
       " {'d': 50, 'k': 2},\n",
       " {'d': 50, 'k': 1},\n",
       " {'d': 40, 'k': 4},\n",
       " {'d': 40, 'k': 3},\n",
       " {'d': 40, 'k': 2},\n",
       " {'d': 40, 'k': 1},\n",
       " {'d': 30, 'k': 4},\n",
       " {'d': 30, 'k': 3},\n",
       " {'d': 30, 'k': 2},\n",
       " {'d': 30, 'k': 1},\n",
       " {'d': 20, 'k': 4},\n",
       " {'d': 20, 'k': 3},\n",
       " {'d': 20, 'k': 2},\n",
       " {'d': 20, 'k': 1},\n",
       " {'d': 10, 'k': 4},\n",
       " {'d': 10, 'k': 3},\n",
       " {'d': 10, 'k': 2},\n",
       " {'d': 10, 'k': 1}]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def expand_params(params):\n",
    "    return [dict(zip(params.keys(), v)) for v in product(*params.values())]\n",
    "\n",
    "d_vals = [50, 40, 30, 20, 10]\n",
    "k_vals = [4, 3, 2, 1]\n",
    "param_vals = {'d' : d_vals, 'k' : k_vals}\n",
    "param_grid = expand_params(param_vals)\n",
    "param_grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-23T00:42:42.500703Z",
     "start_time": "2018-07-22T14:35:00.255253Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'d': 50, 'k': 4}\n",
      "{'d': 50, 'k': 3}\n",
      "{'d': 50, 'k': 2}\n",
      "{'d': 50, 'k': 1}\n",
      "{'d': 40, 'k': 4}\n",
      "{'d': 40, 'k': 3}\n",
      "{'d': 40, 'k': 2}\n",
      "{'d': 40, 'k': 1}\n",
      "{'d': 30, 'k': 4}\n",
      "{'d': 30, 'k': 3}\n",
      "{'d': 30, 'k': 2}\n",
      "{'d': 30, 'k': 1}\n",
      "{'d': 20, 'k': 4}\n",
      "{'d': 20, 'k': 3}\n",
      "{'d': 20, 'k': 2}\n",
      "{'d': 20, 'k': 1}\n",
      "{'d': 10, 'k': 4}\n",
      "{'d': 10, 'k': 3}\n",
      "{'d': 10, 'k': 2}\n",
      "{'d': 10, 'k': 1}\n"
     ]
    }
   ],
   "source": [
    "r = 0.35\n",
    "n_steps = 10**6\n",
    "\n",
    "for param in param_grid:\n",
    "    # get params\n",
    "    print(param)\n",
    "    d = param['d']\n",
    "    k = param['k']\n",
    "    kappa = (d + 1)**k\n",
    "    \n",
    "    # get moments and transformation\n",
    "    a = np.ones(d + 1)\n",
    "    a0 = a.sum()\n",
    "    c = (1/(a0**2 * (a0 + 1)))\n",
    "    sigma_full = np.outer(a=-a*c, b=a)\n",
    "    sigma_full[np.diag_indices(d + 1)] = a*(a0 - a)*c\n",
    "    sigma = sigma_full[0:d, 0:d]\n",
    "    vals, vecs = np.linalg.eigh(sigma)\n",
    "    U = vecs.copy()\n",
    "    D = np.diag([np.sqrt(a0/kappa)] + [1]*(d - 1))\n",
    "    L = np.matmul(U, np.matmul(D, U.T))\n",
    "\n",
    "    # make new body and walk on that body\n",
    "    P = make_standard_simplex(dim=d)\n",
    "    P_new = P.copy()\n",
    "    P_new.A = np.matmul(P.A, np.linalg.inv(L))\n",
    "    x0 = np.matmul(L, np.random.dirichlet(a)[:d])\n",
    "    walk = DikinWalk(body=P_new, x0=x0, r=r, proposal='gaussian', \n",
    "                     lazy_steps=False)\n",
    "\n",
    "    # steps and step types\n",
    "    X = np.full((n_steps, d), np.nan)\n",
    "    rejections = np.full(n_steps, 'not avail', dtype='<U10')\n",
    "    for i in range(n_steps):\n",
    "        walk.step()\n",
    "        X[i, :] = walk.x.copy()\n",
    "        rejections[i] = walk.reject_type \n",
    "    result = param.copy()\n",
    "    result.update({'X':X, 'rejections':rejections})\n",
    "    filename = os.path.join(os.path.normpath('E:/Data/Rounding'), \n",
    "                            'simplex_d{}_k{}.pkl'.format(d, k))\n",
    "    with open(filename, 'wb') as file:\n",
    "        pickle.dump(result, file, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  }
 ],
 "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.5"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}