<jupyter_start><jupyter_text>Tabular Q Learning on Frozen Lake game Import required libraries<jupyter_code>'''
These libraries are related with OpenAI Gym and Frozen Lake enviornment.
'''
import gym
from gym.envs.registration import register
from gym.envs.toy_text.frozen_lake import LEFT, RIGHT, DOWN, UP
'''
These libraries are related with algebraic calculations and plotting. 
'''
%matplotlib inline
import seaborn
import matplotlib.pyplot as plt
import numpy as np
import numpy.linalg as LA
import time
from tabulate import tabulate<jupyter_output><empty_output><jupyter_text>Initialize the Froze Lake env  Code Parameters * Size of the world: It is a grid world with `4 x 4` grid size.* The discount factor, $\gamma$ has been set to 0.9.* The environment is slippery, i.e., the transition kernel is stochastic<jupyter_code>register(
    id='D4x4-FrozenLake-v0',
    entry_point='gym.envs.toy_text.frozen_lake:FrozenLakeEnv',
    kwargs={'map_name': '4x4',
            'is_slippery': True})
env = gym.make('D4x4-FrozenLake-v0')
gamma=0.9
action_names = {LEFT: 'LEFT', RIGHT: 'RIGHT', DOWN: 'DOWN', UP: 'UP'}<jupyter_output><empty_output><jupyter_text>The Environment* The environment consists of 16 states.     * `env.nS`: count of state variable of data type int.* The agent can take 4 actions.     * `env.nA`: count of action variable of data type int.* The transition kernel `P` is a dictionary. * `P[state][action]` is a tuple.     * `P[state][action] = [probability, nextstate, reward, terminal]`.* For more information check here [link](https://gym.openai.com/envs/FrozenLake-v0/).<jupyter_code>'''
Print the total no. of states and actions in this env.
'''
print('Number of Actions', env.nA) 
print('Number of States ', env.nS)
print('P[5,0]', env.P[5][0])<jupyter_output>Number of Actions 4
Number of States  16
P[5,0] [(1.0, 5, 0, True)]<jupyter_text>Generating a Heatmap * Pass the optimal value function and policy to the function `fancy_visual` to obtain a heat map.* This function also prints the value function and policy. Hence you do not have to print it.<jupyter_code>def fancy_visual(value_func,policy_int):    
    grid = 4
    f, ax = plt.subplots(figsize=(11, 9))
    cmap = seaborn.diverging_palette(220, 10, as_cmap=True)
    reshaped=np.reshape(value_func,(grid,grid))
    seaborn.heatmap(reshaped, cmap=cmap, vmax=1.1,
                square=True, xticklabels=grid+1, yticklabels=grid+1,
                linewidths=.5, cbar_kws={"shrink": .5}, ax=ax, annot=True, fmt="f")
    counter=0
    for j in range(0, 4):
        for i in range(0, 4):
            if policy_int[counter]==1:
                plt.text(i+0.5, j+0.7, u'\u2193', fontsize=12)
            elif policy_int[counter]==3:
                plt.text(i+0.5, j+0.7, u'\u2191', fontsize=12)
            elif policy_int[counter]==0:
                plt.text(i+0.5, j+0.7, u'\u2190', fontsize=12)
            else:
                plt.text(i+0.5, j+0.7, u'\u2192', fontsize=12)
            counter=counter+1

    plt.title('Heatmap of policy iteration with value function values and directions')
    print('Value Function',value_func)
    print('Policy',policy_int)
    plt.show()<jupyter_output><empty_output><jupyter_text>Part 1 Value IterationWe do Value Iteration to get the optimal Q-values, which will be used to compare with the Q-values that we get from Tabular Q-Learning Algorithm. Define important functionsFunctions have been written based on the pseudo code provided in the Sutton-Barto Chapter 4 and Section 4.4 Value Iteration. All the equations have been taken from that textbook.* `value_function_to_policy`: Function which computes the optimal policy in greedy manner based on the optimal input value function.* `value_function_to_Qvalues`: Function which computes the optimal Q-values based on the optimal input value function.* `value_iteration`: Function for doing the value iteration. * `print_policy`: Helper function to print the optimal policy followed by the agent. * `run_policy`: Function for running our trained policy on the game see how it performs.  Definitions* `value_func_sa`: value function for given state and action taken by following the policy.* `max_value_func_sa`: value function for given state and action taken by following the policy such that it is maximized.* `a`: current action* `max_a`: action corresponding to the maximum value function.* `q_values`: Q-values* `delta`: parameter to measure the convergence. Refer the Sutton-Barto textbook.<jupyter_code>def value_function_to_policy(env, gamma, value_function):
    """
    Output action numbers for each state in value_function.
    Parameters
    ----------
    env: gym.core.Environment
      Environment to compute policy for. Must have nS, nA, and P as
      attributes.
    gamma: float
      Discount factor. Number in range [0, 1)
    value_function: np.ndarray
      Value of each state.
    Returns
    -------
    np.ndarray
      An array of integers. Each integer is the optimal action to take
      in that state according to the environment dynamics and the
      given value function.
    """
    policy = np.zeros(env.nS, dtype='int')
    for s in range(env.nS):
        max_value_func_sa = -1
        max_a = -1
        for a in range(env.nA):
            value_func_sa = 0
            for possible_next_state in env.P[s][a]:
                prob_action = possible_next_state[0]
                cur_reward = possible_next_state[2]
                future_reward = gamma * value_function[possible_next_state[1]]
                value_func_sa += prob_action * (cur_reward + future_reward)
            if value_func_sa > max_value_func_sa:
                max_value_func_sa = value_func_sa
                max_a = a
        policy[s] = max_a
    return policy

def value_function_to_Qvalues(env, gamma, value_function):
    """
    Output optimal Q value for each state action pair in env.
    Parameters
    ----------
    env: gym.core.Environment
      Environment to compute policy for. Must have nS, nA, and P as
      attributes.
    gamma: float
      Discount factor. Number in range [0, 1)
    value_function: np.ndarray
      Value of each state.
    Returns
    -------
    np.ndarray
      An matirx (2D array) of optimal Q-values.
    """
    q_values = np.zeros((env.nS, env.nA))
    for s in range(env.nS):
        for a in range(env.nA):
            q_values[s][a] = 0
            for possible_next_state in env.P[s][a]:
                prob_action = possible_next_state[0]
                cur_reward = possible_next_state[2]
                future_reward = gamma * value_function[possible_next_state[1]]
                q_values[s][a] += prob_action * (cur_reward + future_reward)
    return q_values

def value_iteration(env, gamma, max_iterations=int(1e3), tol=1e-10):
    """
    Runs value iteration for a given gamma and environment.
    Parameters
    ----------
    env: gym.core.Environment
      The environment to compute value iteration for. Must have nS,
      nA, and P as attributes.
    gamma: float
      Discount factor, must be in range [0, 1)
    max_iterations: int
      The maximum number of iterations to run before stopping.
    tol: float
      Determines when value function has converged.
    Returns
    -------
    np.ndarray, iteration
      The value function and the number of iterations it took to converge.
    """
    value_func_old = np.random.rand(env.nS)
    value_func_new = np.zeros(env.nS)
    diff_arr = []
    for iteration in range(max_iterations):
        delta = 0
        for s in range(env.nS):
            max_value_func_sa = -1
            for a in range(env.nA):
                value_func_sa = 0
                for possible_next_state in env.P[s][a]:
                    prob_action = possible_next_state[0]
                    cur_reward=possible_next_state[2]
                    future_reward = gamma * value_func_old[possible_next_state[1]]
                    value_func_sa += prob_action * (cur_reward + future_reward)
                if value_func_sa > max_value_func_sa:
                    max_value_func_sa = value_func_sa
            diff = np.abs(value_func_old[s] - max_value_func_sa)
            delta = max(delta, diff)
            value_func_new[s] = max_value_func_sa
        if delta <= tol: 
            break
        sq_diff = LA.norm(value_func_new - value_func_old)
        diff_arr.append(sq_diff)
        value_func_old = value_func_new
    return value_func_new, iteration, diff_arr


def print_policy(policy, action_names):
    """
    Print the policy in human-readable format.
    Parameters
    ----------
    policy: np.ndarray
      Array of state to action number mappings
    action_names: dict
      Mapping of action numbers to characters representing the action.
    """
    str_policy = policy.astype('str')
    for action_num, action_name in action_names.items():
        np.place(str_policy, policy == action_num, action_name)
    print(str_policy)
    return str_policy

def run_policy(env, gamma, policy):
    initial_state = env.reset()
    total_reward = 0
    num_steps = 0
    current_state = initial_state
    while True:
        nextstate, reward, is_terminal, debug_info = env.step(policy[current_state])
        total_reward += np.power(gamma, num_steps) * reward
        num_steps += 1
        if is_terminal:
            break
        current_state = nextstate
    return total_reward, num_steps<jupyter_output><empty_output><jupyter_text>Run Value Iteration* Start the Value iteration and time it to measure the performance.<jupyter_code>print("Executing Value Iteration")
start_time=time.time()
value_function, value_iters, error = value_iteration(env, gamma)
print("Total time taken: " + str((time.time()-start_time)))
print("Total Value Iteration Steps: "+str(value_iters))<jupyter_output>Executing Value Iteration
Total time taken: 0.08605813980102539
Total Value Iteration Steps: 197<jupyter_text>Rate of convergence* We convergence performance of the value iteration by measuring following qunatity $\lVert V_{k+1} - V_{k} \rVert_{2}$ and see how it changes with total number of iterations.* Tolerance which is the termination criteria for the code tells us how close should $V_{k+1}$ and  $V_{k}$ should be. By default the tolerance value has been set to 1e-10 i.e, we want the values to be very tight, and we achieve this tolerance in less than 200 iterations.<jupyter_code>plt.figure(figsize=(8, 6))
plt.plot(np.arange(value_iters), error)
plt.xlabel(r'Number of iterations $\rightarrow$')
plt.ylabel(r'$\|V_{k+1} - V_{k}\|_{2} \rightarrow$')
plt.title('Convergence of the value iteration with respect to iterations')
plt.show()<jupyter_output><empty_output><jupyter_text>Optimal Value function and Optimal Q-values<jupyter_code>print('Optimal value function is:')
print(value_function)
print('\n')
q_val = value_function_to_Qvalues(env, gamma, value_function)
print('Optimal Q-values are:')
print(q_val)<jupyter_output>Optimal value function is:
[6.88909050e-02 6.14145717e-02 7.44097623e-02 5.58073218e-02
 9.18545400e-02 4.19797313e-10 1.12208207e-01 8.46112458e-10
 1.45436355e-01 2.47496955e-01 2.99617593e-01 2.87889594e-10
 4.62412369e-10 3.79935901e-01 6.39020148e-01 1.90230436e-10]


Optimal Q-values are:
[[6.88909050e-02 6.66480050e-02 6.66480050e-02 5.97589145e-02]
 [3.90916432e-02 4.29902003e-02 4.07473003e-02 6.14145717e-02]
 [7.44097623e-02 6.88290301e-02 7.27275873e-02 5.74894968e-02]
 [3.90651255e-02 3.90651255e-02 3.34843933e-02 5.58073218e-02]
 [9.18545400e-02 7.11872686e-02 6.42981781e-02 4.82236336e-02]
 [3.77817582e-10 3.77817582e-10 3.77817582e-10 3.77817582e-10]
 [1.12208207e-01 8.98852783e-02 1.12208207e-01 2.23229291e-02]
 [7.61501212e-10 7.61501212e-10 7.61501212e-10 7.61501212e-10]
 [7.11872686e-02 1.17879993e-01 1.01805449e-01 1.45436355e-01]
 [1.57611677e-01 2.47496955e-01 2.03866048e-01 1.33516185e-01]
 [2.99617593e-01 2.65955131e-01 2.25368507e-01 1.07911549e-01]
 [2.5910063[...]<jupyter_text>Optimal Policy* Print the optimal policy which tells the agents what action take based on the current position of the game and optimal value function. * We will aslo plot the output in fancy manner to get more better idea.<jupyter_code># Print the policy
print("Policy:") 
policy = value_function_to_policy(env, gamma, value_function)
policy_str = print_policy(policy, action_names)
# Plot the policy
fancy_visual(value_function, policy)
# Plot
ps=[]
for elem in policy_str:
    ps.append(elem[0])
reshaped_policy=np.reshape(ps,(4,4))
# print(tabulate(reshaped_policy,tablefmt='latex'))
print('Above policy can also be expressed in a tabular form as follows:')
print(tabulate(reshaped_policy))<jupyter_output>Above policy can also be expressed in a tabular form as follows:
-  -  -  -
L  U  L  U
L  L  R  L
U  D  L  L
L  R  D  L
-  -  -  -<jupyter_text>Performance measure<jupyter_code>total_cum_reward=0
maxn = 5
start_time=time.time()
for n in range(maxn):
    cum_reward, nsteps = run_policy(env, gamma, policy)
    total_cum_reward += cum_reward
    if n%1 == 0: 
        print("Done " + str(n)) 
print("Time: " + str((time.time()-start_time)/60))
print("No. of steps for the last run: " + str(nsteps)) 
print("Performance of the agent over " + str(maxn) + " episodes")
print("Average Cumulative Reward: " + str((total_cum_reward/maxn)))<jupyter_output>Done 0
Done 1
Done 2
Done 3
Done 4
Time: 0.0001902341842651367
No. of steps for the last run: 99
Performance of the agent over 5 episodes
Average Cumulative Reward: 0.07826099076725651<jupyter_text>Part 2 Tabular Q-Learning Hyperparameters* In this part I have implemented the Q-learning algorithm. This code has been written based on the pseudo code given the Sutton-Barto textbook in chapter 6 and section 6.5. This version of the tabular Q-learning is known as the Off-policy TD Control Q-learning.* Below I have defined some hyperparameters based on the algorithm and these required for learning the optimal Q-values.<jupyter_code>total_episodes = 150000        # Total episodes
learning_rate = 0.8            # Learning rate
gamma = 0.95                   # Discounting rate

# Exploration parameters
epsilon = 1.0                 # Exploration rate
max_epsilon = 1.0             # Exploration probability at start
min_epsilon = 0.01            # Minimum exploration probability 
decay_rate = 0.005            # Exponential decay rate for exploration prob

# Q-value table
qtable = np.zeros((env.nS, env.nA))
# qtable = np.random.rand(env.nS, env.nA)
# qtable[15:] = 0

# Book keeping parameters
# List of rewards
rewards = []
episode_counter = 0
sq_diff_qval = []

# Debug
verbose = False<jupyter_output><empty_output><jupyter_text>The Q-Learning Alogrithm* Here is the main loop for doing the tabular Q-learning. Note that we during learning we decay our learning rate and epsilon greedy parameter after each episode so that satisfy the Robbins–Monro conditions.<jupyter_code># For total no. of episodes or until learning is stopped
for episode in range(total_episodes):
    # Reset the environment
    state = env.reset()
    is_terminal = False
    total_rewards = 0
    episode_counter += 1
    
    # We break the epsiode and go to next one as soon as we die or reach the goal
    while not is_terminal:
        # Choose an action A in the current world state S
        # First we randomize a number, which help in inducing the 
        # exploration and exploitation effect
        exp_exp_tradeoff = np.random.uniform(0, 1)
        
        # If this number > greater than epsilon 
        # then exploitation (taking the biggest Q value for this state)
        if exp_exp_tradeoff > epsilon:
            action = np.argmax(qtable[state,:])

        # Else doing a random choice which means exploration
        else:
            action = env.action_space.sample() # random action

        # Take the action A and observe the outcome state S' and reward R
        new_state, reward, is_terminal, info = env.step(action) # take a step

        # Update Q(S,A):= Q(S,A) + eta * [R(S,A) + gamma * max{Q(S',A')} - Q(S,A)]
        # qtable[new_state,:] : all the actions we can take from new state
        qtable[state, action] = qtable[state, action] + learning_rate * (reward + gamma * np.max(qtable[new_state, :]) 
                                                                         - qtable[state, action])
        
        total_rewards += reward
        
        # Our new state is state
        state = new_state
    
    # Do some book keeping: store the values of reward and difference between optimal q-value
    # that we obtained from the dynamic programming and q-values we got from tabular q-learning
    # Reduce epsilon because we need less and less exploration and 
    # reduce the learning rate
    epsilon = min_epsilon + (max_epsilon - min_epsilon) * np.exp(-decay_rate * episode) 
    learning_rate = learning_rate * np.exp(-0.00005)
    rewards.append(total_rewards)
    sq_diff_qval.append(LA.norm(qtable - q_val))
    
    # Debugging
    if verbose:
        if episode_counter%1000 == 0:
            print('In episode: ', episode_counter, ' and delta: ', LA.norm(qtable-q_val))
print ("Score over time: " +  str(sum(rewards)/total_episodes))<jupyter_output>Score over time: 0.6799<jupyter_text>Optimal Q-values and Optimal Value function<jupyter_code>print('Optimal Q-values and Optimal value function obtained from Tabular Q-learning algoritm')
value_function = np.max(qtable, axis=1)
print('Optimal value function is:')
print(value_function)
print('\n')
print('Optimal Q-values are:')
print(qtable)<jupyter_output>Optimal Q-values and Optimal value function obtained from Tabular Q-learning algoritm
Optimal value function is:
[0.18150054 0.15530499 0.15446939 0.12390759 0.20928477 0.
 0.17694596 0.         0.26999698 0.37248143 0.39958396 0.
 0.         0.50904335 0.72518004 0.        ]


Optimal Q-values are:
[[0.18150054 0.17191766 0.17251848 0.16328259]
 [0.09560385 0.09666919 0.10011551 0.15530499]
 [0.15446939 0.13183677 0.1331532  0.12670612]
 [0.05071176 0.03421913 0.04270463 0.12390759]
 [0.20928477 0.15732286 0.14271295 0.12187214]
 [0.         0.         0.         0.        ]
 [0.12920402 0.11321215 0.17694596 0.04646206]
 [0.         0.         0.         0.        ]
 [0.16205737 0.20395685 0.18533047 0.26999698]
 [0.24233804 0.37248143 0.29521572 0.22029204]
 [0.39958396 0.36394531 0.28458796 0.1820341 ]
 [0.         0.         0.         0.        ]
 [0.         0.         0.         0.        ]
 [0.2898651  0.38452165 0.50904335 0.35035077]
 [0.51117287 0.72518004 0.68191074 0.6240[...]<jupyter_text>Optimal Policy* Print the optimal policy which tells the agents what action take based on the current position of the game and optimal value function. Note that optimal policy can be computed directly from the Q-values as shown below. These values are computed from the Q-values we got from Tabular Q-learning. * We will aslo plot the output in fancy manner to get more better idea.<jupyter_code># Compute the policy directly from the Q-values
policy = np.argmax(qtable, axis=1)

# Print the policy
print("Policy:") 
policy = value_function_to_policy(env, gamma, value_function)
policy_str = print_policy(policy, action_names)
# Plot the policy
fancy_visual(value_function, policy)
# Plot
ps=[]
for elem in policy_str:
    ps.append(elem[0])
reshaped_policy=np.reshape(ps,(4,4))
print('Above policy can also be expressed in a tabular form as follows:')
print(tabulate(reshaped_policy))<jupyter_output>Above policy can also be expressed in a tabular form as follows:
-  -  -  -
L  U  L  U
L  L  L  L
U  D  L  L
L  R  D  L
-  -  -  -<jupyter_text>Convergence* Here we observe how our Q-values obtained from Tabular Q-learning algorithm converges to the optimal Q-values obtained through dynamic programming. We see this trend over mutliple episodes.* We also do a moving window average to remove the noisy zig-zag patterns and get a smooth trend for understanding.<jupyter_code>plt.figure(figsize=(8, 6))
plt.plot(np.arange(episode_counter), sq_diff_qval)
plt.xlabel(r'No. of episodes $\rightarrow$')
plt.ylabel(r'$\vert\vert Q_{k} - Q^{*}\vert\vert_{2} \rightarrow$')
plt.title('Convergence of the Tabular Q-learning to optimal value over episodes')
plt.show()
def moving_average(a, n):
    '''
    Function for doing moving window average
    '''
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n-1:]/n
# Moving avg window has been set to 10000
sq_diff_avg = moving_average(np.array(sq_diff_qval), 10000)
plt.figure(figsize=(8, 6))
plt.plot(np.arange(len(sq_diff_avg)), sq_diff_avg)
plt.xlabel(r'No. of episodes $\rightarrow$')
plt.ylabel(r'$\vert\vert Q_{k} - Q^{*}\vert\vert_{2} \rightarrow$')
plt.title('Moving average convergence of the Tabular Q-learning to optimal value over episodes')
plt.show()<jupyter_output><empty_output><jupyter_text>Cumulative Rewards* We see how did the rewards improved over episodes of learning. Note since output of each output will be 0 or 1 so to get a more proper understanding of the improvement we take cummulative sum of the rewards over episodes and then normalize them with total no. of episodes. * We also moving window average over the rewards and cumulative sum of rewards to see the trend.<jupyter_code>total_cum_reward = np.cumsum(np.array(rewards)/total_episodes)
plt.figure(figsize=(8, 6))
plt.plot(np.arange(total_cum_reward.shape[0]), total_cum_reward)
plt.xlabel(r'No. of episodes $\rightarrow$')
plt.ylabel(r'Cumulative sum of the rewards $\rightarrow$')
plt.title('Improvement of rewards over episodes')
plt.show()
# Moving avg window has been set to 10000
avg_reward = moving_average(total_cum_reward, 10000)
plt.figure(figsize=(8, 6))
plt.plot(np.arange(len(avg_reward)), avg_reward)
plt.xlabel(r'No. of episodes $\rightarrow$')
plt.ylabel(r'Cumulative sum of the rewards $\rightarrow$')
plt.title('Moving average improvement of rewards over episodes')
plt.show()
# Note: we are only considering rewards not the cumulative sum
# Moving avg window has been set to 1000
avg_reward = moving_average(np.array(rewards), 1000)
plt.figure(figsize=(8, 6))
plt.plot(np.arange(len(avg_reward)), avg_reward)
plt.xlabel(r'No. of episodes $\rightarrow$')
plt.ylabel(r'Rewards $\rightarrow$')
plt.title('Moving average improvement of rewards over episodes')
plt.show()<jupyter_output><empty_output><jupyter_text>Performance measure<jupyter_code>total_cum_reward=0
maxn = 5
start_time=time.time()
for n in range(maxn):
    cum_reward, nsteps = run_policy(env, gamma, policy)
    total_cum_reward += cum_reward
    if n%1 == 0: 
        print("Done " + str(n)) 
print("Time: " + str((time.time()-start_time)/60))
print("No. of steps for the last run: " + str(nsteps)) 
print("Performance of the agent over " + str(maxn) + " episodes")
print("Average Cumulative Reward: " + str((total_cum_reward/maxn)))<jupyter_output>Done 0
Done 1
Done 2
Done 3
Done 4
Time: 0.0001345078150431315
No. of steps for the last run: 13
Performance of the agent over 5 episodes
Average Cumulative Reward: 0.3108728782061132
StarcoderdataJupyterStructuredCleanDedup