Bellman_Error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$Belman \\quad Error = BE = \\Delta_R  + \\gamma \\Delta_P W_{\\phi} $$\n",
    "\n",
    "\n",
    "Считаю, что $Ф = [U_1, r]$ \n",
    "\n",
    "in tabular case we need add extra column: reward, because otherwise reward_error will be very bad\n",
    "\n",
    "# Что если ортогонализовать K(P, r)???\n",
    "## Странности:\n",
    "1. Бустится kry если его ортогонализовать\n",
    "2. Формула в общем случае для вычисления Bellman Error is not valid, because F @ reward_A is not valid by shape\n",
    "3. \n",
    "\n",
    "\n",
    "\n",
    "$\\textbf{Case 1: Tabular raw features}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 100\n",
    "states = [i for i in range(n)]\n",
    "transition_matrix_rank = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class mdp_low_rank():\n",
    "    def __init__(self, A, A_next=None, reward=None, gamma=0.95, method=\"svd\", P=None, r_start=10):\n",
    "        self.A = A\n",
    "        self.A_next = A_next\n",
    "        self.n = A.shape[0]\n",
    "        self.reward = reward\n",
    "        if (self.reward is None):\n",
    "            self.reward = np.random.random(self.n)\n",
    "        self.num_iter = 10\n",
    "        self.gamma = gamma\n",
    "        self.BE = []\n",
    "        self.epsilon = 1e-8\n",
    "#        self.r_list = [i for i in range(10, self.n - 1)]\n",
    "        self.r_list = [i for i in range(r_start, self.n + 1)]\n",
    "        self.transition_matrix_rank = 40\n",
    "        self.P = None\n",
    "        self.method = method\n",
    "        self.max_iter = 10\n",
    "        self.states = np.arange(0, self.n)\n",
    "        self.P = P\n",
    "        if (self.A is not None and self.A_next is not None):\n",
    "            self.P_A = np.linalg.pinv(self.A) @ self.A_next\n",
    "        if (self.A is not None and self.reward is not None):\n",
    "            self.reward_A = np.linalg.pinv(self.A) @ self.reward\n",
    "        \n",
    "        \n",
    "    def check_transition_matrix(self, P):\n",
    "        r = self.transition_matrix_rank\n",
    "        assert(np.linalg.norm(P.sum(axis=1) - np.ones(P.shape[0])) < self.epsilon)\n",
    "\n",
    "        U, S, V = np.linalg.svd(P, full_matrices=False)\n",
    "        P = (U * S) @ V\n",
    "\n",
    "        assert(np.linalg.norm(S[r:] - np.zeros(P.shape[0] - r)) < self.epsilon)\n",
    "\n",
    "       \n",
    "    def create_transition_matrix(self):\n",
    "        \n",
    "        tmp_matr1 = np.random.random((self.n, self.transition_matrix_rank))\n",
    "        tmp_matr2 = np.random.random((self.transition_matrix_rank, self.n))\n",
    "        P = tmp_matr1 @ tmp_matr2\n",
    "        \n",
    "        sum_rows = P.sum(axis=1)\n",
    "        for j in range(n):\n",
    "            P[j, :]  = P[j, :] / sum_rows[j]\n",
    "\n",
    "        self.check_transition_matrix(P)\n",
    "        \n",
    "        return P\n",
    "        \n",
    "        \n",
    "    @staticmethod \n",
    "    def sampling_states(P, size, init_state = None):\n",
    "        \"\"\"\n",
    "        init_states\n",
    "        return list of sampled states\n",
    "        states in {0, ... , n}\n",
    "        \"\"\"\n",
    "        \n",
    "        states_list = []\n",
    "        if (init_state == None):\n",
    "            init_states = np.random.choice(self.states)\n",
    "            \n",
    "        cur_state = init_state\n",
    "            \n",
    "        for i in range(size):\n",
    "            p = P[cur_state]\n",
    "            next_state = np.random.choice(self.states, size=1, p=p)\n",
    "            states_list.append(next_state)\n",
    "            cur_state = next_state\n",
    "            \n",
    "            \n",
    "        return states_list\n",
    "\n",
    "     \n",
    "    \n",
    "    def choose_feature_matrix(self, method, r):\n",
    "        if (method == \"svd\"):\n",
    "            U, S, V = np.linalg.svd(self.P, full_matrices=True)\n",
    "#           F = np.c_[U[:, :r] @ np.diag(S)[:r, :r], self.reward]\n",
    "            F = np.c_[U[:, :r] @ np.diag(S)[:r, :r], self.reward]\n",
    "        elif (method == \"eig-r\"):\n",
    "            w, vr = np.linalg.eig(self.P)\n",
    "            F = np.c_[vr[:, :r], self.reward]\n",
    "        elif (method == \"kry\"): \n",
    "            F = np.zeros((self.n, r))\n",
    "            F[:, 0] = self.reward\n",
    "            for k in range(1, r):\n",
    "                F[:, k] = self.P @ F[:, k - 1]\n",
    "            F, _ = np.linalg.qr(F) # cheat sheet\n",
    "        elif (method == \"lfd\"):\n",
    "            D = np.random.rand(r, A.shape[1])\n",
    "            E = np.random.rand(A.shape[1], r)\n",
    "            diff, i = np.inf, 0\n",
    "            while (np.linalg.norm(E @ D - A) < diff and i < self.max_iter):\n",
    "                diff = np.linalg.norm(E @ D - A)\n",
    "                E = A @ np.linalg.pinv(D)\n",
    "                D = np.linalg.pinv(E) @ A\n",
    "                i += 1\n",
    "            return E\n",
    "\n",
    "            \n",
    "        else:\n",
    "            raise ValueError(\"This method is not implemented\")\n",
    "        \n",
    "        \n",
    "        return F\n",
    "        \n",
    "    def get_approx_transition_matrix(self):\n",
    "        if (self.P_A is None):\n",
    "            self.P_A = np.linalg.pinv(self.A) @ self.A_next\n",
    "            \n",
    "        return self.P_A\n",
    "    \n",
    "    \n",
    "    def get_approx_reward(self):\n",
    "        if (self.reward_A is None):\n",
    "            self.reward_A = np.linalg.pinv(self.A) @ self.reward\n",
    "            \n",
    "        return self.reward_A\n",
    " \n",
    "\n",
    "    def get_approx_feature_matrix(self, method, r):\n",
    "        if (self.method == \"svd\"):\n",
    "            P_A = self.get_approx_transition_matrix()\n",
    "            reward_A = self.get_approx_reward()\n",
    "            U, S, V = np.linalg.svd(P_A, full_matrices=True)\n",
    "            F = np.c_[U[:l, :(r - 1)] @ np.diag(S)[:(r - 1), :(r - 1)], reward_A]           \n",
    "            \n",
    "#            import pdb; pdb.set_trace();\n",
    "#            pdb.set_trace()\n",
    "            \n",
    "            \n",
    "            return self.A @ F, F\n",
    "#           return F\n",
    "        \n",
    "        elif (method == \"kry\"): \n",
    "            F = np.zeros((self.n, r))\n",
    "            F[:, 0] = self.reward\n",
    "            for k in range(1, r):\n",
    "                F[:, k] = self.P @ F[:, k - 1]\n",
    "            F, _ = np.linalg.qr(F) # cheat sheet\n",
    "#            reward_A = \n",
    "            \n",
    "            return self.A @ F, F\n",
    "        \n",
    "        \n",
    "   \n",
    "    \n",
    "    def calculate_BE(self, r, method):\n",
    "        if (method == \"tabular\"):\n",
    "            F = self.choose_feature_matrix(method=self.method, r=r)\n",
    "            P_F = np.linalg.pinv(F) @ self.P @ F\n",
    "            reward_F = np.linalg.pinv(F) @ self.reward\n",
    "            w_F = np.linalg.pinv(np.eye(P_F.shape[0]) - self.gamma * P_F) @ reward_F\n",
    "            reward_error = np.linalg.norm(self.reward - F @ reward_F)\n",
    "            transition_error = np.linalg.norm(self.P @ F - F @ P_F)\n",
    "            all_error = np.linalg.norm((self.reward - F @ reward_F) + self.gamma * (self.P @ F - F @ P_F) @ w_F)\n",
    "        elif (method == \"approx\"):\n",
    "            F, F_hat = self.get_approx_feature_matrix(method=self.method, r=r)\n",
    "            P_F = np.linalg.pinv(F) @ self.P_A @ F\n",
    "#            P_F = np.linalg.pinv(F_hat) @ self.P_A @ F_hat ## smth strange\n",
    "#            reward_F = np.linalg.pinv(F) @ self.reward_A ## smth strange\n",
    "#            reward_F = np.linalg.pinv(F_hat) @ self.reward_A \n",
    "            reward_F = np.linalg.pinv(F) @ self.reward_A \n",
    "            w_F = np.linalg.pinv(np.eye(P_F.shape[0]) - self.gamma * P_F) @ reward_F\n",
    "            reward_error = np.linalg.norm(self.reward - F @ reward_F)\n",
    "            transition_error = np.linalg.norm(self.P @ F - F @ P_F)\n",
    "            all_error = np.linalg.norm((self.reward - F @ reward_F) + self.gamma * (self.P @ F - F @ P_F) @ w_F)\n",
    "        \n",
    "        \n",
    "        return all_error\n",
    "    \n",
    "\n",
    "    def calculate_BE_list(self, method=\"tabular\"):\n",
    "        \n",
    "            \n",
    "        if (method == \"tabular\"):\n",
    "            self.P = self.create_transition_matrix()\n",
    "            for r in self.r_list:\n",
    "                self.BE.append(self.calculate_BE(r, method=method))\n",
    "\n",
    "            return self.BE, self.r_list\n",
    "        \n",
    "        if (method == \"approx\"):\n",
    "            if (self.A is None or self.A_next is None or self.reward is None or self.P is None):\n",
    "                raise ValueError(\"This method is not supported!\")\n",
    "                \n",
    "            for r in self.r_list:\n",
    "                self.BE.append(self.calculate_BE(r, method=method))\n",
    "\n",
    "            return self.BE, self.r_list\n",
    "\n",
    "                \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n = 100 # size of observed states\n",
    "l = 50  # size of features for every state\n",
    "\n",
    "# A = np.eye(n)\n",
    "A = np.random.binomial(n=1, p=0.5, size=(n, n))\n",
    "mdp = mdp_low_rank(A=A, method=\"eig-r\", r_start=1)\n",
    "BE_eig_r, r_list = mdp.calculate_BE_list()\n",
    "mdp = mdp_low_rank(A=A, method=\"kry\", r_start=1)\n",
    "BE_kry, r_list = mdp.calculate_BE_list()\n",
    "mdp = mdp_low_rank(A=A, method=\"svd\", r_start=1)\n",
    "BE_svd, r_list = mdp.calculate_BE_list()\n",
    "mdp = mdp_low_rank(A=A, method=\"lfd\", r_start=1)\n",
    "BE_lfd, r_list = mdp.calculate_BE_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f414d960d30>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABFOUlEQVR4nO3deZwcdZ34/9e7qu+5j54ck0wmFzlIAjkIpwgCEgEBFTWg4CJ+s+wufl1XWfWrLrio6O7+XFFYlBVEQTkED+QQOeS+QghnEnJMrsk1k7mvvqo+vz+qZzJJZpLJZHq6Z+b9fDz6MdPV3VXvmpqud33OEmMMSimlxiYr2wEopZTKHk0CSik1hmkSUEqpMUyTgFJKjWGaBJRSagzzZTuAI1FeXm6qq6uzHYZSSo0oq1at2muMifb12ohKAtXV1bz++uvZDkMppUYUEdna32taHaSUUmOYJgGllBrDNAkopdQYpklAKaXGME0CSik1hmU1CYjINBG5XUQeyGYcSik1Vg15EhCRO0SkTkTePWD5MhF5X0Q2isjXAYwxNcaYq4Y6BqWUUgOTiXECdwI3A7/uXiAiNnALcA5QC6wUkYeMMWsysP2DJLZto+WPf+zjFQGR/Zf4/UgoiBUKYVwXk0hgkknEspFgECsURPx+sH2IbYFl9/wEMMkkJpkEY7DyIliRCOL396zHpFIHbB/EssCykEAQKxhAAgGMMeAaMK63DZ+N2PZB8e5blfTx+wHvFRCfz1uPZUH3NOIi3rLu11wX47rgumAM3dONi4j3Ocvy/k7dD0trFZUaqYY8CRhjnhOR6gMWLwU2GmNqAETkXuAi4LBJQERWACsAqqqqBhXT7o1v037rrUjvk6LeR2HodCcFn29fIrEtSKYwiQRuOiki3hGw8vKwy8rwlZSA34eJxXHjMXBcsC3E8j4vkk44tpeAJRhEbBvjOt57AfF1J699SRnb8hKo42CMC4b9E146kWFMerkB0gmue7vdcfSV4CSdvEW8z9k2EvAjgQDi8/fkXhHx/i6BANh2T1I3jrtfPFY4hJWXhxUOYxw3fbGQRHx+76IjGPT+trbt/ez5P06vo2ddlnex0H0cuo+JiLfORBIwXpz+gLfuUMh73t/FhRr1hmvEcCWwvdfzWuBEESkDvgcsFJFvGGNuPPCDxpjbgNsAlixZMqgzd93c8Vz5dR8/O/tnnFp56oHr7/kCGGO8L0sshhuLIZaV/sL4vS9nIo6Jxbz3OC44Ke+K2XG859DzfgC3qxPT2YmbSGAFAumTgW/fd9hNn4CMwTgOJpH0thGPeycfSZdUXINJJcFx+t7BXgmt5yZB3Se+3l9u42IcJx2v450QRfbtQ8oB14H0yRfLe53ukx0m/V4Xk0qlTywJL7bu5ynHO0mnHMTvQ/wBJOD31pk+4Trt7TgNjaQaG6ArhoRD+IuKvBOl4+w7yRvjncSTKZy2Vsxeb1tekvBKXqRSXizd+9VdgrGs/Uou3sFJ/827SzqCd8IX8bbV/Xcw3j7iOBz0D9edOHoSCOljl/C2OxKJeMnA5+t5HFSyTCdPKz8fu7AQq7BwX4lSxPsuJL1jgW15/+/+AFZRIb5oFF806pWuHe/YuvEYbnsHblsbYLCLirCKirDz8rxEbokXiz+QTqT+fSVWsfaLzYqEsSIRb/3G9JRgrbw8b780wR1SVqeNMMY0AFdnejvVRdUAbGndclAS6P0PIiLeiToQwC4s7GNNeRmMUo10Jp2QerhuT/WgSaW8K3Xb2nfiBC/5dHXhdnbidnZ61Yx+P+KzvfXF47ixmJfsHOeA6kR61tOdzE0q5SWvZDJ9Uk6CSV+cdJcKEgnvEY/hxuK4sS5MLJ6O30voPfvUk/AAx8Fpb8dtbSW5c6d34k9XW4ptexc/Ph84DqlkEhOP4zQ347S09P9HS18wkUwOwRHog8+HnZ/vJYN09SVWr+98IIAVCmOFQvjGjSM4YzqB6dPxjxuHFQ4j4QhWOOSVxkZpiWm4ksAOYHKv55PSy4ZFWaiMPH8eW1v7nT5DqaPWcxWt9uMmEjj19bjxhJcEbRsrGMQqKECCQQBMVxdOczNuZ2dPac6kHEwy4ZWQk8l9JVnX9KoRM7jpJGq6uvCq57wkazo7cVrbvFJk3CtFmkRyX8nZGK+6Mh7D6Wgn9sJ6Wv7wh0PuiwQCSCiE1V09mU4OVijklZIK8rHy8rzqOJ9/Xzted/LorkbsKV2DBIOUfPpT+KJ9zu+WccP1H7sSmCkiU/FO/suBy4Zp24gI1YXVbGnZMlybVEqlWYEAVmXlId8jEa8TRbY5LS3EN20i1dCQLqF1eSWleMIrlcVj+/+e8Eo8blcXTlMTyW3bcDo7IOXsq6bslXQOrEoErzNJ+9NPM+Xuu7LyNxjyJCAi9wBnAOUiUgtcZ4y5XUSuAR4HbOAOY8x7Q73tQ5lSOIU3694czk0qpUYYu6iIyKJFw7rNtmeeofYf/4mdX/s6lTf9eNh72w351owxlxpjJhhj/MaYScaY29PLHzXGHGOMmW6M+d5Qb/dwqgur2dWxi1gqNtybVkqpfhWccQYV/3otbU88Qf2Pbxr27Y+ZDt7VRdUYDNvatmU7FKWU2k/p5z5H8ac+RcNtt9Hypz8N67bHTBKYUjgFQBuHlVI5R0QY/+1vETnhBHb923XE1gzLOFpAk4BSSuUE8fup/PF/Y5eUUPvF/0uqqWlYtjtmkkCeP4+KcAWbWzZnOxSllOqTr6yMST+5iVRdHTu/8lWvd1GGjZkkADClaIqWBJRSOS28YAHjr/s3Ol56iab77sv49sZWEijUJKCUyn3Fl1yChMMkazM/pnZMJYHqwmqa4800x5qzHYpSSh2SBALenFQZNuaSAHhzCCmlVC4Tv9+bLiPDxlQS0B5CSqmRQpNABlQWVOITnyYBpVTO0ySQAX7Lz6SCSVodpJTKeRLwa5tAJkwpnKJJQCmV88Qf0JJAJkwpnMK21m24ZoTeBUopNSZodVCGTCmcQtyJs7tjd7ZDUUqpfolfq4MyYmrRVEC7iSqlcpsEtCSQEeMj4wGo66zLciRKKdU/rQ7KkLJwGQANXQ1ZjkQppfo3JhqGRWSaiNwuIg8M1zYj/ghhX5i9XXuHa5NKKXXEcr5NQETuEJE6EXn3gOXLROR9EdkoIl8/1DqMMTXGmKsGG8NglYXKaIhpSUAplbuGq03gaG40fydwM/Dr7gUiYgO3AOcAtcBKEXkI7+byNx7w+c8bY7JSMV8WLqOxqzEbm1ZKqQEZrjaBQScBY8xzIlJ9wOKlwEZjTA2AiNwLXGSMuRG4YDDbEZEVwAqAqqqqwYa7n/JwuU4doZTKaSO1YbgS2N7reW16WZ9EpExEfgYsFJFv9PUeY8xtxpglxpgl0Wh0SIIsC5Vpm4BSKqcN11TSR1MddNSMMQ3A1cO93bJwGc3xZpJuEr/lH+7NK6XUYY3UksAOYHKv55PSy3JKWcjrJtoUG54bOSul1JHK+TaBfqwEZorIVLyT/3LgsiHexpGrWwsv/DdYPhCLspTXKNzwys1UFM8Eu1dpwBeCYD4ECiC/AgorwRfIUuBKqbFK/H5wHIzjILadse0MOgmIyD3AGUC5iNQC1xljbheRa4DH8XoE3WGMeW9IIj0asRbY9goYF1yHcsuB0iANK38OXbHDfFigYAJMWgyzzodjzoVI6bCErZQauyTgXXyaZDI3k4Ax5tJ+lj8KPDroiDKh6iT457d7npa1boc/nMfeC2+C8SeDcbwXjIFUDOLtEG+Dtl3Qsh2atsLmZ2Htn0FsKJoE/rBXasgfB9FjIDobCsYD4q1LxPtdLLBssPxeicP2gx3wfvrCkBcFO6tNM0qpHCR+r4bCJJMQCmVsO2Py7NMzdUSqA0qmDOxDrgu7VsP7j0HzNkh2QjIGrTug5hlw4oMLRizIq4DiKpj+IZh9HoxfkE4iSqmxar8kkEFjMgl0Tx1xRKOGLQsqF3uPA7kONG2Bjl7dTo0LmJ4qKNwkOCkvWThJ75HsgLY90LYT9m6AZ38Iz/4A8sdDuNhbj1hegiifCeXHQMlUL3EVVnolDKXUqNSTBDLcTXRMJgFITx0xVJPIWTaUTfceR6O9Htb/BTY/t69k4aS8BLPpb/uXNiyflxxKqr3EEJ0F4+bBuGP3JRCl1IjVu00gk8ZuEgjn4PxB+VFYdLn3OJDreNVQTVugeav3s/uxY5XX+N3NDu5rkyibAccthwWf8no7KaVGBK0OyrCyUBnb2rZlO4yBs2woneo9DmQMtO2GPe/C7ncg1uwtMy5sfQn++k144t+8xmvb5zVuW3avn5aXNMTq9Ui/ZtleqcMX9qqhSqd7pY9wMQTyIVjgdalVSg0pTQJDKJ5yaOlKEvLbBH0WAduiLFzG6rrV2Q5taIhA4QTvMfOcg1+vWwdv3wv173slCuP0+umCm/KSBia93N3/NTflNYS/tdN7z4GCRVBSBaXT4IxvQMWcTO+xUqOetgkMoVVbmrjsF6/utyxvXDNWSTMf+I8nyQ8GyQ/a5AV9FIb8FIZ9FIX9RAI+Qn6bsN8m4LO8h23hs8S7cBYQESwRBLBEsCzvp20JAdvCb1sEfILPsvDZ3k8RryOpZQl5AR8hv4VksjdQxWw4+/qjX08ytq86KtYKiXavGqql1lu2+TnY9ip84QmvvUIpNWjaJjCEpkbzuOHiecSTDvGUSzzp8HbrBlZ1GI6dbJNKhOmIp2jsSLBlbwctXUlaYykct4+r3gywLSEvYBMJ+IgEbEJ+G58tiAi2gM+ysC0vsViWt8y2hJDfJj/oIy/96E5kkYCXuLySj51OPpJOSF5iKssPUBg6wnmT/CEvoVTM7vv1PWvgjmVw18fhqr/qoDqljsK+koAmgaM2oSjM5SftPx7gya3bWfUMfOnD45lTdnD1hTEmnTBcupIOiZRLwvGSiOMar8odcE36d2NwjffcNYaUY0g6LknHJeEYUo5LyjGkXINJV6k4rqEj7tAeT9IeS9GVdOhMOMSSDo5rcAy4rvF+dw3xlLPfsljSoS2eoiOeojPhHNHfRATmTSzilOllnDS9jMVTSo48KRxo3Fy49B6462Pw20/Bud/32hPsgDe4LpAPgTzwR7wut0qpfmmbQIb1DBjrp4eQiHelHfLbFJH7M406rqEzkaIj7tCVdOhKOHQlUz1JK+UYEumklHRctuzt5OWaBu54cTM/f64GS2D2+ELmVxYRCXr7nRewiRYEqSgMUZ4X7DlvC4LP3lflFS0IEvKnxyxUnwqf+AXcfwXc3kf7RHoNXoNyAfiC+xqjkX2D5HwhbzR1XtRrhLb9Xq+n6CyYf0mG/5pKZd++6iBtE8iI8lA5MHpuOG9bQkHIT8ERXM1/GehKOKze1sRrWxpZuaWRp9bVEU86xFIOSWfg1WEVBUEmFocJ+iwgSmX05xyX38KsaJgZ5UHKgi6S6PDaERId3rQcsVZwEl4jtHHTA+zSEp3QUe9N/hdr8d7XPU4i3gZLrhxwbEqNRFodlGGHKwmMFeGAzSkzyjllRvlBr8VTDvVtcfa0xmloj/f0CzLG4LiQcr3qsl0tMXY0d7KrJUYi5Z3It1qT+WtNMe1rUwCURPzMmVDNnAmFnDWngpOnlR15Y7jrwG8+CY9eCxVzoerEo9l1pXKaVgdlWM/UEaOkJJAJQZ/NpJIIk0oig/q84xrW7W7lja1NvLujlbW7W7n7la3c/sJm5lUWsuL06Zw3bzw+e4DtA5YNl9wOt50J918OK56BwomDik2pXCd+7R2UcaWhUr3NZAbZlnDsxCKOnVjUsyyWdPjD6h3873M1/N97VnPLuAL+85MLWDCpeGArDZd4jc+/ONtrd7jqCZ1sT41KEhiecQJjuotGebh8zFcHDbeQ3+bSpVU8+S8f5H8+s4jmrgQX3/IiP3hsHbHkAHs4VcyB074MtSuhszGzASuVJVodNAxG3NQRo4hlCefNn8CpM8r5/iNr+dmzm7j9hRoqi8NMLo0ws6KAE6eVcuLUUoojfdzZLVzi/XQzf/s9pbJBq4OGQVm4jDfr38x2GGNaUdjPDy9ZwMULK3luQz3bGzvZ3tjJb1/byh0vbgZgaXUp/37xscweX7jvg923BHU0CajRabiqg7KaBERkDvAloBx4yhhz63BuvyxcRlOsiZSbwmeN6XyYdSdPL+Pk6WU9zxMpl7drm3l5UwN3vrSFj/70Bf7xjBn805kzCPgsbwAaeF1HlRqFhqs6aNBtAiJyh4jUici7ByxfJiLvi8hGEfn6odZhjFlrjLka+BRw6mBjGazyUDkGQ1Osabg3rQ4j4LNYUl3KF8+ayRP/8kHOmz+Bm57awGX/+wrGGG8kMniT2yk1CokI+P25mwSAO4FlvReIiA3cAnwEmAtcKiJzRWS+iDx8wKMi/ZkLgUfIwn2JdazAyFCaF+Cm5Qv54odm8PrWJho7Er2qg7QkoEYvGYYkcDQ3mn9ORKoPWLwU2GiMqQEQkXuBi4wxNwIX9LOeh4CHROQR4LeDjWcwepKAjhUYEY6fXAzAloZOynqqg7RNQI1e4vePuDaBSmB7r+e1QL/DOkXkDODjQJB+SgIisgJYAVBVNbTTE5eFvCSgYwVGhilleQBsa+xgcb42DKvRTwI5XBIYCsaYZ4BnDvOe24DbAJYsWTKkczuXh72pEjQJjAyTS8OIwJa9nVCUTgLaRVSNYsNRHTTUg8V2AJN7PZ+UXpaTIv4Ief48TQIjRNBnM7EozNaGDu0iqsaEkZgEVgIzRWSqiASA5cBDQ7yNIRUNR6nvqs92GGqAqkojbG3s7NVFVJOAGr2Go03gaLqI3gO8DMwSkVoRucoYkwKuAR4H1gL3G2PeG5pQM6MsXEZ9pyaBkaK6PMLWhs5eXUQ1CajRSwKB3G0TMMZc2s/yR8lCd8/BioajrGlYk+0w1ABNKcujsSNBu2ORD9pFVI1qI7E6aMQpD5drddAIMqXUm9Z6Z2t6kJijg8XU6KVJYBhEI1G6Ul10JDuyHYoagO5uott7koCWBNToZQUCudsmMFpEw1EAbRcYIarKvJLA9pb01ZG2CajRTEsCmdc9VkCrhEaG/KCP8vwgW5vTXwztHaRGMa0OGgbdJQEdKzByVJdF2NqULiJrElCjWE53ER0tohFNAiNNVVmEmsbukoC2CajRyxqGLqJjPgkUBgoJWAGtDhpBqsvy2NGWbhjWNgE1muXyLKKjhYhQHi5nb6eWBEaKKWURktjekz6qg5LJJLW1tcRisWGObHiFQiEmTZqEP33zETX65PRU0qNJeUTHCowkXjdRwRUfVh9JoLa2loKCAqqrq70bc4xCxhgaGhqora1l6tSp2Q5HZYi2CQyTaDiqbQIjSHW6m6gjvj7bBGKxGGVlZaM2AYBXgi0rKxv1pZ2xTtsEhomOGh5ZiiMBCkM+Uvj6vb3kaE4A3cbCPo552kV0eETDUVriLSS0p8mIUV2e57ULjLBjdt5559Hc3JztMNQIIX4/uC7GcTK2DU0C6M1lRqIpZXnEjT3ixgk8+uijFBcXD+qzxhhc1x3agFROk3SjfybbBTQJsG+sgFYJjRxTSiPEXBs3lbslgbvvvpulS5dy/PHH8/d///c4jkN1dTV793oXGzfccAOzZs3itNNO49JLL+W//uu/DlrHli1bmDVrFldccQXz5s1j+/btB71HjV5WwLtvRiarhLR3EL1KAtpNdMQozw+QNDbJZILgId73nT+/x5qdrUO67bkTC7nuo8ce8j1r167lvvvu48UXX8Tv9/OP//iP/OY3v+l5feXKlTz44IO89dZbJJNJFi1axOLFi/tc14YNG/jVr37FSSedNKT7oUaA7pKAJoHM6plETksCI0bAZ5Mid0sCTz31FKtWreKEE04AoKuri4qKip7XX3zxRS666CJCoRChUIiPfvSj/a5rypQpmgDGqOGoDtIkAJSGSrHE0iQwggR8Fkl8mMM0DB/uij1TjDF87nOf48Ybb9xv+Z133nnIz23fvr0nIVx99dUsW7aMvLy8TIWpctxwVAdpmwBgWzaloVJtGB5BvCRgY1K52TB81lln8cADD1BXVwdAY2MjW7du7Xn91FNP5c9//jOxWIz29nYefvhhACZPnsybb77Jm2++ydVXX52V2FXukNFeHSQiZwA3AO8B9xpjnslWLNFwVO8pMIIE7IGVBLJl7ty5fPe73+XDH/4wruvi9/u55ZZbel4/4YQTuPDCC1mwYAHjxo1j/vz5FBUVZTFilZNyOQmIyB3ABUCdMWZer+XLgJsAG/iFMeYHh1iNAdqBEFA72FiGQnm4XEsCI0jQb5EyPkwO317y05/+NJ/+9Kf3W7Zly5ae37/61a9y/fXX09nZyemnn95nw3B1dTXvvvtupkNVOSrX2wTuBG4Gft29QERs4BbgHLyT+koReQgvIdx4wOc/DzxvjHlWRMYBPwI+cxTxHJVoJMq6xnXZ2rw6QkHbGpGDxXpbsWIFa9asIRaL8bnPfY5FixZlOySVY3K6i6gx5jkRqT5g8VJgozGmBkBE7gUuMsbciFdq6E8T9N3TT0RWACsAqqqqBhvuYZWHy2mINeC4DrZlZ2w7amgEfBZd9D130Ejx29/+NtshqBw3HG0CQ90wXAn0Hs1Sm17WJxH5uIj8HLgLr1RxEGPMbcaYJcaYJdFodEiD7S0ajuIal6Z4U8a2oYZOd+8g0fsJqFFs1DcMG2N+D/w+mzF0633D+e7BYyp3BXwWKUbetBFKHZEROG3EDmByr+eT0styXnlEbzg/kgRsi4SWBNQoNxLHCawEZorIVBEJAMuBh4Z4GxmhN5wfWQI+i5SxNQmoUS2n2wRE5B7gZWCWiNSKyFXGmBRwDfA4sBa43xjz3tCEmlndVUB1nXVZjkQNxGhtE7jzzju55pprsh2GyhH7uojmYJuAMebSfpY/Cjw66IiyJGAHKA4Wa0lghAjaNklspJ+byig1KvSUBEZOm8CIFo1EtSQwQnSXBOwcLQl0dHRw/vnnc9xxxzFv3jx+9atf8clPfrLn9WeeeYYLLvB6Tf/yl7/kmGOOYenSpbz44ovZClnloJ42gVwsCYxGFeEKnTpihOjuHSTmMCWBx74Ou98Z2o2Pnw8fOdRAePjLX/7CxIkTeeSRRwBoaWnh29/+Nh0dHeTl5XHfffexfPlydu3axXXXXceqVasoKirizDPPZOHChUMbrxqxcrpNYDSKRqLUdWlJYCSwLSElfnwmCcZkO5yDzJ8/nyeeeIKvfe1rPP/88xQVFbFs2TL+/Oc/k0qleOSRR7jooot49dVXOeOMM4hGowQCgYOmmVBjmwxDdZCWBHqJhqM0dOmo4ZHCSPrf13XA7udf+TBX7JlyzDHH8MYbb/Doo4/yrW99i7POOovly5dz8803U1paypIlSygoKMhKbGoEGYaGYS0J9FIRqcAxjo4aHiFc2/uC5OLUETt37iQSifDZz36Wa6+9ljfeeIMPfvCDvPHGG/zv//4vy5cvB+DEE0/k2WefpaGhgWQyye9+97ssR65yiYggfv/oHTGca3ruNayjhkcG8Xnz0OZg4/A777zDtddei2VZ+P1+br31Vmzb5oILLuDOO+/kV7/6FQATJkzg+uuv5+STT6a4uJjjjz8+u4GrnKNJYBhVhL3b/9V31TOHOVmORh2OawfAJSenjjj33HM599xzD1p+8803c/PN+0+TdeWVV3LllVcOV2hqhBG/f0RNGzGidZcEtJvoCGGlr2FyMAkoNVQkENDeQcOlLFyGINpNdIQwlteHOhfbBJQaKpmuDtIk0Ivf8lMSKtFuoiOEdPcI0lHDahTTJDDMKiI6YGykkJ7eQVodpEYvCWibwLCKhnXqiBHDp9VBavQTv7YJDKuKSIXeU2CEECtdEsjBLqJKDRWtDhpm0Yg3ajil9cw5T3pKArmXBLZs2cK8efOyHYYaBTQJDLNoOIrB0NDVkO1Q1GFYOZwEDieV0osMNTCZbhPQwWIHqIjsGzA2Lm9clqNRhyIjpE2gpqaGT3ziE1x22WU8//zztLe34zgOU6ZM4eMf/zgXX3wxAJ/5zGf41Kc+xUUXXZTdgFVO8doEWjO2fk0CB9ABYyOH5etuE+j/qvqHr/2QdY3rhnS7s0tn87WlXxvQe99//32WL1/OnXfeyerVq3njjTd4++23KS0t5dlnn+W///u/ufjii2lpaeGll17qmU5CqW4SGMXVQSLyARH5mYj8QkReymYs3XqmjtBuojnPsnO7JFBfX89FF13Eb37zG4477jgAzjnnHEpLSwH44Ac/yIYNG6ivr+eee+7hE5/4BD6fXpep/WV62ohB/8eJyB3ABUCdMWZer+XLgJsAG/iFMabfuXyNMc8Dz4vIxXg3qc+60lApllg6YGwEsAfQJjDQK/ZMKCoqoqqqihdeeIG5c+cCkJeXt997rrjiCu6++27uvfdefvnLX2YjTJXjcnkCuTuBm4Ffdy8QERu4BTgHqAVWishDeAnhxgM+/3ljTPeZ9jLgqqOIZcjYlk15qFxLAiOA5U/fes9JIFmOpS+BQIA//OEPnHvuueTn5/f5nr/7u79j6dKljB8/vidRKNVbpucOOpobzT8nItUHLF4KbDTG1ACIyL3ARcaYG/FKDQcRkSqgxRjT1s/rK4AVAFVVVYMN94joHcZGhu6SgJNK5mzjVl5eHg8//DDnnHMOl19++UGvjxs3jjlz5vQ0Dit1oFwuCfSlEtje63ktcOJhPnMV0G852BhzG3AbwJIlS4blPoLRSJRd7buGY1PqKNj+IABOMp5zSaC6upp3330XgOLiYlau7Lu2s7Ozkw0bNnDppZcOZ3hqBBn1U0kbY64zxuREo3C3irCOGh4J7PSt95wM3n81k5588knmzJnDF7/4RYqKirIdjspRI60ksAOY3Ov5pPSyEaU8Uk5jrJGkk8TfPUmZyjm+dEkgNUKTwNlnn83WrVuzHYbKcd1tAsYYRIa+9WuoSwIrgZkiMlVEAsBy4KEh3kbGdXcT3du1N8uRqEPx+bvbBEZmElBqIMTvB2PAcTKy/kEnARG5B3gZmCUitSJylTEmBVwDPA6sBe43xrw3NKEOn54BY9o4nNP8gRAAbiqe5UiUyhxJV3tmql3gaHoH9dmSZYx5FHh00BHlgJ6pI7SbaE7zp0cMO6mRN3eQUgMlgXRX6Ay1C2S9YTgXRcM6dcRIEPTbxI0Pk9SSgBq9ekoCmgSGT0moBJ/4NAnkuIDPIoWNm4OziPYeHHbttddy7LHHcu211x7yfUr1JdNJINe6V+cESyzKwmXaMJzjAj6LJD5MjjcM33bbbTQ2NmLbdrZDUSNQptsEtCTQj/JwOXtjmgRyWcC2SOZoSaDbhRdeSHt7O4sXL+a+++5j8+bNnHzyycyfP59vfetb2Q5PjQCZbhPQkkA/ouEouzp01HAu6ykJHCIJ7P7+94mvHdqppINzZjP+//2/Ab33oYceIj8/nzfffBPwksI//MM/cMUVV3DLLbcMaVxqdNI2gSwpC5fpqOEcF/BZpIwNOV4d1NuLL77YM0VEX3MJKXWgnO0iOtpFI1GaYk2k3BQ+S/9MuShgWyTw4T/EjeYHesU+nDIx6lONXloSyJLyUDkGQ2OsMduhqH4E09VBI+kew6eeeir33nsvAL/5zW+yHI0aCXScQJaUR8oBnToil3V3ER1JSeCmm27illtuYf78+ezYMeKm1VJZoF1Es6R7wJgmgdwV9NnpkkDutQm0t7f3+fvUqVN5+eWXe55/97vfHda41MijXUSzpDysJYFc5/UOspFD3GheqZFO2wSypDsJ6PxBucu2hBQ+xM29koBSQ0XbBLIkYAcoDBRqSSDHOeLTkoAa1bQkkEXRcFSTQI5zxIfVRxdRY4blTqRZNRb2UWmbQFaVh8t1wFiOc/soCYRCIRoaGkb1SdIYQ0NDA6FQKNuhqAzTaSOyqDxSzpt1b2Y7DHUIjvixzP5fjkmTJlFbW0t9/ehO4KFQiEmTJmU7DJVh2kU0i7qrgzJ1b0919IzlwzqgJOD3+5k6dWqWIlJqaI3q6iARmSsi94vIrSJySTZj6Ut5uJy4E6c92X74N6uscC0/thk5g8WUOmI+71o95xqGReQOEakTkXcPWL5MRN4XkY0i8vXDrOYjwE+NMf8AXDHYWDKlp5uotgvkLGP5sU1mbsCtVC4QESQQyL0kANwJLOu9QERs4Ba8k/tc4NL01f58EXn4gEcFcBewXET+Eyg7ilgyojsJNHQ1ZDkS1R9XtCSgRj/x+zGJHGsTMMY8JyLVByxeCmw0xtQAiMi9wEXGmBuBC/pZ1T+lk8fvBxtLpnRPHaEDxnKY7cM2Ok5AjW7i92OSI2Mq6Upge6/ntcCJ/b05nUT+H5AH/Gc/71kBrACoqqoaqjgHpCzsFU50rEDuMpYfH5oE1OjmJYEcKwkMBWPMFtIn+EO85zbgNoAlS5YMa8fvwkAhASugSSCX2X58OGAMaA8uNUpV33cvViSSkXUPdRLYAUzu9XxSetmIJCI6YCzHGdsbSIOTBF8gu8EolSH+iRMztu6h7iK6EpgpIlNFJAAsBx4a4m0Mq/JIuZYEcphYXh9qDnF3MaVU/46mi+g9wMvALBGpFZGrjDEp4BrgcWAtcL8x5r2hCTU7dP6gHNe7JKCUOmJH0zvo0n6WPwo8OuiIckx5uJxVe1ZlOwzVD/GlSwKaBJQaFJ1A7jDKw+U0x5tJ6kkmJ4mt1UFKHQ1NAofRM2AspgPGcpH0VAfpjWWUGgxNAoehA8Zym5WuDkplaCCNUqOdJoHD0HsN5zZJdwtNJuJZjkSpkUmTwGHoJHK5zbK7SwKaBJQaDE0Ch1EaLkUQLQnkKNsfBLQkoNRgaRI4DL/lpyRUoiWBHNXdJpDUNgGlBkWTwACUh3XUcK7qLgmktCSg1KBoEhiA8nC53lMgR/m0d5BSR0WTwADoJHK5yw54JQEnpUlAqcHQJDAA3dVBxgzrTNZqAHy+dBLQkoBSg6JJYACi4SgpN0VLvCXboagDdLcJONpFVKlB0SQwADpgLHf5At5gMVerg5QaFE0CA9B9m0ltF8g9Pr+2CSh1NDQJDED3/EFaEsg9fr+WBJQ6GpoEBkCrg3JXIN07SJOAUoOjSWAA8vx5hH1hTQI5qCcJ6P0elBoUTQIDICKUhcq0TSAH+YNeEjBaElBqUIYtCYjINBG5XUQeONSyXBWNRHXUcA7qLgmYlJYElBqMASUBEblDROpE5N0Dli8TkfdFZKOIfP1Q6zDG1Bhjrjrcslylo4ZzU8DvxzGC0eogpQZloCWBO4FlvReIiA3cAnwEmAtcKiJzRWS+iDx8wKNiSKPOAp1ELjf5bIsUPr29pFKD5BvIm4wxz4lI9QGLlwIbjTE1ACJyL3CRMeZG4IKhClBEVgArAKqqqoZqtUesPFxOW6KNWCpGyBfKWhzqYEl8WhJQapCOpk2gEtje63ltelmfRKRMRH4GLBSRb/S37EDGmNuMMUuMMUui0ehRhHt0uscK6A3nc09KbMTVJKDUYAyoJDAUjDENwNWHW5arekYNd9ZTmd9vrlNZkNKSgFKDdjQlgR3A5F7PJ6WXjUo9JQHtIZRzUvgQTQJKDcrRJIGVwEwRmSoiAWA58NDQhJV79IbzucsRn1YHKTVIA+0ieg/wMjBLRGpF5CpjTAq4BngcWAvcb4x5L3OhZldpqBRLLO0hlIM0CSg1eAPtHXRpP8sfBR4d0ohylG3ZlARLNAnkIFd8iJvKdhhKjUg6bcQRiEaimgRykCN+LKMlAaUGQ5PAESgL6/xBucgVH5ZWByk1KJoEjkA0rCWBXORYPiytDlJqUDQJHIHycDmNXY24xs12KKoXY/mxjCYBpQZDk8ARKA+XkzIpmuPN2Q5F9eKKH1vbBJQaFE0CR6BnrECntgvkEmP5sLUkoNSgaBI4AjpqODcZ269JQKlB0iRwBHTUcG4ylh8bTQJKDYYmgSOgN5zPTcby49eSgFKDokngCET8ESK+iCaBXGOPrpKA6xruf307G+vash2KGgOGbSrp0aI8XE5dZ122w1C92QH8pHBcg21JtqM5arf8bSP/3xPrsQQuXljJP591DFVlkWyHpQbAcQ3NnQnK8oMHvea6hi0NHbyzo4VtDZ1cvLCSyaX7H9euhMOqrU28XLOXt7a3cO688Vx+0pSMxqxJ4AhNLZrKpuZN2Q5D9SK2Dx8OiZRLOGBnO5wj0hZLUhDy9zx/bn09P3pyPecvmMDEohC/fnkrD725k8+dUs2XzzmG/OC+r+y2hk6KIn6Kwv6+Vj0oXQmHV2oaOG1mOX77yCoKXNfw1zV7GF8U4vjJxUcVR2NHgtauJNXleUe1niO1fk8b1/3pPRKOy5WnVrPs2PH4DvN32Ly3g2ffr+PFTQ28UtNAWyzF1PI8TptRzrzKQjbWtfPOjhbe29FKW3xfifXWZzfx9Y/M5rMnTqGxM8Evnt/MXS9voSPhYFvChKIQ3/7ju2yqa+fbF8zN2AWOJoEjNKdsDs/veJ7OZCcRv16d5QQ7gD+dBLY3dfLQmztZ8cFpFIaO7uQYSzq4xhAJDO5r4roGq58v7t72ONc99B6PvL2LU6aX8a/LZlOeH+BL967mmIoC/vOSBUQCPq46bRo/fnI9t7+wmUfe3sW3LphD0nG5+5VtrNraRFlegB9+YgFnzx034Lg21bfz82c3UVPfwY0fn8/McQUAtHQmufLO13hjWzPTo3l864K5nDnLuz34juYu3t3Rgt8WCkJ+CkN+pkfzek6QO5u7+Mr9b/Fyjddzbml1KStOn0a0IMh7O1tZu6sV2xJmjsvnmHEFtMWSvFLTyCs1DSQdw+nHlHPGMRVYAr99bRuPvbObhOPyycWT+MZ5cyjNCxy0H9saOlm/p42Tp5eR1ys5uq5h7e5WArZFSV6Akoj32aTjknRc/LZF0Gchsu/YpByXnz9Xw01PbiA/5KMw5OOa366msjjMsnnjKUkn26DfJuUYko7LloYO/rauji0NnQBMLg1z/vwJTCnLY+WWRh58o5a7XnEI+CzmTCjkooUTWTCpmAWTioj4fXzzj+/wb396j3te286WvR3EUg4XLJjIxxdVckJ1KWG/zfceWcsdL25mZ3MXNy1fmJGLHDHGDPlKM2XJkiXm9ddfz2oMT297mi/97Uvc9ZG7OL7i+KzGojzv/PorzNl0B9cvfI77Vm4n6RjOPXYcP/vs4v2+6Efi8fd28/UH36Y1lmLWuAIWVhUza3wBFQVBogUhKgqClOUH+kwQ7+9u4wePreXZ9fUURwKU5wcYVxhizoRCjp1YSCzp8IPH1tERd/jYwkqeWLuHxo4EpXkBkimXh754GlMPuAJ+Y1sT3/zDu6zd1QpAdVmETy6ZzMNv72LtrlYuO7GKFR+YRtJxiadcmjoT7GqJsas5RmciRdBvE/JbvL29hcfX7CZgW4QDNomUy40fn8/J08q44o7XqKnv4JoPzeCPq3dQs7eD4yYVUd8WZ2dL7KD9LAj5OG1GOXMmFPKL52tIuYZvnj+HWNLljhc2s6O5a997gz4cY+hMOD3LArbF8VXF2CKs3NJIyjU96/3EokkEfBZ3vLCZ/JCP//OBaUQCNrGky57WGM+tr6dmb0fPuj+xeBIfmTeeFzfu5cE3duy37b6IQMhn47MF2xIc19AWS/GReeO54eJ5lEYCPLl2D7e/sJm3apuJJQ+eJSDgszh5Whkfml3BmbMqDqqyS6Rcaps6mVwa6bNUZYzhd6tq+Z+/bWRRVQn/9KEZTI/mH/S+X764mX9/eA0XLJjITy9deMj96n9/ZZUxZkmfr2kSODK7O3ZzzgPn8I2l3+CyOZdlNRblefe332De+v9hauxuPrZoMpNKIvzkqQ18+4K5XHXaVMA7Mf/4yfWce+x4Ljp+Yr/JoTOR4oaH13LPa9uYV1nImbMqeHN7M29ua96vKN8t7LeZUBxi9vgC5owvZFtjJw++UUt+0McliyeTcBzq2+LsaO5i/e52Eo53MllYVcx/fGIBM8cV0B5Pcfvzm7lv5Ta+c9E8zunnqj7luDzyzi5K8wKcOr0cyxLiKYcf/XU9tz1fQ39f5aDPIp7ytlsY8nHFydX83anVpBzDF+95g5VbmiiO+EmkXP73iiWcOqOcRMrlrle28sCqWqZH8zihupQFk4oQEdpjKfa2x3mlpoFn19ezqyXG8ZOL+fGnj++pvkk6Lk+trcMYw7ETi5hcGsYY2NnSxYY97QT9FouqSgj5vSvbtliSFzc2EE85fHju+J4r3vV72vjWH97ltS2N++3PSdPK+OAxUaZF8/jj6h08mi45iMBpM8q58LiJBP02TR0JmjoTCILfJ/gti6TrEks4dCUdko7BNQbHNXxgZjnL5k3o828YTzm0dCWJJ72ShM8W8oO+nvgz7ck1e5gWzWNaH0liIDQJDCFjDGfcfwanTzqdG069IauxKE/TX26k5JUf8NKlazhlViXGGFbctYq/ravjvr8/iXdqW/j+Y+twXUPKNZxQXcL1Fx7LsROL9lvPzuYu/u6Xr7Ghrp0Vp0/jK+fMIuDzruBc17C3I05da5z6tjj17XEaOxI0tMfZ1tjJut1tbG3oJGBbXHHyFK750AyKI/tXYSQdl0317extS3Dy9LIhreN9u7aZ9XvaCfq8qo6isJ+JxWEqCoMEfTbGGOIpF9uS/a5Kk47Lfz7+Po++s4ufXrqQhVUlR7RdYwy7WmKMKwxlrM7aGEN9e5yAbRH02QR91kHVbHvb47y0qYETqkuYUBTOSBwjmSaBIXb1E1ezt2svD1z4QLZDUQAv/gSe+DZ8YwcEvSullq4kF/z0eXa3xEg6hjNnRfnhJQt4em0d//H4+zR3JvjMiVP4yoePoTgSYGNdO1fc/iptsRT/89lFfGBm9IjD6IinSDmGosjQNdSOBHEnTtA+uDeMGjrGGFJuCr89uP+tQyWBYWsYFpFpwDeBImPMJellc4AvAeXAU8aYW4crnqMxt2wuv3z3l/rPnyu6vxhOomdRUdjPrZ9ZzD/f9yaXnzSFK06egoiwfGkVH5k3gR898T53vbKVh9/eyedPncodL27GtoR7VpzEvMqifjZ0aL0bJ0eLWCqGJRZ+y4/BsLtjNzUtNWxq3sR7e9/jnb3vUNtey5mTz+Q7p3yHklDfJYmEk6A92Y5rXBzXwW/7KQwU4rN8tCfaWV23mpV7VtKeaOe46HEsqljEpIJJg27TAWiONfPEtid4etvTtCZacV0Xg+GsqrO44tgr9vvubmrexPa27cSdOHEnTnGwmDmlc4hG9l0MGGMwGCzZv37fGENropXGWCNNsSaaYk00xBpo6GqgK9XF9OLpzCmbw7Siafisff8jLfEWtrRuYWf7TqoLqzmm5Bhsy95vvRubN/L4lsd5fMvjfKjqQ3x58ZcH/ffoz4D+a0XkDuACoM4YM6/X8mXATYAN/MIY84P+1mGMqQGuEpEHei1bC1wtIhbwa2BEJIE5ZXNImRQbmjYwr3ze4T+gMqs7CRxwT4F5lUU8+S8fPOjtRRE/37loHsuXVnHdQ+/x/z2xnkklYe6+6sRh75KYq9oT7fz4jR9z//v3YzAIgi02qV4jsyfkTWBe+TxOn3Q696+/n0seuoTvf+D7HFt2LGsb1/Le3vdY17SO9xvfZ3PLZhzjHLSdPH8esVQMxzj4LB9hO8zv1v8OgIgvQn4gn4gvQmGwkOlF05lRPIPJBZMxGJJuks5kJzUtNWxs3khtWy1BO0hBoAARYXXdalJuiurCaibkTcC2bDqTnfxk9U94cMODfGXJV4ilYtz3/n28Vf9Wn3+H8nA5xcFimuPNNMeacXEpDZVSHi4nYAXY27WX+q56kv3c1Mhn+Uil/y994iNgB3pO9G2J/QcDFvgLOK7iOCyxaIo1Ud9Vz+6O3VhicUJ4InMldOQHcgAGeulyJ3Az3okaABGxgVuAc4BaYKWIPISXEG484POfN8b0OcJKRC4E/gG464giz6I5pXMAWNOwRpNALrAOLgkMxJwJhdy34iRe3NjA7AkFlPce4LNzNTz2dSicAOf/CCKlQxhwduzp2ENtey0LKxYedDXbzRjDs7XPcsMrN1DfWc/HZ36cyvxK4k6clJuisqCSaUXTmFY0bb+r/otnXMy/PvevfOGvX9hvfeMi45hVOoszJ59JebgcW2wsyyLhJGhNtNIabyXPn8eS8Us4LnocQTtITXMNb9S9weaWzXSmOulMdtIQa+DZ2mf5w8Y/HBRzwAowtWgqM0tmknSTtCXa6Ep28dk5n+X8aeczq2TWfiWKV3a9wg9f+yH/8sy/ADClcApfXfJVFo9bTNAOErJD1HXVsa5xHWsa1vSUTkpCJQhCY6yRvV17iTtxFhcupjxSTnmonNJwKSXBEkpCJZSFyigNl2JhsbV1K2sa11DTXNPzd3SNS2V+JdVFXoLa0LyB13e/zlv1b+GzfJQES6gqrOL46PGc3d5B+cP/AoHpsGgo/hP2N+A2ARGpBh7uLgmIyMnA9caYc9PPvwFgjDkwARy4nge6q4MOWP6IMeb8PpavAFYAVFVVLd66deuA4s0kYwyn3XsaH67+MNedfF22w8kYY8xRFceHzVv3wR9WwIU3w6LLj25d8Xb42/fg1Z9BuBRiLZBXDh/7GUw746C3p9wUKTdFyJe+SnNdTFcTr2x7mpfrVyORMnx2gMJAIedWn8v4vPEDCsM17n4n6s54OyvX3MPKmseYnTeZ8xZdjRWd7fV1xDtWwMHHy3XAuDy+/Wmuf+l62pPtTC+azufnXclHpnwYbK86ZmvrVp7e+iR/2/oEWzp2MrN4Jtefcj0LogsODq5hE2x53vt7lFTvF+OvX/4uEirm2MmnMbdsLqWhI0yergsNG6B2JaTisODTPe08JLtoePo77Nz8NPas8/HPv4RQsJAJeRP2q2bpc53tu6FpKzRvAxFS/jBPtm6kODqXpVPP6TcpDkh7HbTthpIpEBpAVeKut+DlW6DqZG//AocYb1TzDNx9CVSdBJ99EHyDq34ekobhPpLAJcAyY8wX0s8vB040xlzTz+fLgO/hlRx+YYy5UUTOAD4OBIG3jTG3HCqGXGkYBvjCX79AW6KN+y64L9uhDDljDD9d/VMe3PAgN5x6A6dPOn1Yt7+peRPvN75PUbCI4mAxE/Mn9lvXDEBXE+7dn8DasQrmfwrO+09MqIjGWCO7O3ZT4S+gvHkHsuddHCfB5lQba+INbIjVszG2l83xRoLGZUIyyfjOVmbGu1g4bRnHnPMD7JZa6n5/FZvba/HNOJvZCz9P/tQzaUt18OD6B7l77d00dDVwcmg857S1EmnczC+L8nkvGMRnDAKkEIyAAKcQ4dyUjzafj522RaNlMbVgMvPGLaZ6whJW1b3BUzue56XmdVgIpVaAAiPUOO0kRbCNwRFhVjzBlxMB8vLH8ahp56/SRQeGSeKjSoJMcYWpnW1Ut+3lobwwDxTksUDCXByYwG9j29goKXzGkOqVNHzGsLQrxtmdnVxMIf4Fn4a5F3kn49Yd0LAR1j4Me97xPmD54YQvwOlfha0vwTM3Qt0aQGD2+XDKF2HcsdDVBJ2N3gm4YQPs3QjteyDeBol2b/3gJbT2eoi37Du2kTI49UtQcSw8di001kB0DtSvhaLJcPq1MOsjkO8NaqN5O7x+O7zzgJfAnYT3ONTdAMfNg6mnQ14UYs3e5+LtkOyERAf4w1A8xUt4wXxvXzoboGkz7FgNrbW94i2H6GyYfgbMOAfGLwArnWBcF16+GZ76dxALnDiES2Dx38GMs704wsX71rXnPbhjGRRNgisf2/+1I5QTSWAo5FIS+NHrP+LutXfz6mWvDrrF/lBc4/Jc7XP8dctfOS56HOdPO5/8wOD6CB+JhJPg2y9+m0c3P0ppqJTmeDNfWfwVLp97OSJCU6yJdY3rmFM6h+JQ8YDXW9tWyzPbn+GFHS9QGCzk7KqzOa3ytP1GXde01HDrm7fy+JbHMez7vxSEhRULOXvK2Zw04SQKAgWEfWGa4808te0pntr6FG/vfZt88VOS6CKCxU6fTZvsW0ee6zIhlWKHz0dX+ksZcA1Tk0mmJZMk7QC7QhF2WUKj8aqVIj4vts5U5377MsWBvbZFBy5L4ymOiXXyt0iEHX7varTKX8SVEz7AhSXzCex6C2pXUtu8iT8WFPLHkMUe8U5IEQMljstOWzC9TsYTUinO6OgigKHR9tHkDzI9Mo7Tpi7j+OOu5Oktf+Enb/+cHUlv4FjQCKdLmPHGZrtJsI0k28Ul2atQcJVdwT+1tOGPt2NKp/JcYQmvS4I816UglSKKxUll8ymMzgHL551EN/wVDqzHn7QUjr0YppwCr98Bq+/2lhsXymbAB77iJYuVt3sn1L4UTPAewQLvYQd6SjSEiqByCUw6wfv8MzfCpqe910qq4aM3eSWQTX+Dp77jVdsBlEyF4smw5QXv+cwPQ+k0r73IDkDBeO/zxel5eBLt3sl+xyqoeRa2vwqpmPfeUBEE8iGQB/6Ilwyatnif6eYLQeFEmLgQKhdDYSU0b/X2fddb3gO80mRJtXcib6+D7a/A7Avgwp9C/Tp45X9g3SP7klThJLBsL0HGWiB/HHzhCe/zRyFTSWBQ1UFHI5eSwGObH+Nfn/tXfvfR3zG7dPZ+r72882V2deziYzM+dkTVKSk3RU1LDa/teo171t3DtrZtRHwROlOdhH1hzpt6HovHLWZywWQq8ytpijdR01zD5pbNjM8bz6mVp1IRqTjk+t9vfJ836t7ANS7nTT1vv94PO9p38O0Xv83K3Sv50qIvcdnsy/jmC9/kyW1PcsakM2iMNfLO3nd6ekjML5/PKRNPoTRUit/yY4nFro5dbG3dyva27XQmO0mZFF2prp5J96YWTaU51kxTvImgHWRyweSeovjG5o0E7SCfnfNZlk1dRkeyg5Z4C2sb1vLktidZ37S+z/2aWzaXkyacRNyJ09i4iY76NUzEppoA460Qe/LL2BoKsxOHyvxKji2eydzCqUwpnILPF/S+dKHinhPRrvZdrK5bzVv1b2EwTCuaRnVRNYlYK2vXP8Ta3avIS8b4THgyc4uPgYrZmJnLWJNqoSnexMkTTt6vl0dvjuuwpXULZaEyioLe4KuOllrWbP4rm+rfYX7RdOZWLEQKxntXwcHCfVeSvSScBI/UPILP8nHm5DMPukBwXIed7TupaamhPFzOseXH9v+P15/2Otj8nHcFWljpnfQOrO6of9874U9cCPM/CXa6WibRAe8+CF3N3tVuuASKKr1EESw4sji2veKVMBYs37/qxBiv2mjbK1D7GtSvh1nLYMlVXtXMkUilSwv+fhpfjfFKAIl2r3rQH9mXuPrSXgcbn4JtL3mlk9Yd3t/kg1+DRVfs/9mOvbDrTdj9DuxZ470WLPD+1sd/BsqmH9m+9CFTScAHrAfOAnYAK4HLjDHvHXXE/cilJLC1dSsX/OEC/v2Uf+djMz8GePcZ+I+V/8Fjmx8DvMayfzv53/BbfhJOglvf8q5yJ+ZNpLqomnGRcTTGGqnrrGNn+042NG8g7nhF4wXlC7h87uWcNeUs1jWs43frf8djmx8j5hw8fL+3OaVzmFQwiaSTJGmSJJwE8VScLqeL2rZaulL7htPbYnP6pNOZVTqLZ7c/y9rGtfgsH9899bucP81rnnGNy82rb+ZX7/2K2aWzOa3yNOaVz+Odve/wfO3zvNvw7n7bF4SJ+ROZXDCZgkABPsuH3/Izq8RrHJxcOJmUm2J13Wqe2vYUuzt2Y4zBxWVG8Qwun3t5v/XI21q38c7ed4ilYnSluvBbfk6bdBqV+ZWDO4hKjRFHnQRE5B7gDLz+/HuA64wxt4vIecCP8XoE3WGM+d5QBd2XXEoCrnE55Z5TOL7ieE6ecDI723fy55o/E0vF+D/z/w+Ocfj52z/n1Imn8oX5X+B7r36Pjc0bOWXiKbQn29ncspm2RBt5/jyi4Sjj88Yzq2QWc8rmMLdsLlOLph60zYSToLa9lto271ESKmFa0TSmFE5ha+tWXtjxAi/seIHmeDM+y9fTJS3kCxGyQ1REKlg0bhELKxYSS8X4/cbf86eNf6Ix1shx0eM4q+oszplyDpMKDi56HthQ2a0z2UlXqoukmyTlpohGojp2QqkcoyOGM+Tvn/h7Xtr5EuDVHy8ct5CvnfC1nhP4g+sf5IZXbsAxDtFwlOtPub6nkdUYQ9yJ7+tVkiXdfa2LgoMbIKWUyn2aBDKkJd7C7o7dTMifQIG/oM/6/5d2vMTLu17mC/O/oCdapVRW5MS0EaNRUbDosCf2UypP4ZTKU4YpIqWUOjJ6j2GllBrDNAkopdQYpklAKaXGME0CSik1hmkSUEqpMUyTgFJKjWGaBJRSagzTJKCUUmPYiBoxLCL1wJHcVaYc2JuhcHLZWNzvsbjPMDb3eyzuMxzdfk8xxkT7emFEJYEjJSKv9zdUejQbi/s9FvcZxuZ+j8V9hsztt1YHKaXUGKZJQCmlxrDRngRuy3YAWTIW93ss7jOMzf0ei/sMGdrvUd0moJRS6tBGe0lAKaXUIWgSUEqpMWzUJgERWSYi74vIRhH5erbjyQQRmSwifxORNSLynoh8Kb28VESeEJEN6Z8l2Y41E0TEFpHVIvJw+vlUEXk1fczvE5FAtmMcSiJSLCIPiMg6EVkrIiePhWMtIl9O/3+/KyL3iEhoNB5rEblDROpE5N1ey/o8vuL5SXr/3xaRRYPd7qhMAiJiA7cAHwHmApeKyNzsRpURKeArxpi5wEnAP6X38+vAU8aYmcBT6eej0ZeAtb2e/xD4b2PMDKAJuCorUWXOTcBfjDGzgePw9n1UH2sRqQT+L7DEGDMPsIHljM5jfSew7IBl/R3fjwAz048VwK2D3eioTALAUmCjMabGGJMA7gUuynJMQ84Ys8sY80b69za8k0Il3r7+Kv22XwEXZyXADBKRScD5wC/SzwX4EPBA+i2jar9FpAg4HbgdwBiTMMY0MwaONd5tcMMi4gMiwC5G4bE2xjwHNB6wuL/jexHwa+N5BSgWkQmD2e5oTQKVwPZez2vTy0YtEakGFgKvAuOMMbvSL+0GxmUrrgz6MfCvgJt+XgY0G2NS6eej7ZhPBeqBX6arwH4hInmM8mNtjNkB/BewDe/k3wKsYnQf6976O75Ddo4brUlgTBGRfOBB4J+NMa29XzNeH+BR1Q9YRC4A6owxq7IdyzDyAYuAW40xC4EODqj6GaXHugTvqncqMBHI4+AqkzEhU8d3tCaBHcDkXs8npZeNOiLix0sAvzHG/D69eE930TD9sy5b8WXIqcCFIrIFr6rvQ3j15cXpKgMYfce8Fqg1xryafv4AXlIY7cf6bGCzMabeGJMEfo93/Efzse6tv+M7ZOe40ZoEVgIz0z0IAngNSQ9lOaYhl64Hvx1Ya4z5Ua+XHgI+l/79c8Cfhju2TDLGfMMYM8kYU413bJ82xnwG+BtwSfpto2q/jTG7ge0iMiu96CxgDaP8WONVA50kIpH0/3v3fo/aY32A/o7vQ8AV6V5CJwEtvaqNjowxZlQ+gPOA9cAm4JvZjidD+3gaXvHwbeDN9OM8vPrxp4ANwJNAabZjzeDf4Azg4fTv04DXgI3A74BgtuMb4n09Hng9fbz/CJSMhWMNfAdYB7wL3AUER+OxBu7Ba/dI4pX8rurv+AKC1wNyE/AOXu+pQW1Xp41QSqkxbLRWBymllBoATQJKKTWGaRJQSqkxTJOAUkqNYZoElFJqDNMkoJRSY5gmAaWUGsP+f+JbVlSP1PLxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(r_list, BE_eig_r, label=\"eig-r\")\n",
    "plt.semilogy(r_list, BE_svd, label=\"svd\")\n",
    "plt.semilogy(r_list, BE_kry, label=\"kry\")\n",
    "plt.semilogy(r_list, BE_lfd, label=\"lfd\")\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\textbf{Case 2: Image-based raw features}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampling_states(P, size, init_state = None):\n",
    "    \"\"\"\n",
    "    init_states\n",
    "    return list of sampled states\n",
    "    states in {0, ... , n}\n",
    "    \"\"\"\n",
    "\n",
    "    states_list = []\n",
    "    if (init_state == None):\n",
    "#        import pdb; pdb.set_trace()\n",
    "        init_state = np.random.choice(states)\n",
    "\n",
    "    cur_state = init_state\n",
    "\n",
    "    for i in range(size):\n",
    "        p = P[cur_state, :]\n",
    "        next_state = np.random.choice(states, size=1, p=p)[0]\n",
    "        states_list.append(next_state)\n",
    "        cur_state = next_state\n",
    "\n",
    "\n",
    "    return states_list\n",
    "\n",
    "\n",
    "\n",
    "def check_transition_matrix(P, transition_matrix_rank):\n",
    "    epsilon = 1e-6\n",
    "    r = transition_matrix_rank\n",
    "    assert(np.linalg.norm(P.sum(axis=1) - np.ones(P.shape[0])) < epsilon)\n",
    "\n",
    "    U, S, V = np.linalg.svd(P, full_matrices=False)\n",
    "    P = (U * S) @ V\n",
    "\n",
    "    assert(np.linalg.norm(S[r:] - np.zeros(P.shape[0] - r)) < epsilon)\n",
    "\n",
    "    return \n",
    "\n",
    "\n",
    "\n",
    "def get_transition_matrix(transition_matrix_rank, n):\n",
    "    tmp_matr1 = np.random.random((n, transition_matrix_rank))\n",
    "    tmp_matr2 = np.random.random((transition_matrix_rank, n))\n",
    "    P = tmp_matr1 @ tmp_matr2\n",
    "\n",
    "    sum_rows = P.sum(axis=1)\n",
    "    for j in range(n):\n",
    "        P[j, :]  = P[j, :] / sum_rows[j]\n",
    "\n",
    "    check_transition_matrix(P, transition_matrix_rank)\n",
    "    \n",
    "    return P\n",
    "\n",
    "\n",
    "\n",
    "def create_raw_features(n=100, l=50, transition_matrix_rank=40):\n",
    "    states = [i for i in range(n)]\n",
    "    P = get_transition_matrix(transition_matrix_rank, n)\n",
    "    A = np.zeros((n, l))\n",
    "    A_next = np.zeros((n, l))\n",
    "    state_to_vec = dict()\n",
    "    for state in states:\n",
    "        state_to_vec[state] = np.random.binomial(n=1, p=0.5, size=l)\n",
    "    sampled_states = sampling_states(P=P, size=n+1)\n",
    "#    print(sampled_states)\n",
    "    for i in range(len(sampled_states) - 1):\n",
    "        A[i, :] = state_to_vec[sampled_states[i]]\n",
    "        A_next[i, :] = state_to_vec[sampled_states[i + 1]]\n",
    "    \n",
    "    reward = np.random.random(A.shape[0])\n",
    "\n",
    "    return A, A_next, P, reward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A, A_next, P, reward = create_raw_features()\n",
    "#P_A = np.linalg.pinv(A) @ A_next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 50)\n",
      "(100, 50)\n",
      "> \u001b[0;32m<ipython-input-8-b6bebc85314e>\u001b[0m(161)\u001b[0;36mcalculate_BE\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m    159 \u001b[0;31m        \u001b[0;32melif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmethod\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"approx\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m    160 \u001b[0;31m            \u001b[0mF\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mF_hat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_approx_feature_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m--> 161 \u001b[0;31m            \u001b[0mP_F\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpinv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mF\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mP_A\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mF\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m    162 \u001b[0;31m\u001b[0;31m#            P_F = np.linalg.pinv(F_hat) @ self.P_A @ F_hat ## smth strange\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\u001b[0;32m    163 \u001b[0;31m\u001b[0;31m#            reward_F = np.linalg.pinv(F) @ self.reward_A ## smth strange\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0m\n",
      "ipdb> F.shape\n",
      "(100, 10)\n",
      "ipdb> F_hat.shape\n",
      "(50, 10)\n",
      "ipdb> A.shape\n",
      "(100, 50)\n",
      "ipdb> self.P_A.shape\n",
      "(50, 50)\n",
      "ipdb> F.shape\n",
      "(100, 10)\n"
     ]
    }
   ],
   "source": [
    "A, A_next, P, reward = create_raw_features()\n",
    "print(A.shape)\n",
    "print(A_next.shape)\n",
    "mdp = mdp_low_rank(A=A, A_next=A_next, P=P, reward=reward, method=\"svd\", r_start=10)\n",
    "%debug\n",
    "BE_svd, r_list = mdp.calculate_BE_list(method=\"approx\")\n",
    "print(\"svd\")\n",
    "#mdp = mdp_low_rank(A=A, A_next=A_next, P=P, reward=reward, method=\"kry\", r_start=10)\n",
    "#%debug\n",
    "#BE_svd, r_list = mdp.calculate_BE_list(method=\"approx\")\n",
    "#print(\"kry\")\n",
    "#mdp = mdp_low_rank(A=A, method=\"kry\")\n",
    "#BE_kry, r_list = mdp.calculate_BE_list()\n",
    "#print(\"kry\")\n",
    "#mdp = mdp_low_rank(A=A, method=\"svd\")\n",
    "#BE_svd, r_list = mdp.calculate_BE_list()\n",
    "#print(\"svd\")\n",
    "#mdp = mdp_low_rank(A=A, method=\"lfd\")\n",
    "#BE_lfd, r_list = mdp.calculate_BE_list()\n",
    "#print(\"lfd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fc3891a7f70>"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlK0lEQVR4nO3de3RdZZ3/8fc390vTJG1JL0lLS2mrpdgWOwXECgpKGZHOqGAZWDgFyujIjM6anyP+xlnoGpezhnF0ZkTrD1SKIwJSAemIXEYFFKHQC0JbWugF2vSStE1zz0lyku/vj3PSpuEkzTk57d7n5PNaiwVn7332/uasEz55nmfv5zF3R0REZKCcoAsQEZFwUkCIiEhCCggREUlIASEiIgkpIEREJKG8oAtIlwkTJvj06dODLkNEJKNs2LDhsLufkWhf1gTE9OnTWb9+fdBliIhkFDN7e7B96mISEZGEFBAiIpKQAkJERBLKmjEIEZF06e7upra2lkgkEnQpaVNUVERNTQ35+fnDfo8CQkRkgNraWsrKypg+fTpmFnQ5I+buHDlyhNraWmbMmDHs96mLSURkgEgkwvjx47MiHADMjPHjxyfdIlJAiIgkkC3h0CeVn0cBIRkl2tPL/S/toSvaG3QpIllPASEZ5aXdDXz54dd49JV9QZciEkqrV6/m1ltvTcu5FBCSUY62dwPwxOaDAVcikv0UEJJRmiOxgPj9m4dpif+3SDZqa2vjox/9KPPnz2fevHnce++9XH311cf2P/PMM1x55ZUA3HPPPcyePZvFixfz/PPPp60G3eYqGaWpIxYKXT29/GZbPcsWVAdckWS7r63dwtb9zWk959wpY7n9Y+cMecwTTzzBlClT+OUvfwlAU1MT//RP/0RbWxulpaU8+OCDLF++nAMHDnD77bezYcMGysvL+eAHP8jChQvTUqdaEJJRmju6yc0xqsoK+dVr6maS7HXuuefy9NNP86UvfYnf/e53lJeXs3TpUtauXUs0GuWXv/wly5YtY926dVxyySWcccYZFBQU8KlPfSptNagFIRmlOdLN2KI8Lj9nEg9t2Et7V5SSAn2N5dQ52V/6p8rs2bPZuHEjjz/+OF/5yle49NJLWb58OXfeeSfjxo1j0aJFlJWVndIa1IKQjNLUEaW8OJ8r5k0i0t3Ls9sPBV2SyCmxf/9+SkpKuP766/niF7/Ixo0bufjii9m4cSN33303y5cvB+D888/n2Wef5ciRI3R3d/PQQw+lrQb96SUZpbmjm7HF+SyeMY7Kknx+tfkgV5w7OeiyRNLutdde44tf/CI5OTnk5+ezatUqcnNzufLKK1m9ejX33nsvAJMnT+arX/0qF154IRUVFSxYsCBtNSggJKPEupjyycvN4SNzJ/HL1w7QGe2hMC836NJE0uryyy/n8ssvf8f2O++8kzvvvPOEbStWrGDFihVpr0FdTJJRmjq6KS+OzUa59NxJtHZG+f2bhwOuSiQ7KSAkozR3RBlbHGv4XjRzAmWFefx6W33AVYlkJwWEZJTmSGwMAqAgL4fZk8rYfagt4KokG7l70CWkVSo/jwJCMkaku4euaC9ji44veFJdUUxtY3uAVUk2Kioq4siRI1kTEn3rQRQVFSX1Pg1SS8Zojj9F3deCAKipLObx1w7Q0+vk5mTX9MwSnJqaGmprazl0KHtuo+5bUS4ZoQ4IMysFvgd0Ac+4+30BlyQB6puHqfyEgCgh2uvUNUeYUlEcVGmSZfLz85NaeS1bDauLycw+b2abzWyLmX0h1YuZ2Y/MrN7MNifYt9TMtpvZDjO7Lb7548Aad18JXJXqdSU79M3DNLbo+N811ZWxUKg92hFITSLZ7KQBYWbzgJXAYmA+cKWZnT3gmCozKxuw7YRj4lYDSxNcIxf4LnAFMBe41szmAjXA3vhhPSerVbJbc0cUeGcXE8A+jUOIpN1wWhDvBta5e7u7R4Fnif1l39/FwKNmVghgZiuB7ww8kbs/BzQkuMZiYIe773L3LuABYBlQSywkhlurZLFEXUzV8W6l2ga1IETSbTj/090MLDGz8WZWAvwpMLX/Ae7+EPAk8KCZXQfcCFz9jjMNrprjLQWIBUM18DDwCTNbBaxN9EYz+5iZ3dXU1JTE5SQTHe9iOh4QRfm5TBhTqC4mkVPgpIPU7v66mf0r8BTQBrxCgu4ed7/DzB4AVgEz3b11pMW5exsw5PPj7r4WWLto0aKVI72ehNvxu5hO/NrWVBazr1EBIZJuw+q2cfcfuvt73f0DwFHgjYHHmNkSYB7wCHB7knXs48RWSU18m8gxzZEoRfk575h3qaaymNqjGoMQSbfh3sVUFf/3NGLjDz8dsH8hcBexcYMVwHgz+3oSdbwMzDKzGWZWACwHHkvi/TIKNLV3n9C91KemsoT9jRF6e7PjoSaRsBjuwO/PzWwrsXGAz7l744D9JcA17r7T3XuBG4C3B57EzO4HXgDmmFmtmd0EEB/8vpXYOMbrwM/cfUsqP5Bkr/7TbPRXXVlMV08vh1o7A6hKJHsN60E5d19ykv3PD3jdDdyd4LhrhzjH48Djw6lHRqfmSPcJdzD1qTn2LEQ7E8cmN5WAiAxOt45Kxmjq6D7hIbk+U/WwnMgpoYCQjNEcX250oL4pNhQQIumlgJCMMdgYRElBHuNLCxQQImmmgJCM0NvrsfWoE9zFBLrVVeRUUEBIRmjritLrJOxigtidTHpYTiS9FBCSEZojfRP1Jb7xrqayhH1HO44t8HK0rYvLv/0cL+46ctpqFMk2CgjJCE3t75yHqb+aymI6o8efhfjZ+r1sr2vhj3sbT1eJIllHASEZIdFMrv31zeq672gHvb3Ofev2AHA0HiwikrxQrygn0ifRcqP91VSWALFbXZs6utnTEBuwPtrWdXoKFMlCCgjJCImm+u6v/8pyG95uYMKYQsYW59HQroAQSZW6mCQj9A1SD9bFNKYwj8qSfF7afYTfbKvnU39SQ1VZIY0KCJGUKSAkI/R1MY1JMNVGn+rKYn67/RAA1y6eRmVJAQ3qYhJJmQJCMkJTRzdlhXnk5tigx9RUxMYhPvSuKmoqS6gsLaBRg9QiKVNASEYYbJqN/vpmdb3ugjMBGFdSwNH2Lq0TIZIiDVJLRmjuiJ40IK6cP4XOaC8XzzoDgIqSfHodWiJRykuGfq+IvJMCQjJCc0c35YM8Rd1nwdQKFkytOPZ6XGkBAA3tXQoIkRSoi0kyQnNk8In6BlNZEg8IDVSLpEQBIRmhuePkYxADVcZbELrVVSQ1CgjJCE0diZcbHUplvFtJLQiR1CggJPSiPb20dfUk38V0rAWhW11FUqGAkNBrOclU34MpK8wjL8c03YZIihQQEnp98zAl28VkZlSUFGgMQiRFCggJvb6pvpPtYgIYV5qvMQiRFCkgJPSaO/q6mJIPiIqSAq0JIZIiBYSEXqpdTBCfbkMtCJGUKCAk9I51MSU5SA1QWZrPUY1BiKREASGh13ySxYKGUhnvYnLXhH0iyQp1QJhZqZnda2Z3m9l1QdcjwWjq6CYvxygpyE36vZUlBfT0+rEFh0Rk+IYVEGb2d2a2xcw2m9n9ZlaUysXM7EdmVm9mmxPsW2pm281sh5ndFt/8cWCNu68ErkrlmpL5+qb6Nht8LYjBaLoNkdSdNCDMrBr4W2CRu88DcoHlA46pMrOyAdvOTnC61cDSBNfIBb4LXAHMBa41s7lADbA3fljPyWqV7NQSiVI2xEpyQxlXquk2RFI13C6mPKDYzPKAEmD/gP0XA4+aWSGAma0EvjPwJO7+HNCQ4PyLgR3uvsvdu4AHgGVALbGQSKZWyTKtkShjClMLiIoSTbchkqqT/k/X3fcB3wT2AAeAJnd/asAxDwFPAg/GxwpuBK5Ooo5qjrcUIBYM1cDDwCfMbBWwNtEbzexjZnZXU1NTEpeTTNLSmXpAjNOU3yIpG04XUyWxv+ZnAFOAUjO7fuBx7n4HEAFWAVe5e+tIi3P3Nndf4e6fdff7BjlmrbvfUl5ePtLLSUi1dabexdS3JoRudRVJ3nC6bS4Ddrv7IXfvJvZX/fsGHmRmS4B5wCPA7UnWsQ+Y2u91TXybCK2dUUpTbEGUFeWRm2MKCJEUDCcg9gAXmFmJxW4juRR4vf8BZrYQuItYS2MFMN7Mvp5EHS8Ds8xshpkVEBsEfyyJ90sWG8kYRE6OUVGcT0ObxiBEkjWcMYh1wBpgI/Ba/D13DTisBLjG3Xe6ey9wA/D2wHOZ2f3AC8AcM6s1s5vi14gCtxIbx3gd+Jm7b0n5p5Ks0toZZUyKXUwQu9VVt7mKJG9Yv3XufjtDdBu5+/MDXncDdyc47tohzvE48Phw6pHRo7unl85oL2MKRhAQJZrRVSQVunVUQq2tM/YE9IhaECUFus1VJAUKCAm1vtXkUh2DABhXWqBV5URSoICQUGvtHHlAVMSn/NaEfSLJUUBIqKWji2lcaT7RXj8WNiIyPAoICbWW+P/UU30OAo5Pt3FUt7qKJEUBIaHWGh+DKBvJGISephZJiQJCQi0tdzH1zeg6jIB4+a0GvvPrNzVeIcIwn4MQCUprGrqYKkuGtyaEu/OPj7zGG3Wt5Ofl8JmLZ6Z8TZFsoBaEhFrfba6lI3pQrm9G16HHIF7YeYQ36lqpqSzmjie28eKuIylfUyQbKCAk1No6o5QW5JKbk/xqcn3GFueTYydvQaz+w1uMKy3gF5+7iOkTSrn1p5uob46kfF2RTKcuJgm1kczk2ic3x6goKRhyuo29De387+t1fObimYwfU8j3r38vy+58nlt/uonPXHLWiK4vcqrNmTSW6oritJ9XASGhNtKJ+vpUlOQPeRfTT158GzPj+gvOBGD2xDL+5ePn8oUHX+Gl1YkWQRQJj2/8+bn8xfnT0n5eBYSEWusIVpPrb9wQLYiOrh4eeHkvl58zkSn9/gr7s4XVLJxWoXmcJPSqK9PfegAFhITcSNaC6G9yRTGb9hxNuO8Xr+yjqaObT184/R37zhxfypnjR3x5kYykQWoJtXS1IOZMHEPt0Y6E023c+8LbvHvyWBbPGDfi64hkEwWEhFq6AmL2xDIA3qxrOWH70bYuXj/QzFXzpxBbMFFE+iggJNTSNUg9Z1IsIN4YEBBbDzQDMK967IivIZJtFBASWu5OW5paEFMrSyjKz+GNutYTtm/dHwuIc6aUj/gaItlGASGh1RntpbvHR/wcBEBOjjGrquwdLYgt+5uYXF7EuNKCEV9DJNsoICS0+gaUy9LQxQSxcYjtBwcGRDNzJ6t7SSQRBYSEVlsaVpPrb86kMdS3dHI0/jxEpLuHnYdaOWeKAkIkEQWEhNaxifrSFBB9dzL1dTNtO9hCr8NcjT+IJKSAkNDqa0GMZLGg/o7dyVQfG6jesr8JQC0IkUEoICS00rEWRH+TxhZRVpjHG/FxiC37mxlblEfNKZqmQCTTKSAktFrTsJpcf2bG7EllbI93MW3d38zcKWP1gJzIIBQQElqtae5igtg4xBt1LUR7etl2sJm5kzX+IDIYBYSEVmuaB6khNidTY3s3L+1uINLdq/EHkSEoICS0WjujmEFJQW7azjk7PlD9yKZ9AMxVQIgMSgEhodU3UV86xwj6bnX91eaDFOTlcHbVmLSdWyTbKCAktNK1FkR/E8YUMr60gNbOKHMmlpGfq18BkcHot0NCK11TfQ/U14rQFBsiQ1NASGila6rvgfoemDtHU3yLDEkBIaF1qlsQuoNJZGihXpPazEqB7wFdwDPufl/AJclp1NYZZdLYorSfd9mCKfS6s3BqZdrPLZJNTtqCMLM5ZvZKv3+azewLqVzMzH5kZvVmtjnBvqVmtt3MdpjZbfHNHwfWuPtK4KpUrimZqzUSTeszEH1KC/O4/oIzycnRE9QiQzlpQLj7dndf4O4LgPcC7cAj/Y8xsyozKxuw7ewEp1sNLB240cxyge8CVwBzgWvNbC5QA+yNH9Zzslolu7Scoi4mERmeZMcgLgV2uvvbA7ZfDDxqZoUAZrYS+M7AN7v7c0BDgvMuBna4+y537wIeAJYBtcRCYtBazexjZnZXU1NTkj+KhFnfcqPpWixIRJKXbEAsB+4fuNHdHwKeBB40s+uAG4GrkzhvNcdbChALhmrgYeATZrYKWJvoje6+1t1vKS/XnDrZpKO7h15P7zQbIpKcYf/2mVkBsXGALyfa7+53mNkDwCpgpru3JjouGe7eBqwY6Xkk8/TNw6QuJpHgJNOCuALY6O51iXaa2RJgHrHxiduTrGMfMLXf65r4Nhml0r0etYgkL5mAuJYE3UsAZrYQuIvYuMEKYLyZfT2Jc78MzDKzGfGWynLgsSTeL1nm2GJBBQoIkaAMKyDizyN8mNiYQCIlwDXuvtPde4EbgIED2ZjZ/cALwBwzqzWzmwDcPQrcSmwc43XgZ+6+JdkfRrLHsS4mtSBEAjOs3774WMD4IfY/P+B1N3B3guOuHeIcjwOPD6ceyX7HVpPTGIRIYDTVhoSSAkIkeAoICaW2NK9HLSLJU0BIKLWoBSESOAWEhFJrJEpejlGYp6+oSFD02yeh1BZfCyKdy42KSHIUEBJKLZ1RPQMhEjAFhIRSa0QT9YkETQEhodTWpam+RYKmgJBQOlWLBYnI8CkgJJRa44PUIhIcBYSEUmtnlDEapBYJlAJCQqk1ohaESNAUEBI6ke4e2rp6NEgtEjAFhITOCzuPALBgWkWwhYiMcgoICZ2nth6ktCCX980cdIZ5ETkNFBASKr29ztNb67lkThWFeblBlyMyqikgJFQ27W3kcGsnHzlnYtCliIx6CggJlae2HiQvx7hkTlXQpYiMegoICZWnt9ZxwVnjKS/OD7oUkVFPASGhsaO+lV2H2tS9JBISCggJjae2HgTgsncrIETCQAEhofH01jrOrS5nSkVx0KWICAoICYn65gib9jTykblqPYiEhQJCQuF3bx4G4FJ1L4mEhgJCQmFfYwcAM6tKA65ERPooICQU6pojjCst0NPTIiGigJBQqGvupKqsMOgyRKQfBYSEQn1LhIlji4IuQ0T6UUBIKNQ1R9SCEAkZBYQErqfXOdzapRaESMgoICRwR9o66el1Jo5VC0IkTBQQErj65k4AqtSCEAkVBYQErq45AqAuJpGQUUBI4OriLQh1MYmEiwJCAlfXHMEMJoxRQIiEiQJCAlffEmF8aSH5ufo6ioSJfiMlcHXNnepeEgkhBYQErq5ZT1GLhFGoA8LMSs3sXjO728yuC7oeOTU0D5NIOA0rIMyswszWmNk2M3vdzC5M5WJm9iMzqzezzQn2LTWz7Wa2w8xui2/+OLDG3VcCV6VyTQm37p5ejrR16hkIkRAabgviP4En3P1dwHzg9f47zazKzMoGbDs7wXlWA0sHbjSzXOC7wBXAXOBaM5sL1AB744f1DLNWySCHWztx1y2uImF00oAws3LgA8APAdy9y90bBxx2MfComRXG37MS+M7Ac7n7c0BDgsssBna4+y537wIeAJYBtcRCYtBazexjZnZXU1PTyX4UCaFjz0CUqQUhEjbDaUHMAA4B95jZJjP7gZmdsOyXuz8EPAk8GB8ruBG4Ook6qjneUoBYMFQDDwOfMLNVwNpEb3T3te5+S3l5eRKXk7DQU9Qi4TWcgMgDzgNWuftCoA24beBB7n4HEAFWAVe5e+tIi3P3Nndf4e6fdff7Rno+CZ/6Fj1FLRJWwwmIWqDW3dfFX68hFhgnMLMlwDzgEeD2JOvYB0zt97omvk2yXH1zhByD8XqKWiR0ThoQ7n4Q2Gtmc+KbLgW29j/GzBYCdxEbN1gBjDezrydRx8vALDObYWYFwHLgsSTeLxmqrjnCGWWF5OZY0KWIyADDvYvpb4D7zOxVYAHwjQH7S4Br3H2nu/cCNwBvDzyJmd0PvADMMbNaM7sJwN2jwK3ExjFeB37m7ltS+Hkkw8Seotb4g0gY5Q3nIHd/BVg0xP7nB7zuBu5OcNy1Q5zjceDx4dQj2aOuOUJNZUnQZYhIAqF+klqyX31LJ1UaoBYJJQWEBKYz2kNDW5eegRAJKQWEBOaQbnEVCTUFhATm+EpyakGIhJECQgJTH3+KWmMQIuGkgJDAaJoNkXBTQEhg6lo6ycsxxpUUBF2KiCSggJDA1DVHqCorJEdPUYuEkgJCAlPfrIWCRMJMASGBqT3aTnVFcdBliMggFBASiK5oL3uPdjBjQunJDxaRQCggJBB7j7bT0+sKCJEQU0BIIHYfagNgxhkKCJGwUkBIIHYfjgXEWWpBiISWAkICsetwG5Ul+VToGQiR0FJASCB2H27V+INIyCkghrC/sYNXaxtx96BLyTpvHW5nugJCJNSGtaLcaNTeFWX5XS+yp6GdaeNKWLZgClfMm8yUiiLGFOaRl5tDU0c3b9a1sL2uhV2H2th3tIN9jR20dkb5xz99N5fNnRj0jxFKbZ1RDjZHNP4gEnIKiEF888k32NPQzt9dNpuX32rgu7/dwXd+s+PY/uL8XDq6e469LsrPobqimJrKEiLdPXzupxv58Y2LOf+s8UGUH2pvHYnfwTRhTMCViMhQFBAJbHi7gXv+sJsbLjyTz182C4hNTf38zsM0tnfT3BGlOdLN+DEFzJlYxpxJZVRXFGMWm1Oooa2Lq7//B27+8XoevOVC5k4ZG+SPEzp9dzBpDEIk3EZ9QBxt66K7p/fYnECR7h6+uOZVppQX8w9L33XsuKqxRfz5wpphnXNcaQE/vul8PrnqD3z6npf4+Wfex7TxJaek/kzU9wzE9An6TETCbNQHxM/W7+VffrWNd08eywdmT6ChtYtdh9r475sWM6Yw9Y+nuqKY/75pMZ/8/gt86eevcv8tF6Sx6sy2+3Abk8uLKCkY9V8/kVAb9b+hH547kV6H5944xI9+v5vuHueaRTUsmXXGiM99dlUZK5ecxb89uZ29De1MHae/mCH2DIS6l0TCb9QHxFlnjOGzl4zhs5fMpLUzyqt7GznvzMq0nf/PF1bzzae28/ONtXzhstlpO28me+tIGx89d3LQZYjISeg5iH7GFObxvrMnUJSfm7ZzTqko5v1nT2DNhlp6e/U8xdG2Lhrbu9WCEMkACojT4JPvraH2aAfrdjcEXUrgdukOJpGMoYA4DT4ydxJlhXms2VCb9nPvqG+lJ4NaJrrFVSRzKCBOg+KCXK6cP5lfbT5AW2c0qff29jpffWwL//Pq/nfse622iQ9/+1ke3pj+4DlVdh9uJTfHNGAvkgEUEKfJJ99bQ3tXD4+/diCp9616dier//AWX31sC5F+T24D3PP8btxh/VtH01nqKbX7cBvTxpWQn6uvnkjY6bf0NDlvWiUzJpTy05f2sG7XEf6w4zC/f/Mw7V2Dtyhe2t3Avz+1nffUlHO4tYuf92sp1LdEWBtvVfyxtvFUl582uw7pFleRTKGAOE3MjKsX1bBpTyOfuutF/uIH67j+h+u47N+f5YnNB98xY2xDWxd/e/8mpo0r4b6bz2d+TTl3P7fr2HjDT9ftobvHuWr+FN6oaxkyaMKit9d5+0i7AkIkQ4z65yBOp5vffxbzayoAyMsxWiJRvvnUdj7zkw186F1VfObimeTmQLTH+d4zO2lo6+Lhv34fZUX5/NXFM/nr+zby1JaDfOjdVfzkxT18cM4Z/NnCKTz2x/1s3tfM4hnjgv0BE/j+szv596e2U5CbQ1F8gkNN8y2SGRQQp1FBXg4XnT3hhG2XzDmD1X94i289/Qa/2VZ/wr5/XnYO86rLAbj8nEmcOb6E7z+7k47uHg63dvKXF83gnPhEgK/WNoYuIOqbI/zXr9/k3OpyFk6rJNLdQ6/D5ZoGXSQjKCAClpebw81LzuKq+VPYvL+J3Jwc8nKMypKCE2aBzc0xVi45i688uplvPP46M88o5QOzJmBmVFcU88rexsB+hq5oLxvePsr5M8aRk2PHtv/Hr9+kK9rLt65ZoFaDSAZSQIRE1dgiPhSfUXYwn3xvDd9++g0Ot3bx+ctmH5tefP7U8sAGqiPdPfz1fRv5zbZ6brjwTL521TmYGTvqW3nw5b1cf/40hYNIhtIgdQYpys/ls5fMZEp5EZ84r/rY9vk1Fext6KChreu01tPR1cPKH6/nN9vqWTJrAj9+4W3+7cntANzxxDaK83P5m0tnndaaRCR91ILIMDcvOYubl5x1wrb5UyuA2O2uH5xTdVrqaO+KctPq9by4+wh3fPI9XP3eGr7y6Ga+98xODjZFeGprHX//4dlMGFN4WuoRkfRTQGSBedXlmMEf956egDjU0skt/72eP+5t5FvXzD+2kNI/L5tHa2eUhzfto6qskJuWzDjltYjIqaOAyAJjCvOYVTWGV2ubTvm1th9s4cbVL3OkrZPvXXceS+cdn7Y7J8f45tXzmVRexEUzJ2hBIJEMp9/gLPGemgp+u60edz82eJ0OrZ1RWiNR2ruibD3QzG0/f43Swlwe+qv3cW5N+TuOz8/N4ctXvDtt1xeR4CggssT8qRWs2VDLvsYOaipHPhHevsYOvvbYFp7aWnfC9nnVY/nBDX/CpPKh77gSkcyngMgS8+N/zf9xb9NJA6K+OcLvdxwmPzeHwrzYE87jSguYMKaQ8uJ87n3hLf7zf98E4LOXzGRqZQklBbmUFeXxvpkTKC5I34JKIhJeCogs8a5JYynIzeGPtY189D2DL+e5t6Gda/7fCxxoigx5vg/PncjtH5ubltaIiGQmBUSWKMjLYV71WO55fjeb9zVxyZwz+NC7qji7quzYMXXNEa7/4TraOqPcv/ICzigrINLdS0d3D0dauzjc2smR1i7eM7X8tN0uKyLhZQNnEc1UixYt8vXr1wddRqB21Lfw0Ppantl+iO11LUCs62n54mm8/+wJ3Lj6ZfY3dvCTm89n4bTKgKsVkTAwsw3uvijhPgVEdtrf2MGvNh/kwZf38EZdKwCFeTnce+NiLjhrfMDViUhYDBUQ6mLKUlMqirnp/TO48aLpbNrbyC827ePycyYpHERk2BQQWc7MOG9aJeepS0lEkhTqgDCzUuB7QBfwjLvfF3BJIiKjxrBmczWzt8zsNTN7xcxS7ug3sx+ZWb2ZbU6wb6mZbTezHWZ2W3zzx4E17r4SuCrV64qISPKSme77g+6+INFghplVmVnZgG1nJzjHamBpgvfnAt8FrgDmAtea2VygBtgbP6wniVpFRGSE0rUexMXAo2ZWCGBmK4HvDDzI3Z8DGhK8fzGww913uXsX8ACwDKglFhKD1mpmHzOzu5qaTv1EdSIio8lwA8KBp8xsg5nd8o6d7g8BTwIPmtl1wI3A1UnUUc3xlgLEgqEaeBj4hJmtAtYmLMx9rbvfUl7+zonjREQkdcMdpH6/u+8zsyrgaTPbFm8NHOPud5jZA8AqYKa7t460OHdvA1aM9DwiIpK8YbUg3H1f/N/1wCPEuoROYGZLgHnx/bcnWcc+YGq/1zXxbSIiEpCTBoSZlfYNQMdvO/0IsHnAMQuBu4iNG6wAxpvZ15Oo42VglpnNMLMCYDnwWBLvFxGRNBtOF9NE4JH4IjR5wE/d/YkBx5QA17j7TgAzuwH4y4EnMrP7gUuACWZWC9zu7j9096iZ3UpsHCMX+JG7b0nmB9mwYcNhM3s7mfeE0ATgcNBFhIg+j+P0WZxIn8eJRvJ5nDnYjqyZiykbmNn6weZEGY30eRynz+JE+jxOdKo+j3Td5ioiIllGASEiIgkpIMLlrqALCBl9HsfpsziRPo8TnZLPQ2MQIiKSkFoQIiKSkAJCREQSUkAExMymmtlvzWyrmW0xs8/Ht48zs6fN7M34v0fNSj9mlmtmm8zsf+KvZ5jZuvgU8A/GH6IcFcyswszWmNk2M3vdzC4crd8NM/u7+O/IZjO738yKRtN3I9EyCYN9Fyzmv+Kfy6tmdt5Irq2ACE4U+Ht3nwtcAHwuPsX5bcCv3X0W8Ov469Hi88Dr/V7/K/Btdz8bOArcFEhVwfhP4Al3fxcwn9jnMuq+G2ZWDfwtsMjd5xF7kHY5o+u7sZp3LpMw2HfhCmBW/J9biM2NlzIFREDc/YC7b4z/dwux/wFUE5uu5N74YfcCfxZIgaeZmdUAHwV+EH9twIeANfFDRtNnUQ58APghgLt3uXsjo/S7QWwGh2IzyyM2a8MBRtF3Y5BlEgb7LiwDfuwxLwIVZjY51WsrIELAzKYDC4F1wER3PxDfdZDYVCejwX8A/wD0xl+PBxrdPRp/3TcF/GgwAzgE3BPvcvtBfB60UffdiE8U+k1gD7FgaAI2MHq/G30G+y4MtnRCShQQATOzMcDPgS+4e3P/fR67Bznr70M2syuBenffEHQtIZEHnAescveFQBsDupNG0XejkthfxTOAKUApCValHM1O5XdBAREgM8snFg73ufvD8c11fU3C+L/rg6rvNLoIuMrM3iK2muCHiPXBV8S7FWB0TQFfC9S6+7r46zXEAmM0fjcuA3a7+yF37ya2iNhFjN7vRp/BvgtpXTpBARGQeB/7D4HX3f1b/XY9Bnw6/t+fBn5xums73dz9y+5e4+7TiQ1A/sbdrwN+C3wyftio+CwA3P0gsNfM5sQ3XQpsZRR+N4h1LV1gZiXx35m+z2JUfjf6Gey78BhwQ/xupguApn5dUUnTk9QBMbP3A78DXuN4v/v/JTYO8TNgGvA2sWnUE63jnZXM7BLg/7j7lWZ2FrEWxThgE3C9u3cGWN5pY2YLiA3YFwC7iK2zksMo/G6Y2deATxG7828TcDOxfvVR8d3ov0wCUEdsQbZHSfBdiIfoncS64dqBFe6+PuVrKyBERCQRdTGJiEhCCggREUlIASEiIgkpIEREJCEFhIiIJKSAEBGRhBQQIiKS0P8HeA9W4IRiiU4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.semilogy(r_list, BE_eig_r, label=\"eig-r\")\n",
    "plt.semilogy(r_list, BE_svd, label=\"svd\")\n",
    "#plt.semilogy(r_list, BE_kry, label=\"kry\")\n",
    "#plt.semilogy(r_list, BE_lfd, label=\"lfd\")\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}