RL1 - MDPs and model-based policy search.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Class 1: MDPs and model-based policy search**\n",
    "\n",
    "Notation: MDP = Markov Decision Process\n",
    "1. [Everything you need to know](#everything)\n",
    "2. [Modeling sequential decision problems with MDPs](#modeling)\n",
    "    1. [My first Markov Decision Processes](#mdp)\n",
    "    2. [Some Gym notions](#gym)\n",
    "    3. [Policies](#policies)\n",
    "    4. [Value Functions](#values)\n",
    "    5. [Optimal policies](#optimal)\n",
    "    6. [Evaluation equation](#eval)\n",
    "    7. [Optimality equation](#optimality)\n",
    "    8. [Summary](#summary)\n",
    "3. [Optimizing policies - algorithms](#algos)\n",
    "  1. [Value Iteration](#vi)\n",
    "  2. [Policy Iteration](#pi)\n",
    "  3. [Policy optimization with linear programming](#lp)\n",
    "  4. [Asynchronous Dynamic Programming](#adp)\n",
    "4. [Model-free Reinforcement Learning](#rl)\n",
    "  1. [Taxonomy of Machine Learning](#taxonomy)\n",
    "  2. [Introduction to model-free Reinforcement Learning](#introrl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"everything\"></a>Everything you need to know\n",
    "\n",
    "\n",
    "Everything you should remember after this session.\n",
    "<div class=\"alert alert-success\">\n",
    "<ul>\n",
    "<li> A Markov Decision Process (MDP) is a 4-tuple $\\langle S,A,p,r \\rangle$. $S$ is the state space, $A$ is the action space, $p(s'|s,a)$ the transition model and $r(s,a)$ the reward model.\n",
    "<li> A policy is a mapping $\\pi:S\\rightarrow A$ that specifies what to do in a given state.\n",
    "<li> Value function of a policy: $V^\\pi(s)=\\mathbb{E}\\left( \\lim\\limits_{H\\rightarrow\\infty} \\sum\\limits_{t = 0}^H \\gamma^t r_t \\bigg| s_0 = s, \\pi \\right)$\n",
    "<li> State-action value function of a policy $Q^\\pi(s,a) = \\mathbb{E}\\left( \\lim\\limits_{H\\rightarrow\\infty} \\sum\\limits_{t = 0}^H \\gamma^t r_t \\bigg| s_0 = s, a_0=a, \\pi \\right)$\n",
    "<li> Optimal policy $\\pi^*$ dominates every other possible policy in every state: $\\forall s, \\ V^*(s) = V^{\\pi^*}(s) \\geq V^\\pi(s)$\n",
    "<li> Evaluation equation: $V^\\pi \\in \\mathcal{F}\\left(S,\\mathbb{R}\\right)$ is the only solution to $V\\left(s\\right) = (T^\\pi V)(s) = r\\left(s,\\pi\\left(s\\right)\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,\\pi\\left(s\\right)\\right) V\\left(s'\\right)$<br>\n",
    "Similarly, $Q^\\pi \\in \\mathcal{F}\\left(S\\times A,\\mathbb{R}\\right)$ is the only solution to $Q\\left(s,a\\right) = (T^\\pi Q)(s,a) = r\\left(s,a\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,a\\right) Q\\left(s', \\pi\\left(s'\\right)\\right)$\n",
    "<li> Properties of $T^\\pi$ 1) $V=T^\\pi V$ is a linear system, 2) $T^\\pi$ is a contraction mapping over $\\mathcal{F}(S,\\mathbb{R})$ or $\\mathcal{F}(S\\times A,\\mathbb{R})$.\n",
    "<li> Finding $V^\\pi$ or $Q^\\pi$: by matrix inversion or by repeatedly applying $T^\\pi$ to any initial function.\n",
    "<li> Optimality equation: $V^* \\in \\mathcal{F}\\left(S,\\mathbb{R}\\right)$ is the only solution to $V(s) = (T^* V) (s) = \\max\\limits_{a\\in A} \\left\\{ r(s,a) + \\gamma \\sum\\limits_{s'\\in S} p(s'|s,a) V(s') \\right\\}$<br>\n",
    "Similarly, $Q^* \\in \\mathcal{F}\\left(S\\times A,\\mathbb{R}\\right)$ is the only solution to $Q(s,a) = (T^* Q) (s,a) = r(s,a) + \\gamma \\sum\\limits_{s'\\in S}p(s'|s,a) \\max\\limits_{a'\\in A} Q(s',a')$\n",
    "<li> $T^*$ is a contraction mapping.\n",
    "<li> Value Iteration: finding $V^*$ or $Q^*$ by repeatedly applying $T^*$ to any initial function.\n",
    "<li> Policy Iteration: finding $V^*$ or $Q^*$ and $\\pi^*$ by building the sequence $\\pi_{n+1}(s) = \\arg\\max_{a\\in A} Q^{\\pi_n}(s,a)$ that converges to $\\pi^*$. Repeatedly alternates an evaluation and an improvement phase.\n",
    "<li> An OpenAI Gym environment has a standard API with 5 methods (step, reset, render, close, seed) and 3 attributes (action_space, observation_space, reward_range). Environment-specific features can be accessed via the \"unwrapped\" attribute.\n",
    "</ul>\n",
    "</div>\n",
    "\n",
    "Of course, all this seems very obscure right now and the block above will only serve as a reminder when you re-open the notebook later. We will introduce every concept intuitively and progessively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"modeling\"></a>Modeling sequential decision problems with MDPs\n",
    "\n",
    "## <a id=\"mdp\"></a>My first Markov Decision Process\n",
    "\n",
    "A few `import` first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gym\n",
    "import gym.envs.toy_text.frozen_lake as fl\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "np.set_printoptions(precision=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shall start with the [FrozenLake](https://gym.openai.com/envs/FrozenLake-v0/) game and work our way to more general concepts. It's also the occasion to familiarize with OpenAI Gym."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "env = gym.make('FrozenLake-v0')\n",
    "_=env.render()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at this problem's description (using for example `help(fl.FrozenLakeEnv)`). We read:\n",
    "\n",
    "`|  Winter is here. You and your friends were tossing around a frisbee at the park\n",
    "|  when you made a wild throw that left the frisbee out in the middle of the lake.\n",
    "|  The water is mostly frozen, but there are a few holes where the ice has melted.\n",
    "|  If you step into one of those holes, you'll fall into the freezing water.\n",
    "|  At this time, there's an international frisbee shortage, so it's absolutely imperative that\n",
    "|  you navigate across the lake and retrieve the disc.\n",
    "|  However, the ice is slippery, so you won't always move in the direction you intend.\n",
    "|  The surface is described using a grid like the following\n",
    "|  \n",
    "|      SFFF\n",
    "|      FHFH\n",
    "|      FFFH\n",
    "|      HFFG\n",
    "|  \n",
    "|  S : starting point, safe\n",
    "|  F : frozen surface, safe\n",
    "|  H : hole, fall to your doom\n",
    "|  G : goal, where the frisbee is located\n",
    "|  \n",
    "|  The episode ends when you reach the goal or fall in a hole.\n",
    "|  You receive a reward of 1 if you reach the goal, and zero otherwise.`\n",
    "\n",
    "So it's a game of navigation.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Questions:</b><br>What are the possible states of an agent in this game?<br> What are its possible actions?<br>How would you describe the result of action $a$ in state $s$?<br> How would you formalize mathematically the goal of an agent?</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\"><a href=\"#answers1\" data-toggle=\"collapse\"><b>Answers:</b></a><br>\n",
    "<div id=\"answers1\" class=\"collapse\">\n",
    "States set: the 16 positions on the map.<br>\n",
    "Actions set: the 4 actions $\\{$N,S,E,W$\\}$<br>\n",
    "$s'$ resulting from $(s,a)$ follows a distribution $P(s'|s,a)$<br>\n",
    "Since each time step yields a reward, reaching the goal $\\Leftrightarrow$ $\\max \\sum r_t$\n",
    "</div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's develop a general theory for describing environments such as the FrozenLake game. Here is the main result we will develop on the white board or the presentation:<br>\n",
    "<br>\n",
    "<div class=\"alert alert-success\"><b>Markov Decision Process (MDP)</b><br>\n",
    "A Markov Decision Process is given by:\n",
    "<ul>\n",
    "<li> A set of states $S$\n",
    "<li> A set of actions $A$\n",
    "<li> A (Markovian) transition model $p(s'|s,a)$\n",
    "<li> A reward model $r(s,a)$ or $r(s,a,s')$\n",
    "<li> A set of discrete decision epochs $T=\\{0,1,\\ldots,H\\}$\n",
    "</ul>\n",
    "</div>\n",
    "\n",
    "<a href=\"#moreMDP\" data-toggle=\"collapse\"> More?</a><br>\n",
    "<div id=\"moreMDP\" class=\"collapse\">\n",
    "<ol>\n",
    "<li> 1. Illustration of an MDP's components:\n",
    "<img src=\"images/transitions.png\">\n",
    "<li> 2. And the resulting dynamics:\n",
    "<img src=\"images/mdp2.png\">\n",
    "<li> 3. If $H\\rightarrow\\infty$ we have an infinite horizon planning problem.<br>\n",
    "<li> 4. We shall see that the two formulations $r(s,a)$ and $r(s,a,s')$ are actually equivalent. We will admit that for now and will only use $r(s,a)$.\n",
    "<li> 5. $S$ and $A$ may each be either:\n",
    "    <ul>\n",
    "    <li> arbitrary finite sets,\n",
    "    <li> arbitrary countable infinite sets,\n",
    "    <li> compact subsets of finite dimensional Euclidean space, or\n",
    "    <li> non-empty Borel subsets of complete, separable metric spaces.\n",
    "    </ul>\n",
    "<li> 6. Markov transition model: $p(s_{t+1}|s_t,a_t,s_{t-1},a_{t-1},\\ldots) = p(s_{t+1}|s_t,a_t)$\n",
    "</ol>\n",
    "</div>\n",
    "Since we will only work with infinite horizon problems, we shall identify the MDP with the 4-tuple $\\langle S,A,p,r\\rangle$.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Question</b><br>\n",
    "Now, can you identify the MDP underlying the FrozenLake game?\n",
    "</div>\n",
    "<br>\n",
    "Let's confirm that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(env.action_space)\n",
    "print(env.observation_space)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `print(help(gym.envs.toy_text.discrete))` one sees that all discrete state environments (such as FrozenLake) have the same inner attributes:\n",
    "\n",
    "`|  - nS: number of states\n",
    "|  - nA: number of actions\n",
    "|  - P: transitions (*)\n",
    "|  - isd: initial state distribution (**)\n",
    "|  \n",
    "|  (*) dictionary dict of dicts of lists, where\n",
    "|    P[s][a] == [(probability, nextstate, reward, done), ...]\n",
    "|  (**) list or array of length nS`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(env.unwrapped.P[1][0])\n",
    "print(env.unwrapped.P[1][0][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just a few utility stuff from the FrozenLake environment (so that you don't spend time on useless things):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "actions = {fl.LEFT: '\\u2190', fl.DOWN: '\\u2193', fl.RIGHT: '\\u2192', fl.UP: '\\u2191'}\n",
    "print(actions)\n",
    "\n",
    "def to_s(row,col):\n",
    "    return row*env.unwrapped.ncol+col\n",
    "\n",
    "def to_row_col(s):\n",
    "    col = s%env.unwrapped.ncol\n",
    "    row = int((s-col)/env.unwrapped.ncol)\n",
    "    return row,col"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"gym\"></a>Some Gym notions\n",
    "\n",
    "Gym environments all are encapsulated so that they have the same API. This can be described by calling `help(gym.Env)` Namely, they all have 5 methods and 3 attributes.\n",
    "\n",
    "Methods:\n",
    "- `step` --- Runs one timestep of the environment's dynamics. When end of episode is reached, you are responsible for calling `reset()`  to reset this environment's state. Accepts an action and returns a tuple (observation, reward, done, info).\n",
    "- `reset` --- Resets the state of the environment and returns an initial observation.\n",
    "- `render` --- Renders the environment. The set of supported modes varies per environment. See `help(gym.Env.render)` for details and arguments.\n",
    "- `close` --- Cleans everything up, shuts the light and closes the door.\n",
    "- `seed` --- Sets the seed for this env's random number generator(s).\n",
    "\n",
    "Attributes:\n",
    "- `action_space` --- The Space(*) object corresponding to valid actions.\n",
    "- `observation_space` --- The Space(*) object corresponding to valid observations.\n",
    "- `reward_range` --- A tuple corresponding to the min and max possible rewards.\n",
    "\n",
    "(*) `Space` is a class (`help(gym.Space)`) that defines a set, be it continuous or discrete, with some useful methods. It is used specifically as the type of the action and observation spaces.\n",
    "\n",
    "Beyond those common API methods, some class-specific methods or attributes can be accessed by using the `unwrapped` attribute of each environment (if present) as in the example above.\n",
    "\n",
    "General references on OpenAI Gym:\n",
    "- [The short description paper](https://arxiv.org/pdf/1606.01540.pdf)\n",
    "- [The github page](https://github.com/openai/gym/blob/master/gym/envs/toy_text/discrete.py) (the most useful resource!)\n",
    "- [The general introduction to Gym](https://gym.openai.com/docs/)\n",
    "- [The deprecated list of environments](https://gym.openai.com/envs/) (refer to the github page for the latest version)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"policies\"></a>Policies\n",
    "\n",
    "Now the question is: \"formally, how does one write the behaviour of an agent in the game?\". This behaviour provides the agent with an action given everything the agent has seen since the game started. Let's call such a behaviour a *policy* and define the broadest possible family of policies.\n",
    "\n",
    "At time step $t$, action $a$ is picked with probability $\\delta_t(a|h)$, where $h$ is the history of state and actions seen between the initial time step and $t$.\n",
    "\n",
    "<div class=\"alert alert-success\"><b>Policy $\\pi$</b><br>\n",
    "A policy $\\pi$ is a sequence of decision rules $\\delta_t$: $\\pi = \\{\\delta_t\\}_{t\\in T}$,<br>\n",
    "with $\\delta_t : \\left\\{\\begin{array}{ccc}\n",
    "S^{t+1}\\times A^{t} & \\rightarrow & \\mathcal{P}'(A)\\\\h & \\mapsto & \\delta_t(a|h)\n",
    "\\end{array}\\right.$\n",
    "</div>\n",
    "In other words:\n",
    "$\\delta_t(a|h)$ indicates<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; the distribution over action $a$<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; to undertake at time $t$, given<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; the history of states/actions $h$.\n",
    "\n",
    "This is a *history-dependent, non-stationary, stochastic* policy.\n",
    "\n",
    "Ok, that's quite a complicated definition but we said we wanted to be generic so we will keep that for now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"values\"></a>Value functions\n",
    "\n",
    "Now we want to find the best strategy for our game, ie. the best policy, so the question is \"how can one say that one policy is better than another?\". Consequently, we need to define criteria over policies.\n",
    "\n",
    "The baseline idea is to say that the criterion should reflect what the agent gains by applying the policy. Two problems arise:\n",
    "- this gain depends on the initial state\n",
    "- two applications of the policy can result in two different trajectories since the agent's environment is stochastic (it's an MDP, remember)\n",
    "\n",
    "Let us introduce three criteria as examples to give a general intuition of various possibilities. Bear in mind that many more (exotic) criteria are possible, they simply express what we consider an optimal behaviour. The three criteria below are actually functions, that map an **initial state** $s$ to the **expected** return the agent gets by applying the policy from $s$.\n",
    "\n",
    "<table>\n",
    "<tr>\n",
    "    <td> Average reward </td>\n",
    "    <td width=\"300px\">$V(s) = \\mathbb{E}\\left( \\lim\\limits_{H\\rightarrow\\infty}  \\frac{1}{H} \\sum\\limits_{t = 0}^H r_t \\bigg| s_0 = s \\right)$ </td>\n",
    "    <td></td>\n",
    "</tr>\n",
    "<tr>\n",
    "    <td> Total reward </td>\n",
    "    <td>$V(s) = \\mathbb{E}\\left( \\lim\\limits_{H\\rightarrow\\infty} \\sum\\limits_{t = 0}^H r_t \\bigg| s_0 = s \\right)$ </td>\n",
    "    <td></td>\n",
    "</tr>\n",
    "<tr>\n",
    "    <td>Discounted reward </td>\n",
    "    <td>$V(s) = \\mathbb{E}\\left( \\lim\\limits_{H\\rightarrow\\infty} \\sum\\limits_{t = 0}^H \\gamma^t r_t \\bigg| s_0 = s \\right)$ </td>\n",
    "    <td width=\"150px\">with $0\\leq \\gamma<1$</td>\n",
    "</tr>\n",
    "</table>\n",
    "\n",
    "- The average reward criterion characterizes the average reward per time step the agent gets. This can be useful in some control applications. However, in the case of FrozenLake, we don't want to average our rewards, we want to get our frisbee as soon as possible.\n",
    "- The total reward criterion seems more adapted: it maximizes the cumulated rewards obtained during an episode. But it does not disciminate whether they were obtained at the beginning or late in the episode. Additionnaly, it suffers from a major flaw: for infinite horizon problems, even if the reward model is bounded, this sum might diverge. So we need a better formulation for the general case of infinite horizon problems.\n",
    "- The discounted reward criterion suits our needs. The gamma factor ($0\\leq \\gamma<1$) guarantees that with bounded reward models $r$, the sum always converges. Also it has the properties we desire: a reward of 1 obtained at the first time step weights 1 in the final criterion, while a reward of 1 obtained after $t$ time steps only weights $\\gamma^t$; it is *discounted* by $\\gamma^t$ (hence the criterion's name).\n",
    "\n",
    "From now on we shall concentrate on the discounted reward criterion. Given this criterion, we can introduce the notion of a value function.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-success\"><b>Value function $V^\\pi$ of a policy $\\pi$ under a $\\gamma$-discounted criterion</b><br>\n",
    "$$V^\\pi : \\left\\{\\begin{array}{ccl}\n",
    "S & \\rightarrow & \\mathbb{R}\\\\\n",
    "s & \\mapsto & V^\\pi(s)=\\mathbb{E}\\left( \\lim\\limits_{H\\rightarrow\\infty} \\sum\\limits_{t = 0}^H \\gamma^t r_t \\bigg| s_0 = s, \\pi \\right)\\end{array}\\right. $$\n",
    "</div>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Exercice</b><br>\n",
    "Use the FrozenLake environment we've introduced earlier to obtain a Monte-Carlo estimate of $V^\\pi(s_0)$ over 100000 trials, with $s_0$ being the initial state and $\\pi$ being a simple policy that always goes right. Take $\\gamma = 0.9$.\n",
    "</div>\n",
    "Note that $\\gamma^{200} \\sim 10^{-9}$ so any reward obtained after 200 time steps will have a negligible contribution to $V^\\pi(s_0)$, thus rolling an episode out for 200 time steps should be sufficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_episodes = 100000\n",
    "horizon = 200\n",
    "gamma = 0.9\n",
    "Vepisode = np.zeros(nb_episodes)\n",
    "for i in range(nb_episodes):\n",
    "    env.reset()\n",
    "    for t in range(horizon):\n",
    "        next_state, r, done,_ = env.step(fl.RIGHT)\n",
    "        Vepisode[i] += gamma**t * r\n",
    "        if done:\n",
    "            break\n",
    "print(\"value estimate:\", np.mean(Vepisode))\n",
    "print(\"value variance:\", np.std(Vepisode))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"optimal\"></a>Optimal policies\n",
    "\n",
    "The fog clears up a bit: we can now compare policies given an initial state. We can now define what an optimal policy is.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-success\"><b>Optimal policy $\\pi^*$</b><br>\n",
    "$\\pi^*$ is said to be optimal iff $\\pi^* \\in \\arg\\max\\limits_{\\pi} V^\\pi$.<br>\n",
    "<br>\n",
    "    \n",
    "A policy is optimal if it **dominates** over any other policy in every state:\n",
    "$$\\pi^* \\textrm{ is optimal}\\Leftrightarrow \\forall s\\in S, \\ \\forall \\pi, \\ V^{\\pi^*}(s) \\geq V^\\pi(s)$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now get to our first fundamental result. Fortunately for us...<br>\n",
    "<br>\n",
    "<div class=\"alert alert-success\"><b>Optimal policy theorem</b><br>\n",
    "For $\\left\\{\\begin{array}{l}\n",
    "\\gamma\\textrm{-discounted criterion}\\\\\n",
    "\\textrm{infinite horizon}\n",
    "\\end{array}\\right.$, \n",
    "there always exists at least one optimal stationary, deterministic, Markovian policy.\n",
    "</div>\n",
    "\n",
    "Let's explain a little:\n",
    "- Markovian : $\\left\\{\\begin{array}{l}\n",
    "\\forall \\left(s_i,a_i\\right)\\in \\left(S\\times A\\right)^{t-1}\\\\\n",
    "\\forall \\left(s'_i,a'_i\\right)\\in \\left(S\\times A\\right)^{t-1}\n",
    "\\end{array}\\right., \\delta_t\\left(a|s_0, a_0, \\ldots, s_t\\right) = \\delta_t\\left(a|s'_0, a'_0, \\ldots, s_t\\right)$.\\\\\n",
    "One writes $\\delta_t(a|s)$.\n",
    "- Stationary : $\\forall (t,t')\\in \\mathbb{N}^2, \\delta_t = \\delta_t'$.\\\\\n",
    "One writes $\\pi = \\delta_0$.\n",
    "- Deterministic : $\\delta_t(a|h) = \\left\\{\\begin{array}{l}\n",
    "1\\textrm{ for a single }a\\\\\n",
    "0\\textrm{ otherwise}\n",
    "\\end{array}\\right.$.\n",
    "\n",
    "So in simpler words, we know that among all possible optimal policies, at least one is a function $\\pi:S\\rightarrow A$.\n",
    "\n",
    "That helps a lot: we don't have to search for optimal policies in that complex family of history-dependent, stochastic, non-stationary policies, can simply search for a function $\\pi(s)=a$ that maps states to actions.\n",
    "\n",
    "<a href=\"#optPol\" data-toggle=\"collapse\"> Curious to know why?</a><br>\n",
    "<div id=\"optPol\" class=\"collapse\">\n",
    "    The proof (very simple but a little long) is in chapter 6 of <b>Markov Decision Processes</b> (book by Martin L. Puterman).<br>\n",
    "To give you the general flavour:\n",
    "<ul>\n",
    "<li> The infinite horizon leads to the existence of an optimal <b>stationary</b> policy: if the horizon is infinitely far, the optimal decision rule $n$ steps before the end if the same as the one $n+1$ steps before the end (watch out, this intuition can be very false)\n",
    "<li> The <b>Markovian</b> property of $p(s'|s,a)$ allows to get optimal Markovian policies.\n",
    "<li> The <b>deterministic</b> part is somehow more tricky but just note that this result only holds for single-player MDPs. For a two-agents competitive game for example (like poker for instance), there is no deterministic optimal policy.\n",
    "</ul>\n",
    "</div>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"eval\"></a>Evaluation equation\n",
    "\n",
    "Let's play a little bit (mathematically) with our new toys.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Question</b><br>\n",
    "What's the value of \"$a$ then $\\pi$ from $s$\"?\n",
    "</div>\n",
    "\n",
    "Note: the answer uses the discounted criterion's definition and the value function of $\\pi$.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-danger\"><a href=\"#answers2\" data-toggle=\"collapse\"><b>Answer:</b></a><br>\n",
    "<div id=\"answers2\" class=\"collapse\">\n",
    "\\begin{align*}\n",
    "Q^\\pi(s,a) & = \\mathbb{E}\\left( \\sum\\limits_{t=0}^\\infty \\gamma^t r\\left(s_t, a_t\\right) \\bigg| s_0 = s, a_0=a, \\pi \\right)\\\\\n",
    " & = r\\left(s,a\\right) + \\mathbb{E}\\left( \\sum\\limits_{t=1}^\\infty \\gamma^t r\\left(s_t, a_t\\right) \\bigg| s_0 = s, a_0=a, \\pi \\right)\\\\\n",
    " & = r\\left(s,a\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,a\\right) \\mathbb{E}\\left( \\sum\\limits_{t=1}^\\infty \\gamma^{t-1} r\\left(s_t, a_t\\right) \\bigg| s_1 = s', \\pi \\right)\\\\\n",
    " & = r\\left(s,a\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,a\\right) V^\\pi\\left(s'\\right)\n",
    "\\end{align*}\n",
    "</div>\n",
    "</div>\n",
    "\n",
    "Note that we have introduced the very important state-action value function $Q^\\pi$.\n",
    "<div class=\"alert alert-success\"><b>State-action value function</b><br>\n",
    "$$Q^\\pi(s,a) = \\mathbb{E}\\left( \\sum\\limits_{t=0}^\\infty \\gamma^t r\\left(s_t, a_t\\right) \\bigg| s_0 = s, a_0=a, \\pi \\right)$$\n",
    "</div>\n",
    "\n",
    "That's interesting: given a policy $\\pi$, the best one-step lookahead action can be selected by maximizing $Q^\\pi$.\n",
    "- To improve on a policy $\\pi$, it is more useful to know $Q^\\pi$ than $V^\\pi$ and pick the *greedy* action.\n",
    "- Also $V^\\pi(s) = Q^\\pi(s,\\pi(s))$. Let's replace that above and we get an important equation to characterize $V^\\pi$.\n",
    "\n",
    "<div class=\"alert alert-success\"><b>Evaluation equation</b><br>\n",
    "$V^\\pi$ is a solution to the linear system:\n",
    "\\begin{gather*}\n",
    "V^\\pi\\left(s\\right) = r\\left(s,\\pi\\left(s\\right)\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,\\pi\\left(s\\right)\\right) V^\\pi\\left(s'\\right)\\\\\n",
    "V^\\pi = r^\\pi + \\gamma P^\\pi V^\\pi = T^\\pi V^\\pi\n",
    "\\end{gather*}\n",
    "Similarly:\n",
    "\\begin{gather*}\n",
    "Q^\\pi\\left(s,a\\right) = r\\left(s,a\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,a\\right) Q^\\pi\\left(s', \\pi\\left(s'\\right)\\right)\\\\\n",
    "Q^\\pi = r + \\gamma P Q^\\pi = T^\\pi Q^\\pi\n",
    "\\end{gather*}\n",
    "</div>\n",
    "\n",
    "<a href=\"#morePpi\" data-toggle=\"collapse\"> What are $P^\\pi$, $r^\\pi$ and $T^\\pi$ here?</a><br>\n",
    "<div id=\"morePpi\" class=\"collapse\">\n",
    "When the state space is discrete, $V^\\pi$ is vector of size $|S|$, $P^\\pi$ is a matrix containing the values $P^\\pi_{ij} = p\\left(s_j|s_i,\\pi(s_i)\\right)$ and, similarly, $r^\\pi$ is a vector containing the values $r^\\pi_i = r(s_i,\\pi(s_i))$. In better words, $P^\\pi$ is the <i>transition kernel</i> of the Markov chain describing the state dynamics under policy $\\pi$ and $r^\\pi$ is the associated reward model.<br>\n",
    "<br>\n",
    "This generalizes straightforwardly to the continuous action case: $V^\\pi$ is a function in the $\\mathcal{F}(S,\\mathbb{R})$ function space (the generalization of the vector in the previous sentence), $r^\\pi$ becomes the function $s\\mapsto r(s,\\pi(s))$ and  $P^\\pi$ becomes the operator over $\\mathcal{F}(S,\\mathbb{R})$ that maps function $V$ to $s\\mapsto \\int\\limits_{s'\\in S} p\\left(s'|s,\\pi\\left(s\\right)\\right) V\\left(s'\\right)ds'$.<br>\n",
    "<br>\n",
    "In the same fashion, one can define the $T^\\pi$ operator. In the discrete state space case, $T^\\pi$ is the linear operator in $\\mathbb{R}^{|S|}$ that maps $V^\\pi$ to $r^\\pi + \\gamma P^\\pi V^\\pi$.<br>\n",
    "<br>\n",
    "In the continuous state space case, $T^\\pi$ is the linear operator in $\\mathcal{F}(S,\\mathbb{R})$ that maps $V^\\pi$ to the function $s\\mapsto r\\left(s,\\pi\\left(s\\right)\\right) + \\gamma \\sum\\limits_{s'\\in S} p\\left(s'|s,\\pi\\left(s\\right)\\right) V^\\pi\\left(s'\\right)$.\n",
    "</div>\n",
    "\n",
    "<a href=\"#moreEval\" data-toggle=\"collapse\"> A remark</a><br>\n",
    "<div id=\"moreEval\" class=\"collapse\">\n",
    "For stochastic policies:\n",
    "$\\forall s\\in S, \\quad V^\\pi(s) = \\sum\\limits_{a\\in A} \\pi(a|s) \\left(r(s,a) + \\gamma \\sum\\limits_{s'\\in S} p(s'|s,a)V^\\pi(s') \\right)$\n",
    "</div>\n",
    "\n",
    "We've gone far from our original FrozenLake problem. Let's make all this very concrete:\n",
    "- A policy $\\pi$ is an agent's behaviour\n",
    "- In every state $s$, one can expect to gain $V^\\pi(s)$ in the long run by applying $\\pi$\n",
    "- The function $V^\\pi$ actually obeys the linear system of equations above that simply link the value of a state with the values of its successors in an episode.\n",
    "\n",
    "We can stop for a minute on the $T^\\pi$ evaluation operator (that maps a function $S\\rightarrow\\mathbb{R}$ to a function $S\\rightarrow\\mathbb{R}$) and the search for $V^\\pi$.<br>\n",
    "<br>\n",
    "\n",
    "<div class=\"alert alert-success\"><b>Properties of $T^\\pi$</b><br>\n",
    "<ol>\n",
    "<li> $T^\\pi$ is linear.<br>\n",
    "$\\Rightarrow$ Solving $V^\\pi = T^\\pi V^\\pi$ and $Q^\\pi = T^\\pi Q^\\pi$ by matrix inversion?<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;With $\\gamma<1$, $V^\\pi = \\left(I-\\gamma P^\\pi\\right)^{-1}r^\\pi$ and $Q^\\pi = \\left(I-\\gamma P\\right)^{-1}r^\\pi$\n",
    "<li> With $\\gamma<1$, $T^\\pi$ is a $\\| \\cdot \\|_\\infty$-contraction mapping over the $\\mathcal{F}(S,\\mathbb{R})$ (resp. $\\mathcal{F}(S\\times A,\\mathbb{R})$) Banach space.<br>\n",
    "$\\Rightarrow$ With $\\gamma<1$, $V^\\pi$ (resp. $Q^\\pi$) is the unique solution to the (linear) fixed point equation:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$V=T^\\pi V$ (resp. $Q=T^\\pi Q$).\n",
    "</ol>\n",
    "</div>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Exercice</b><br>\n",
    "Use the first property above to compute $V^\\pi$ for the policy that always moves right. To do this, you'll need to compute $r^\\pi$ and $P^\\pi$. Again, $\\gamma = 0.9$.\n",
    "Check if your result for $V^\\pi(s_0)$ is consistent with the Monte Carlo estimate.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 0.9\n",
    "r_pi0 = np.zeros((env.observation_space.n))\n",
    "P_pi0 = np.zeros((env.observation_space.n, env.observation_space.n))\n",
    "for s in range(env.observation_space.n):\n",
    "    outcomes = env.unwrapped.P[s][fl.RIGHT]\n",
    "    for o in outcomes:\n",
    "        p  = o[0]\n",
    "        s2 = o[1]\n",
    "        r  = o[2]\n",
    "        P_pi0[s][s2] += p\n",
    "        r_pi0[s] += p*r\n",
    "I = np.eye(env.observation_space.n)\n",
    "V_pi0 = np.dot(np.linalg.inv(I - gamma*P_pi0), r_pi0)\n",
    "print(V_pi0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could generalize that to a function that takes a policy as input. We'll suppose in this case that the policy is an array of actions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def policy_eval_lin(pi):\n",
    "    # build r and P\n",
    "    r_pi = np.zeros((env.observation_space.n))\n",
    "    P_pi = np.zeros((env.observation_space.n, env.observation_space.n))\n",
    "    for x in range(env.observation_space.n):\n",
    "        outcomes = env.unwrapped.P[x][pi[x]]\n",
    "        for o in outcomes:\n",
    "            p = o[0]\n",
    "            y = o[1]\n",
    "            r = o[2]\n",
    "            P_pi[x,y] += p\n",
    "            r_pi[x] += r*p\n",
    "    # Compute V\n",
    "    I = np.eye(env.observation_space.n)\n",
    "    return np.dot(np.linalg.inv(I - gamma*P_pi), r_pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi0 = fl.RIGHT*np.ones((env.observation_space.n))\n",
    "\n",
    "V_pi0 = policy_eval_lin(pi0)\n",
    "print(V_pi0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\"><b>Exercice</b><br>\n",
    "Use the second property above to compute $V^\\pi$ for the policy that always moves right. To do this, you'll need to remember that since $T^\\pi$ is a contraction mapping, the sequence $V_{n+1}=T^\\pi V_n$ converges to $T^\\pi$'s fixed point (which happens to be $V^\\pi$ according to the property). Again, $\\gamma = 0.9$. For a start, apply $T^\\pi$ for a fixed number of steps `max_iter`.\n",
    "Check if your result is consistent with the previous estimate.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def policy_eval_iter(pi, max_iter):\n",
    "    V = np.zeros((env.observation_space.n))\n",
    "    W = np.zeros((env.observation_space.n))\n",
    "    for i in range(max_iter):\n",
    "        for x in range(env.observation_space.n):\n",
    "            W[x]=0\n",
    "            outcomes = env.unwrapped.P[x][pi[x]]\n",
    "            # W[x] = sum_y P[y] (r_y + gamma V[y])\n",
    "            for o in outcomes:\n",
    "                p = o[0]\n",
    "                y = o[1]\n",
    "                r = o[2]\n",
    "                W[x] += p * (r+gamma*V[y])\n",
    "        V[:] = W\n",
    "    return V\n",
    "\n",
    "V_pi0 = policy_eval_iter(pi0,10000)\n",
    "print(V_pi0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could make a matrix-operations version of the code above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def policy_eval_iter_mat(pi, max_iter):\n",
    "    # build r and P\n",
    "    r_pi = np.zeros((env.observation_space.n))\n",
    "    P_pi = np.zeros((env.observation_space.n, env.observation_space.n))\n",
    "    for x in range(env.observation_space.n):\n",
    "        outcomes = env.unwrapped.P[x][pi[x]]\n",
    "        for o in outcomes:\n",
    "            p = o[0]\n",
    "            y = o[1]\n",
    "            r = o[2]\n",
    "            P_pi[x,y] += p\n",
    "            r_pi[x] += r*p\n",
    "    # Compute V\n",
    "    V = np.zeros((env.observation_space.n))\n",
    "    for i in range(max_iter):\n",
    "        V = r_pi + gamma * np.dot(P_pi, V)\n",
    "    return V\n",
    "\n",
    "V_pi0 = policy_eval_iter_mat(pi0,10000)\n",
    "print(V_pi0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Between two iterations in the algorithm above, the distance between $V_{n+1}$ and $V_n$ decreases as $\\|V_{n+1}-V_n\\| = \\|r^\\pi + \\gamma P^\\pi V_n - V_n\\|$. Since $T^\\pi$ is a contraction mapping, we have $\\|V_{n+1}-V_n\\| \\leq \\|V_{n}-V_{n-1}\\|$. Let's call this distance at time step $n$ the **residual**. Then the successive residuals monotonically tend to zero.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Exercice:</b><br>\n",
    "Now, can you use the property on the residuals to replace `max_iter` by a precision parameter `epsilon` that specifies the maximum error on $V^\\pi$? Advice: still keep `max_iter` to stop the computation in case you specify an `epsilon` that is too small. Return both $V^\\pi$ and the sequence of residuals. Plot the sequence of residuals and display the number of iterations necessary to reach the chose precision `epsilon`.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def policy_eval_iter_mat2(pi, epsilon, max_iter):\n",
    "    # build r and P\n",
    "    r_pi = np.zeros((env.observation_space.n))\n",
    "    P_pi = np.zeros((env.observation_space.n, env.observation_space.n))\n",
    "    for x in range(env.observation_space.n):\n",
    "        outcomes = env.unwrapped.P[x][pi[x]]\n",
    "        for o in outcomes:\n",
    "            p = o[0]\n",
    "            y = o[1]\n",
    "            r = o[2]\n",
    "            P_pi[x,y] += p\n",
    "            r_pi[x] += r*p\n",
    "    # Compute V\n",
    "    V = np.zeros((env.observation_space.n))\n",
    "    W = np.zeros((env.observation_space.n))\n",
    "    residuals = np.zeros((max_iter))\n",
    "    for i in range(max_iter):\n",
    "        W = r_pi + gamma * np.dot(P_pi, V)\n",
    "        residuals[i] = np.max(np.abs(W-V))\n",
    "        V[:] = W\n",
    "        if residuals[i]<epsilon:\n",
    "            residuals = residuals[:i+1]\n",
    "            break\n",
    "    return V, residuals\n",
    "\n",
    "V_pi0, residuals = policy_eval_iter_mat2(pi0,1e-4,10000)\n",
    "print(V_pi0)\n",
    "plt.plot(residuals)\n",
    "plt.figure()\n",
    "plt.semilogy(residuals)\n",
    "print(\"number of iterations:\", residuals.size)\n",
    "print(\"last residual\", residuals[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"optimality\"></a>Optimality equation\n",
    "\n",
    "Ok, so given a policy, we have a way (two actually) to compute its value function. Now we would like to find the optimal policy and we'd rather not enumerate all possible policies.\n",
    "\n",
    "Let's introduce some notations and start with two trivial remarks. We write:\n",
    "$$V^{\\pi^*} = V^*, \\quad Q^{\\pi^*} = Q^*$$\n",
    "\n",
    "Then we can note that:\n",
    "1. $Q^*\\left(s,a\\right) = r\\left(s,a\\right) + \\gamma \\sum\\limits_{s'\\in S}p\\left(s'|s,a\\right) V^*\\left(s'\\right)$\n",
    "2. If $\\pi^*$ is an optimal policy, then $V^*\\left(s\\right) = Q^*\\left(s,\\pi^*\\left(s\\right)\\right)$\n",
    "\n",
    "From the first remark, we note that the maximum expected return one can get from a trajectory initiated in $s$ starts with the action $a$ that maximizes $Q^*(s,a)$. So:\n",
    "<div class=\"alert alert-success\"><b>Optimal greedy policy</b><br>\n",
    "Any policy $\\pi$ defined by $\\pi(s) \\in \\arg\\max\\limits_{a\\in A} Q^*(s,a)$ is an optimal policy.\n",
    "</div>\n",
    "\n",
    "Consequently, finding $\\pi^*$ is equivalent to finding $Q^*$.\n",
    "\n",
    "Here comes the key theorem of this class:\n",
    "<div class=\"alert alert-success\"><b>Bellman optimality equation</b><br>\n",
    "The optimal value function obeys:\n",
    "$$V^*(s) = \\max\\limits_{a\\in A} \\left\\{ r(s,a) + \\gamma \\sum\\limits_{s'\\in S} p(s'|s,a) V^*(s') \\right\\}=(T^* V^*) (s)$$\n",
    "or in terms of $Q$-functions:\n",
    "$$Q^*(s,a) = r(s,a) + \\gamma \\sum\\limits_{s'\\in S}p(s'|s,a) \\max\\limits_{a'\\in A} Q^*(s',a')=(T^* Q^*) (s,a)$$\n",
    "</div>\n",
    "\n",
    "Why does that unfold intuitively from the previous remarks? Simply because of the Dynamic Programming principle: any subpath of an optimal path is itself optimal. That might still sound obscure. So let's say that if you start a trajectory in $s$ and if you know you will get the maximum value from any state $s'$ you reach, then the maximum value you can get from $s$ is precisely $\\max\\limits_{a\\in A} \\left\\{ r(s,a) + \\gamma \\sum\\limits_{s'\\in S} p(s'|s,a) V^*(s') \\right\\}$. The actual proof is quite trickier but the key ideas are here.<br>\n",
    "<br>\n",
    "As for the policy evaluation operator $T^\\pi$, we define the **Bellman optimality operator** $T^*$ as above.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-success\"><b>Properties of $T^*$</b><br>\n",
    "<ol>\n",
    "<li> $T^*$ is non-linear.<br>\n",
    "<li> With $\\gamma<1$, $T^*$ is a $\\| \\cdot \\|_\\infty$-contraction mapping over the $\\mathcal{F}(S,\\mathbb{R})$ (resp. $\\mathcal{F}(S\\times A,\\mathbb{R})$) Banach space.<br>\n",
    "$\\Rightarrow$ With $\\gamma<1$, $V^*$ (resp. $Q^*$) is the unique solution to the fixed point equation:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$V=T^* V$ (resp. $Q=T^* Q$).\n",
    "</ol>\n",
    "</div>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Exercice</b><br>\n",
    "Use the property above to compute $V^*$. To do this, you'll need to remember that since $T^*$ is a contraction mapping, the sequence $V_{n+1}=T^* V_n$ converges to $T^*$'s fixed point (which happens to be $V^*$ according to the property). Again, $\\gamma = 0.9$.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vf_optim(V, epsilon, max_iter):\n",
    "    W = np.copy(V)\n",
    "    residuals = np.zeros((max_iter))\n",
    "    for i in range(max_iter):\n",
    "        for s in range(env.observation_space.n):\n",
    "            Q = np.zeros((env.action_space.n))\n",
    "            for a in range(env.action_space.n):\n",
    "                outcomes = env.unwrapped.P[s][a]\n",
    "                for o in outcomes:\n",
    "                    p  = o[0]\n",
    "                    s2 = o[1]\n",
    "                    r  = o[2]\n",
    "                    Q[a] += p*(r+gamma*V[s2])\n",
    "            W[s] = np.max(Q)\n",
    "            #print(W[s])\n",
    "        residuals[i] = np.max(np.abs(W-V))\n",
    "        #print(\"abs\", np.abs(W-V))\n",
    "        np.copyto(V,W)\n",
    "        if residuals[i]<epsilon:\n",
    "            residuals = residuals[:i+1]\n",
    "            break\n",
    "    return V, residuals\n",
    "\n",
    "Vinit = np.zeros((env.observation_space.n))\n",
    "Vstar, residuals = vf_optim(Vinit, 1e-4, 10000)\n",
    "print(Vstar)\n",
    "plt.plot(residuals)\n",
    "plt.figure()\n",
    "plt.semilogy(residuals)\n",
    "print(\"number of iterations:\", residuals.size)\n",
    "print(\"last residual\", residuals[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"summary\"></a>Summary\n",
    "\n",
    "Let's wrap this whole section up. Our goal was to formally define the search for the best strategy for our game of FrozenLake. This has led us to formalizing the general **discrete-time stochastic optimal control problem**:\n",
    "- Environment (discrete time, non-deterministic, non-linear) $\\leftrightarrow$ MDP.\n",
    "- Behaviour $\\leftrightarrow$ control policy $\\pi : s\\mapsto a$.\n",
    "- Policy evaluation criterion $\\leftrightarrow$ $\\gamma$-discounted criterion.\n",
    "- Goal $\\leftrightarrow$ Maximize value function $V^\\pi(s)$, $Q^\\pi(s,a)$.\n",
    "- Evaluation eq. $\\leftrightarrow$ $V^\\pi = T^\\pi V^\\pi$, $Q^\\pi = T^\\pi Q^\\pi$.\n",
    "- Bellman optimality eq. $\\leftrightarrow$ $V^* = T^* V^*$, $Q^* = T^* Q^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"algos\"></a>Optimizing policies - algorithms\n",
    "\n",
    "From the analysis of the previous section, we introduce several algorithms to find $\\pi^*$.\n",
    "\n",
    "## <a id=\"vi\"></a>Value Iteration\n",
    "\n",
    "Value iteration is actually the algorithm we defined in function `vf_optim` above. It directly exploits the contraction mapping property of $T^*$ and iterates over value functions in order to converge to $V^*$. Once, $V^*$ is found, finding the optimal policy is a matter of writing $Q^*$ and defining a policy that is $Q^*$-greedy.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Exercice:</b><br>\n",
    "Compute and display an optimal policy for the FrozenLake game, using Value Iteration.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def value_iteration(V,epsilon,max_iter):\n",
    "    return vf_optim(V,epsilon,max_iter)\n",
    "\n",
    "def Q_from_V(V):\n",
    "    Q = np.zeros((env.observation_space.n, env.action_space.n))\n",
    "    for s in range(env.observation_space.n):\n",
    "        for a in range(env.action_space.n):\n",
    "            outcomes = env.unwrapped.P[s][a]\n",
    "            for o in outcomes:\n",
    "                p  = o[0]\n",
    "                s2 = o[1]\n",
    "                r  = o[2]\n",
    "                Q[s,a] += p*(r+gamma*V[s2])\n",
    "    return Q\n",
    "\n",
    "def greedyQpolicy(Q):\n",
    "    pi = np.zeros((env.observation_space.n),dtype=np.int)\n",
    "    for s in range(env.observation_space.n):\n",
    "        pi[s] = np.argmax(Q[s,:])\n",
    "    return pi\n",
    "\n",
    "def print_policy(pi):\n",
    "    for row in range(env.unwrapped.nrow):\n",
    "        for col in range(env.unwrapped.ncol):\n",
    "            print(actions[pi[to_s(row,col)]], end='')\n",
    "        print()\n",
    "    return\n",
    "\n",
    "Vinit = np.zeros((env.observation_space.n))\n",
    "Vstar,residuals = value_iteration(Vinit,1e-4,1000)\n",
    "Qstar = Q_from_V(Vstar)\n",
    "print(actions)\n",
    "print(Qstar)\n",
    "pi_star = greedyQpolicy(Qstar)\n",
    "print_policy(pi_star)\n",
    "_=env.render()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those familiar with the principles of Dynamic Programming will note that Value Iteration is a Dynamic Programming algorithm that operates in value function space, monotonically hopping from value function to value function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"pi\"></a>Policy Iteration\n",
    "\n",
    "The Policy Iteration algorithm stems from the following remark. Suppose we have a policy $\\pi$ and know its value function $V^\\pi$ and state-action value function $Q^\\pi$. Then, the non-stationary policy $\\pi'$ that acts greedily with respect to $Q^\\pi$ for the first time step and then follows $\\pi$ has a value function $V^{\\pi'}$ that is greater or equal to $V^\\pi$ (equal if $\\pi$ is optimal, strictly greater otherwise). Actually, the contraction property of $T^*$ insures that the stationary policy $\\pi'$ that is greedy with respect to $Q^\\pi$ is at least as good as $\\pi$, that is $V^{\\pi'}\\geq V^\\pi$. Consequently, the sequence of policies defined by $\\pi_{n+1}(s) = \\arg\\max_{a\\in A} Q^{\\pi_n}(s,a)$ has a monotonically improving corresponding sequence of value functions $V^{\\pi_n}$ and converges to $\\pi^*$.<br>\n",
    "<br>\n",
    "Let's make this more simple with a drawing. Policy iteration alternates two phases:\n",
    "1. Evaluate $\\pi_n$ $\\rightarrow Q^{\\pi_n}$\n",
    "2. Compute $\\pi_{n+1}$ as the $Q^{\\pi_n}$-greedy policy\n",
    "\n",
    "<img src=\"images/policyiteration.png\">\n",
    "\n",
    "<div class=\"alert alert-warning\"><b>Exercice:</b><br>\n",
    "Compute and display an optimal policy for the FrozenLake game, using Policy Iteration.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def policy_iteration(pi0,max_iter):\n",
    "    policies = np.zeros((max_iter, env.observation_space.n))\n",
    "    policies[0,:] = np.copy(pi0)\n",
    "    for i in range(max_iter-1):\n",
    "        Vpi = policy_eval_iter_mat(policies[i],500)\n",
    "        Qpi = Q_from_V(Vpi)\n",
    "        policies[i+1,:] = greedyQpolicy(Qpi)\n",
    "        if np.array_equal(policies[i,:],policies[i+1,:]):\n",
    "            policies = policies[:i,:]\n",
    "            break\n",
    "    return policies\n",
    "\n",
    "pi0 = fl.RIGHT*np.ones((env.observation_space.n))\n",
    "print_policy(pi0)\n",
    "policies = policy_iteration(pi0,10)\n",
    "print(\"number of iterations:\", policies.shape[0])\n",
    "print_policy(policies[-1,:])\n",
    "env.render()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warning:** the stopping criterion $\\pi_{n+1} = \\pi_n$ is good in most cases. However, it is not robust to multiple optimal policies! In this case, Policy Iteration might alternate between two (or more) equivalent optimal policies. This happens also when the computation of $Q^{\\pi_n}$ carries small errors. Even if its a bit more costly computationally, it might be safer to compare $V^{\\pi_{n+1}}$ and $V^{\\pi_n}$ and allow a precision error of $\\epsilon$. To observe this phenomenon, just replace the `policy_eval_iter_mat` call by a call to `policy_eval_lin` (which carries the errors of matrix inversion): the algorithm does not terminate anymore and alternates between equivalent optimal policies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def policy_iteration(pi0,max_iter):\n",
    "    policies = np.zeros((max_iter, env.observation_space.n))\n",
    "    policies[0,:] = np.copy(pi0)\n",
    "    for i in range(max_iter-1):\n",
    "        Vpi = policy_eval_lin(policies[i])\n",
    "        Qpi = Q_from_V(Vpi)\n",
    "        policies[i+1,:] = greedyQpolicy(Qpi)\n",
    "        if np.array_equal(policies[i,:],policies[i+1,:]):\n",
    "            policies = policies[:i,:]\n",
    "            break\n",
    "    return policies\n",
    "\n",
    "pi0 = fl.RIGHT*np.ones((env.observation_space.n))\n",
    "print_policy(pi0)\n",
    "policies = policy_iteration(pi0,10)\n",
    "print(\"number of iterations:\", policies.shape[0])\n",
    "print(policies[-2,:])\n",
    "print(policies[-1,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for Value Iteration, those familiar with Dynamic Programming will remark that Policy Iteration is a Dynamic Programming algorithm in policy space, monotonically hopping from policy to policy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"lp\"></a>Linear Programming\n",
    "\n",
    "If this is the first time you read this notebook, this part can be skipped. [Link to next part](#adp)\n",
    "\n",
    "An alternative way of finding $V^*$ is by casting the optimality equation as a linear optimization problem. This formulation is mainly given for your curiosity and we will not study it any further.<br>\n",
    "<br>\n",
    "Recall the optimality equation:\n",
    "$$\\forall s\\in S, V(s)=\\max\\limits_{a\\in A} \\left\\{r(s,a) + \\gamma \\sum\\limits_{s'\\in S} p(s'|s,a) V(s')\\right\\}$$\n",
    "The key remark to transform this into a linear program is to rephrase it as \"$V^*$ is the smallest value that dominates over all policy values\". This can be written as:\n",
    "$$\\left\\{ \\begin{array}{c}\n",
    "\\min \\sum\\limits_{s\\in S} V(s)\\\\\n",
    "s.t. \\ \\forall \\pi, \\ V \\geq T^\\pi V\n",
    "\\end{array} \\right.$$\n",
    "\"For all $\\pi$\" means for all possible association $s\\leftrightarrow a$, so this can be expanded as:\n",
    "$$\\left\\{ \\begin{array}{c}\n",
    "\\min \\sum\\limits_{s\\in S} V(s)\\\\\n",
    "s.t. \\ \\forall (s,a)\\in S\\times A, \\quad V(s) - \\gamma \\sum\\limits_{s'\\in S} p(s'|s,a)V(s') \\geq r(s,a)\n",
    "\\end{array}\\right.$$\n",
    "Which, finally, is a linear program with $|S|$ variables and $|S||A|$ constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"adp\"></a>Asynchronous Dynamic Programming\n",
    "\n",
    "If this is the first time you read this notebook, this part can be skipped. [Link to next part](#rl)\n",
    "\n",
    "We have seen that Value Iteration and Policy Iteration are Dynamic Programming algorithms. They follow a path, respectively in value function and in policy space that leads to $V^*$ and $\\pi^*$. But we can remark that they both perform *state-wise* operations such as:\n",
    "\n",
    "- $V(s) \\leftarrow \\max_{a} r(s,a) + \\gamma \\sum_{s'} p(s'|s,a) V(s')$ for value iteration\n",
    "- $\\pi(s) \\leftarrow \\arg\\max_a Q^{\\pi}(s,a)$ for policy iteration\n",
    "- $V(s) \\leftarrow r(s,\\pi(s)) + \\gamma \\sum_{s'} p(s'|s,\\pi(s)) V(s')$ for policy evaluation\n",
    "\n",
    "These state-wise operations are called Bellman backups. Let's rename them, respectively:\n",
    "- `BBVopt(V,s): return` $\\max_{a} r(s,a) + \\gamma \\sum_{s'} p(s'|s,a) V(s')$\n",
    "- `BBpiopt(V,s): return` $\\arg\\max_a r(s,a) + \\gamma \\sum_{s'} p(s'|s,a) V(s)$\n",
    "- `BBVval(V,s): return` $r(s,\\pi(s)) + \\gamma \\sum_{s'} p(s'|s,\\pi(s)) V(s')$\n",
    "\n",
    "Then we can rewrite value iteration as:<br>\n",
    "`\n",
    "V(s) = Vinit(s) for all s\n",
    "while error>epsilon\n",
    "  for s in S\n",
    "    W(s) = BBVopt(V,s)\n",
    "  error = norm(W-V)\n",
    "  V = W\n",
    "`\n",
    "\n",
    "And policy iteration as:<br>\n",
    "`\n",
    "while(true)\n",
    "  V(s) = 0 for all s\n",
    "  for k=1 to K\n",
    "    for s in S\n",
    "      V(s) = BBVval(V,s)\n",
    "  for s in S\n",
    "    pi2 BBpiopt(V,s)\n",
    "  if(pi == pi2)\n",
    "    stop\n",
    "  else\n",
    "    pi = pi2\n",
    "`\n",
    "\n",
    "### <a id=\"avi\"></a>Asynchronous Value Iteration\n",
    "\n",
    "Let's take the pseudo-code of value iteration. Why don't we perform directly `V(s) = BBVopt(V,s)`, instead of relying on the intermediate `W` function? Doing so is actually called **Gauss-Seidl Value Iteration** and it opens the door to a much wider class of algorithms called Asynchronous Value Iteration.<br>\n",
    "<br>\n",
    "It is crucial to note that in Gauss-Seidl Value Iteration, the order in which the states are considered for backups greatly affects of rewards are propagated through the state space and how the sequence of value functions converges to $V^*$. But still, in Gauss-Seidl Value Iteration, states are updated once per sweep over the state space.\n",
    "\n",
    "Why wouldn't we update the value of some states more often than others? Would the overall value function still converge to $V^*$? A very powerful theorem actually states that:\n",
    "<div class=\"alert alert-success\">As long as every state is visited infinitely often by the `V(s)` $\\leftarrow$ `BBVopt(V,s)` operation as time tends to $+\\infty$, the value function $V$ converges to $V^*$\n",
    "</div>\n",
    "\n",
    "Consequently, we could pick states totally randomly in order to perform Bellman backups on $V$ and $V$ would still converge to $V^*$. Although picking states randomly for that purpose seems like a bad idea, identifying a good ordering for the backups can lead to drastic improvements in convergence speed. This is the key idea of **Asynchronous Value Iteration** and has justified (among other things) the popular **Prioritized Sweeping** and **Real-Time Dynamic Programming** algorithms.\n",
    "\n",
    "### <a id=\"api\"></a>Asynchronous Policy Iteration\n",
    "\n",
    "Let's now take the pseudo-code of policy iteration. The evaluation step and the improvement step are clearly separated. But one can note that if we require the evaluation step to have infinite precision, `K` needs to tend to $\\infty$. So, realistic implementations of Policy Iteration, as the one we crafted earlier, tolerate some error on $V^\\pi$. But could we take an arbitrary value for `K`? Would this still converge? The answer is yes and this algorithm is known as **Modified Policy Iteration**.\n",
    "\n",
    "But then, can we push it a little further and say that `K=1`? Then in this case, Modified Policy Iteration becomes... Value Iteration (check it by yourself)!\n",
    "\n",
    "As in the value iteration case, can we update the value or policy of a given state in any ordering? The answer is yes again and the theorem states that:\n",
    "<div class=\"alert alert-success\"> As long as every state is visited infinitely often by the `V(s)` $\\leftarrow$ `BBVval(V,s)`  and the $\\pi(s) \\leftarrow$ `BBpiopt(V,s)` operations as time tends to $+\\infty$, the value function $V$ and the policy $\\pi$ converge respectively to $V^*$ and $\\pi^*$\n",
    "</div>\n",
    "\n",
    "That is the most general framework one can give for **Asynchronous Dynamic Programming** in MDP resolution. It is often called **Asynchronous Policy Iteration**.<br>\n",
    "<br>\n",
    "<div class=\"alert alert-warning\"><b>Exercice:</b><br>\n",
    "Instead of reasoning on state value functions $V$, we can work directly with state-action value functions $Q$ and policies $\\pi$. Suppose we maintain a memory of a function $Q$ and a policy $\\pi$, then we can define two Bellman backups:\n",
    "<ul>\n",
    "<li> `BBQ(s,a):` $Q(s,a) \\leftarrow r(s,a) + \\gamma \\sum_{s'} p(s'|s,a) Q(s',\\pi(s'))$\n",
    "<li> `BBpi(s):` $\\pi(s) \\leftarrow \\arg\\max_a Q(s,a)$\n",
    "</ul>\n",
    "We can note the similarity between `BBQ` and the previous `BBVval` and the one between `BBpi` and the previous `BBpiopt`.  The previous `BBVopt` is then the successive application of `BBQ` and `BBpi`.<br>\n",
    "Rewrite Value Iteration and Policy Iteration with these two elementary operations. Rephrase the key property of Asynchronous Policy Iteration using these operations.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\"><a href=\"#answers3\" data-toggle=\"collapse\"><b>Answers:</b></a><br>\n",
    "<div id=\"answers3\" class=\"collapse\">\n",
    "Value Iteration:\n",
    "<code>\n",
    "Q(s,a) = Qinit(s,a) for all s,a\n",
    "while error>epsilon\n",
    "  for s,a in SxA\n",
    "    BBQ(s,a)\n",
    "    BBpi(s)</code>\n",
    "<br><br>\n",
    "Policy Iteration:\n",
    "<code>\n",
    "while(pi not constant)\n",
    "  Q(s,a) = 0 for all s,a\n",
    "  while error>epsilon\n",
    "    for s,a in SxA\n",
    "      BBQ(s,a)\n",
    "  for s in S\n",
    "    BBpi(s)</code>\n",
    "<br><br>\n",
    "Asynchronous Policy Iteration: as long as all $s$ and all $a$ are visited infinitely often for Bellman backups `BBQ` and `BBpi` as time tends to $+\\infty$, $Q$ and $\\pi$ will tend respectively to $Q^*$ and $\\pi^*$ (whatever the ordering of backups).\n",
    "</div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id=\"rl\"></a>Model-free Reinforcement Learning\n",
    "\n",
    "So we have seen quite a lot of things about MDPs and how to find optimal policies for them. That would be helpful if we always knew the MDP model underlying the control problem we wish to solve, but in most cases, we have no idea of the transition probabilities or the reward model. If you are playing FrozenLake for the first time for example, you will have to try different actions to get a sense of the transition model. Moreover, this class is called *Reinforcement Learning* and we haven't seen much learning yet.\n",
    "\n",
    "We're coming to it. Let's first put it all in perspective."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"taxonomy\"></a>Taxonomy of Machine Learning\n",
    "\n",
    "You may have had classes on Machine Learning before. There are three strongly distinct categories of problems in ML:\n",
    "- Supervized Learning\n",
    "- Unsupervized Learning\n",
    "- Reinforcement Learning\n",
    "\n",
    "Let's try to answer the following questions for each category.\n",
    "- What's the abstract problem we are trying to solve?\n",
    "- What's the data provided to the algorithms?\n",
    "- Give examples of algorithms in SL/UL/RL.\n",
    "\n",
    "<table border=\"1\">\n",
    "<tr>\n",
    "  <td> \n",
    "      \n",
    "**Question** </td>\n",
    "  <td> **Supervized** </td>\n",
    "  <td> **Unsupervized** </td>\n",
    "  <td> **Reinforcement** </td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td> Target </td>\n",
    "  <td> $f(x)=y$ </td>\n",
    "  <td> $x\\in X$ </td>\n",
    "  <td> $\\pi(s)=a$ </td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td> Target (rephrased) </td>\n",
    "  <td> Predict outputs given inputs</td>\n",
    "  <td> Discover structure in data </td>\n",
    "  <td> Find an optimal behavior </td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td> Data </td>\n",
    "  <td> $\\left\\{\\left(x,y\\right)\\right\\}$ supervisor's labels </td>\n",
    "  <td> $\\left\\{x\\right\\}$ unlabelled data </td>\n",
    "  <td> $\\left\\{\\left(s,a,r,s'\\right)\\right\\}$ experience samples </td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td> Output </td>\n",
    "  <td> Classifier or regressor</td>\n",
    "  <td> Clusters or dimension reduction </td>\n",
    "  <td> Policies, value functions </td>\n",
    "</tr>\n",
    "<tr>\n",
    "  <td> Algorithms </td>\n",
    "  <td> Neural networks, SVMs</td>\n",
    "  <td> k-means, PCA </td>\n",
    "  <td> Q-learning </td>\n",
    "</tr>\n",
    "</table>\n",
    "\n",
    "This table helps show the different nature of the problems tackled. The RL problem is about finding the optimal policy for a given environment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a id=\"introrl\"></a>Introduction to model-free Reinforcement Learning\n",
    "\n",
    "Reinforcement Learning studies the optimization of an agent's behaviour when interacting with an a-priori unknown environment. This general framework can be depicted as:\n",
    "\n",
    "<img src=\"images/agent4-1.png\">\n",
    "\n",
    "So the problem RL tries to solve is the *evaluation* and the *improvement* of the agent's behavior, based on experience samples:\n",
    "$$(s,a,r,s')$$\n",
    "\n",
    "One can distinguish two subproblems in RL:\n",
    "- Value prediction: what is this policy's value?<br>\n",
    "$\\rightarrow$ Useful for decision support applications: predict the cost of using the water in a hydro-electric reservoir, predict the expected gain of an investment policy.\n",
    "- Policy optimization: what is the best control policy for this system?<br>\n",
    "$\\rightarrow$ Useful for control applications: robotic actuator control, operations planning, elevator scheduling, agro-ecosystems management, etc.\n",
    "\n",
    "The core hypothesis of RL, is that the environment behaves as an MDP, although we do not know its transition and reward function.\n",
    "\n",
    "Finally, there are different contexts in RL that we can discuss via a few keywords:\n",
    "- curse of dimensionnality<br>\n",
    "The more state/actions variables, the larger the state space and the harder it is to efficiently optimize the policy.\n",
    "- exploration/exploitation<br>\n",
    "Suppose we have discovered some region of the state space yields great rewards. Should we rush towards that region or explore the unknown in search for even better rewards?\n",
    "- model-based vs. model-free RL\n",
    "One straightforward approach to RL would be to collect samples $(s,a,r,s')$ and approximate the underlying MDP in order to use the algorithms seen above. Model-free RL tries to infer the optimal policy without this intermediate step.\n",
    "- online vs. interactive vs. non-interactive RL<br>\n",
    "Online: we directly interact with the environment, as a player in a game.<br>\n",
    "Interactive: we control the interaction with the environment, in particular the choice of the current state or the possility of resets, as if we had a simulator that we can query.<br>\n",
    "Offline: somebody collected that $\\left\\{\\left(s,a,r,s'\\right)\\right\\}$ pile of data and left it here without a note. Can we still get something out of it?\n",
    "- on-policy vs. off-policy<br>\n",
    "Can we learn $V^\\pi$ while applying $\\pi'$?"
   ]
  },
  {
   "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.6.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "426.717px",
    "left": "1354.8px",
    "top": "159.567px",
    "width": "486px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}