<jupyter_start><jupyter_text>S&B 2.5 I. Single Example<jupyter_code># Initialization
q_star_i = 0 # can be modified
k = 10 # Number of levers for the bandit
steps = 10000
# The cumulative drift in each q* due to the random walk
# Indices along axis 0 correspond to action, a = 0,...,9
cum_walk = np.cumsum(0.01 * np.random.randn(k, steps), axis=1)
q_star = q_star_i + cum_walk
# Plot the drift of the q*'s to get a sense of expectations
plt.figure(figsize=(20,15))
plt.plot(q_star.T)
plt.title("Cumulative drift of q*'s", size=40)
plt.xlabel("Step Number", size=40)
plt.xticks(size=30)
plt.ylabel("q* Value", size=40)
plt.yticks(size=30)
plt.legend(range(k), prop={'size':30})
plt.tight_layout()
plt.show()
# This represents the samples actually draw by the agent
# Summing q* here sets the mean to q*(action, step)
bandit_samples = q_star + np.random.randn(k, steps)
class Bandit():
    
    def __init__(self, k, epsilon):
        self.estimates = np.zeros(k)
        self.step_count = 0
        self.k = k
        self.epsilon = epsilon
        
    def sample_policy(self):
        is_exploring = np.random.binomial(1, self.epsilon)
        
        if not is_exploring:
            max_Q = np.amax(self.estimates)
            greedy_actions = np.argwhere(max_Q == self.estimates).flatten()
            next_action = np.random.choice(greedy_actions)
            
        else:
            next_action = np.random.randint(0, self.k)
        
        return next_action
        
    def update_estimates(self, lever_num, sample):
        raise NotImplementedError
        
class AverageBandit(Bandit):
    
    def update_estimates(self, lever_num, sample):
        self.step_count += 1
        self.estimates[lever_num] = self.estimates[lever_num] + \
            + (1/self.step_count) * (sample - self.estimates[lever_num])
        
class ConstantStepBandit(Bandit):
    
    def __init__(self, k, epsilon, alpha):
        self.alpha = alpha
        Bandit.__init__(self, k, epsilon)
    
    def update_estimates(self, lever_num, sample):
        self.estimates[lever_num] = self.estimates[lever_num] + \
            + (self.alpha) * (sample - self.estimates[lever_num])
        
epsilon = 0.1
alpha = 0.1
average_bandit = AverageBandit(k, epsilon)
step_bandit = ConstantStepBandit(k, epsilon, alpha)
bandits = [average_bandit, step_bandit]
reward_curves = np.empty((2, steps))
action_taken = np.empty((2, steps))

for step in range(steps):
    if step % 500 == 0:
        print("{} steps have been executed...".format(step))
    
    for bandit_index, bandit in enumerate(bandits):
        next_action = bandit.sample_policy()
        reward = bandit_samples[next_action, step]
        bandit.update_estimates(next_action, reward)
        reward_curves[bandit_index, step] = reward
        action_taken[bandit_index, step] = next_action
optimal_actions = np.argmax(q_star, axis=0)
optimal_action_taken = (optimal_actions == action_taken)
cumulative_rewards = np.cumsum(reward_curves, axis=1)
plt.plot(cumulative_rewards.T)<jupyter_output><empty_output><jupyter_text>II. Multiple Simulations<jupyter_code>def generate_bandit_samples(q_star_i, k ,steps):
    cum_walk = np.cumsum(0.01 * np.random.randn(k, steps), axis=1)
    q_star = q_star_i + cum_walk
    bandit_samples = q_star + np.random.randn(k, steps)
    
    return q_star, bandit_samples
def run_simulations(epsilon, alpha, simulations, q_star_i, k ,steps):
    
    reward_curves = np.empty((2, simulations, steps))
    action_taken = np.empty((2, simulations, steps))
    q_star_sims = np.empty((simulations, k, steps))
    t_avg = 0
    
    for sim in range(simulations):
        t_i = time.time()
        
        average_bandit = AverageBandit(k, epsilon)
        step_bandit = ConstantStepBandit(k, epsilon, alpha)
        bandits = [average_bandit, step_bandit]
        q_star, bandit_samples = generate_bandit_samples(q_star_i, k, 
                                                         steps)
        q_star_sims[sim] = q_star

        for step in range(steps):

            for bandit_index, bandit in enumerate(bandits):
                next_action = bandit.sample_policy()
                reward = bandit_samples[next_action, step]
                bandit.update_estimates(next_action, reward)
                reward_curves[bandit_index, sim, step] = reward
                action_taken[bandit_index, sim, step] = next_action
                
        t_avg = (time.time() - t_i)/(sim + 1) + sim/(sim+1) * t_avg 
        t_left = (simulations - sim - 1) * t_avg
        
        if (sim + 1) % 100 == 0:
            print("\nSimulation #{}".format(sim + 1))
            print("Avg. Sim. Time: {} s".format(round(t_avg, 2)))
            print("Expected Time Remaining: {} mins".format(round(t_left/60, 2)))
        
            
    optimal_actions = np.argmax(q_star_sims, axis=1)
    optimal_action_perc = np.average(optimal_actions == action_taken,
                                     axis=1)
    avg_rewards = np.average(reward_curves, axis=1)
    
    return optimal_action_perc, avg_rewards
# Could be implemented more efficiently if simulations vectorized
# This implementation is enough for our purposes as it is though
optimal_action_perc, avg_rewards = \
    run_simulations(epsilon, alpha, 2000, q_star_i, k, steps)
plt.figure(figsize=(20,15))
plt.plot(optimal_action_perc.T)
plt.title("% Optimal Action Selected, Average Over 2,000 Sims", size=40)
plt.xlabel("Step Number", size=40)
plt.xticks(size=30)
plt.ylabel("% Optimal Action", size=40)
plt.yticks(size=30)
plt.legend(range(k), prop={'size':30})
plt.tight_layout()
plt.show()
plt.figure(figsize=(20,15))
plt.plot(avg_rewards.T)
plt.title("Average Rewards, Average Over 2,000 Sims", size=40)
plt.xlabel("Step Number", size=40)
plt.xticks(size=30)
plt.ylabel("Average Reward", size=40)
plt.yticks(size=30)
plt.legend(range(k), prop={'size':30})
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>So we can see that the constant step averaging performs significantly better since it can capture the non-stationarity of the problem.   S&B  4.7 Jack's Car Rental ModifiedLet's let $V(c_1, c_2)$ be the value function, which depends on the number of cars at the first and second location respectively. Thus $s = (c_1, c_2)$, but it is more convenient to work with the individual car counts.There are a few constraints we need to enforce:+ We can never have more than 20 cars at any location, so we need to sum the tails of the Poisson distributions for returns.+ We can never have less than 0 cars at any location, so we need to sum the tails for the Poisson distributions for the rentals.+ The policy for any pair state $s=(c_1, c_2)$ can only take on certain values. These are determined by the fact that we can move a maximum of $b(=5)$ cars overnight and the fact that we can be no more than $t_1, t_2(=20)$ cars at each location:$$\mathscr{A}(c_1, c_2) = [-min(t_1 - c_1, min(c_2, b)), min(t_2 - c_2, min(c_1, b))] \cap \mathbb{Z}$$Let $req_1$, $req_2$ and $ret_1$, $ret_2$ be the number of cars requested and returned respectively for each of the locations. Let $t_1$, $t_2$ denote the top number of cars we can hold at each location$. Then note that the maximum number of cars that can be rented each day is equal to the number of cars after transfers are accounted for:$$ \begin{align*}maxrent_1 &= C_1 - \pi(C_1, C_2) \\maxrent_2 &= C_2 + \pi(C_1, C_2)\end{align*}$$Here $C_i$ is the state today. And the state the next day is:$$ \begin{align*}C_1' &= min[max[C_1 - Req_1 - \pi(C_1, C_2), \ 0] + Ret_1, \ t_1] \\C_2' &= min[max[C_2 - Req_2 + \pi(C_1, C_2), \ 0] + Ret_2, \ t_2]\end{align*}$$Here $Ret_i$, $Req_i$ are the cars requested and returned tomorrow day, and $C_i'$ is the state tomorrow. We use the capital letters to denote that these are random variables rather than specific values. The max and min functions guarantee that more than the number of cars available are not rented, and that we do not hold more than $t_1, t_2$ cars per location, respectively.Finally, if we let $R_i$ denote the reward we made at each location, $P$ denotes the cost of moving the cars, $p$ the price of moving a car between locations, and $m$ denotes the amount of money we make per rental then we have that:$$ \begin{align*}R_1 &= m \cdot min[Req_1, \ maxrent_1(C_1, C_2)] \\R_2 &= m \cdot min[Req_2, \ maxrent_2(C_1, C_2)] \\ P &= p \cdot \left| \pi(C_1, C_2) \right|\end{align*}$$And the total reward is just $R = R_1 + R_2 - P$. Since our problem is defined in terms of the distributions for $Req_i$ and $Ret_i$ it's easier to write the sum required for the update rule in terms of these variables. This is why we have derived these equations. We can use them now to rewrite the update rule -- note that every term in the expectation is a function of $C_1, C_2, ret_1, ret_2, req_1, ret_2$, so since the former 2 variables are fixed it suffices to take the expectation over the latter four:$$ \begin{align*}V(s) &\leftarrow q_{\pi}(s, \pi(s)) = E_{Ret_1, Ret_2, Req_1, Req_2}[R + \gamma V(C_1', C_2') \ | \ S=(C_1,C_2), A=\pi(C_1, C_2)]\end{align*}$$Here $V$ is our current estimate of $v_{\pi}$, which is a lookup table of values. The expectation is taken over the distribution:$$ p(req_1, req_2, ret_1, ret_2) = p(req_1) \cdot p(req_2) \cdot p(ret_1) \cdot p(ret_2)$$Where each of the indepedent distributions is a Poisson distribution with means $\lambda_{req_1}$, $\lambda_{req_2}$, $\lambda_{ret_1}$, and $\lambda_{ret_2}$ respectively.Given this formula for the expectation we can calculate it in many ways:+ Analytically computing the expectation for each term and subterm within the expectation, to derive a closed formula for the expectation. If we did this we may be able to sum the tail of the probability distributions as aforementioned. However this analytical computation is likely complex (in particular for $C_1', C_2'$) and not worth the increased precision.+ Calculating the value of the probability and of the variables, and calculating the infinite sum in the expectation until the tail of the sum gets sufficiently small. This is not a bad approach but it raises the question of how to order the terms in the sum so that the terms are monotonically decreasing in size.+ Sampling the probability distribution in order to produce a Monte Carlo estimate of $V(s)$. Note that even though we are estimating the expecation with MC in this case, this is NOT the same as the Reinforcement Learning MC approach of sampling trajectories. This would still be a Dynamic Programming approach, but with MC estimation of the environment distribution as a subroutine.The simplest implementation seemed to be the Monte Carlo based one so we attempted this alternative. Unfortunately we did not achieve sufficient accuracy with this approach for the policy to converge. It appears that even with $10^6$ samples we get $O(1)$ error in the policy matrix which means the policy does not converge. Fortunately we overestimated the complexity of motivating the sum truncation approach analytically. We note that each of the Poisson distributions is independent and that if we sum all of the terms for each random variable up to the $N$th term, the the largest the error can be is bounded by $O(1/N!)$. The factorial grows superexponentially, so even at $N=30$ we get an error less than $(10^{-33})$ and we only have to sum less than $10^6$ sum tertms. Thus we attempt this approach next.Before moving on to the implementation, we develop the theory behind the modifications required for Exercise 4.7. There are two:+ An employee is willing to move one car from location 1 to location 2 for free.+ If more than 10 cars are kept at either location an extra $\$4$ dollars must be paid to use an extra parking lot, at  each (for a poissible total of $\$8$).Note that both of these are simply modifications to the cost term $P$ of the reward function. Thus to apply the modifications we can simply define a new cost term:$$ P_{mod} = p \cdot \left| \pi(C_1, C_2) \right| - p \cdot H(\pi(C_1, C_2) - 1) + $_{park, 1} \cdot H(maxrent_1 - 10) + $_{park,2} \cdot H(maxrent_2 - 10) $$Here $H$ is the Heaviside step function with the convention that $H(0) = 1$ and $\$_{park, i}(=4)$ is the parking lot cost. Thus we can simply switch out the code for calculating the cost in the two trials to modify the problem. Monte Carlo estimation of $V(s)$ will work equally well in both cases -- in fact since P_{mod} depends on variables that we condition on the new cost terms doesn't add any additional MC computational cost. Now we are finally ready to implement Policy Iteration on this problem.Let $\theta$ be the policy evaluation error and $N$ the number of samples drawn for estimating each value of $V(s)$ by MC approximation.<jupyter_code># Set parameters
t = [20, 20]
b = 5
lambda_req = [3, 4]
lambda_ret = [3, 2]
m, p, cost_park = 10, 2, [4, 4]
theta, N = 1, 10**6
gamma = 0.9
# Define methods
def action_set(c_1, c_2, t=t, b=b):
    lower_bound = - min(t[0] - c_1, min(c_2, b))
    upper_bound = min(t[1] - c_2, min(c_1, b))
    return np.arange(lower_bound, upper_bound + 1)

def initialize(t=t, b=b):
    V = np.zeros((t[0] + 1, t[1] + 1))
    pi = np.zeros((t[0] + 1, t[1] + 1), dtype=int)

    return V, pi

def H(x):
    return int(x >= 0)

def state_policy_evaluation(c_1, c_2, a, V, gamma=gamma, P_mod=False, t=t, b=b, 
                            labmda_req=lambda_req, lambda_ret=lambda_ret, 
                            m=m , p=p, cost_park=cost_park, N=N):
    
    req_1 = np.random.poisson(lambda_req[0], N)
    req_2 = np.random.poisson(lambda_req[1], N)
    ret_1 = np.random.poisson(lambda_ret[0], N)
    ret_2 = np.random.poisson(lambda_ret[1], N)
    
    maxrent_1 = c_1 - a
    maxrent_2 = c_2 + a
    
    r_1 = m * np.minimum(req_1, maxrent_1)
    r_2 = m * np.minimum(req_2, maxrent_2)
    price = p * abs(a)
    
    if P_mod:
        price += - p * H(a - 1) + cost_park[0] * H(maxrent_1 - 10) \
                 + cost_park[1] * H(maxrent_2 - 10)
        
    r = r_1 + r_2 - price
    
    c_prime_1 = np.minimum(np.maximum(maxrent_1 - req_1, 0) + ret_1, t[0])
    c_prime_2 = np.minimum(np.maximum(maxrent_2 - req_2, 0) + ret_2, t[1])
    
    V_s = np.average(r + gamma * V[(c_prime_1, c_prime_2)])
    
    return V_s

def state_policy_evaluation_sum(c_1, c_2, a, V, gamma=gamma, P_mod=False, t=t, b=b, 
                                labmda_req=lambda_req, lambda_ret=lambda_ret, 
                                m=m , p=p, cost_park=cost_park, N=N):
    
    V_s = 0
    
    poisson_vals = np.arange(N + 1)
    
    req_probs_1 = stats.poisson(lambda_req[0]).pmf(poisson_vals)
    req_probs_2 = stats.poisson(lambda_req[1]).pmf(poisson_vals)
    ret_probs_1 = stats.poisson(lambda_ret[0]).pmf(poisson_vals)
    ret_probs_2 = stats.poisson(lambda_ret[1]).pmf(poisson_vals)
    
    maxrent_1 = c_1 - a
    maxrent_2 = c_2 + a
    
    price = p * abs(a)
    
    if P_mod:
        price += - p * H(a - 1) + cost_park[0] * H(maxrent_1 - 10) \
                 + cost_park[1] * H(maxrent_2 - 10)
    
    # Considered doing this as a tensor operation but it gets tricky.
    # May come back to that approach if this is too inefficient.
    for req_1, req_prob_1 in enumerate(req_probs_1):
        for req_2, req_prob_2 in enumerate(req_probs_2):
            for ret_1, ret_prob_1 in enumerate(ret_probs_1):
                for ret_2, ret_prob_2 in enumerate(ret_probs_2):
                    r_1 = m * min(req_1, maxrent_1)
                    r_2 = m * min(req_2, maxrent_2)
                    r = r_1 + r_2 - price
                    c_prime_1 = min(max(maxrent_1 - req_1, 0) + ret_1, t[0])
                    c_prime_2 = min(max(maxrent_2 - req_2, 0) + ret_2, t[1])
                    
                    V_s += req_prob_1*req_prob_2*ret_prob_1*ret_prob_2* \
                            (r + gamma * V[(c_prime_1, c_prime_2)])
    
    return V_s

def state_policy_evaluation_vec(c_1, c_2, a, V, gamma=gamma, P_mod=False, t=t, b=b, 
                                labmda_req=lambda_req, lambda_ret=lambda_ret, 
                                m=m , p=p, cost_park=cost_park, N=N):
    
    poisson_vals = np.arange(N + 1, dtype=int)
    
    # Probability tensor, shape (N + 1, N + 1, N + 1, N + 1)
    req_probs_1 = stats.poisson(lambda_req[0]).pmf(poisson_vals) # req_1 axis
    req_probs_2 = stats.poisson(lambda_req[1]).pmf(poisson_vals) # req_2 axis
    ret_probs_1 = stats.poisson(lambda_ret[0]).pmf(poisson_vals) # ret_1 axis
    ret_probs_2 = stats.poisson(lambda_ret[1]).pmf(poisson_vals) # ret_2 axis
    # axes ret_1, req_1, ret_2, req_2
    p_tensor = np.multiply.outer(
                   np.multiply.outer(
                        np.multiply.outer(ret_probs_1, req_probs_1), 
                           ret_probs_2), 
                               req_probs_2)
    
    assert p_tensor.shape == (N + 1, N + 1, N + 1, N + 1), "prob tensor is not the correct shape. " + \
                                                           "Current shape is {}".format(p_tensor.shape)
    
    # Price and bound scalars
    maxrent_1 = c_1 - a
    maxrent_2 = c_2 + a
    
    price = p * abs(a)
    
    if P_mod:
        price += - p * H(a - 1) + cost_park[0] * H(maxrent_1 - 10) \
                 + cost_park[1] * H(maxrent_2 - 10)
        
    # Reward matrix, shape (N + 1, N + 1, N + 1, N + 1)
    r_1 = m * np.minimum(poisson_vals, maxrent_1)
    r_2 = m * np.minimum(poisson_vals, maxrent_2)
    # axes ret_1, req_1, ret_2, req_2
    r = np.multiply.outer(
            np.multiply.outer(
                np.multiply.outer(np.ones(N + 1), r_1
                ), np.ones(N + 1) 
        ), r_2
    ) -  price
    
    assert r.shape == (N + 1, N + 1, N + 1, N + 1), "r tensor is not the correct shape. " + \
                                                    "Current shape is {}".format(r.shape)
        
    # c prime matrices, shape (N + 1, N + 1, N + 1, N + 1) each
    # axes ret_1, req_1
    c_prime_1 = np.minimum(np.add.outer(np.maximum(maxrent_1 - poisson_vals, 0), poisson_vals), t[0])
    # axes ret_1, req_1, ret_2, req_2
    c_prime_full_1 = np.multiply.outer(c_prime_1, np.ones((N + 1, N + 1), dtype=int))
    assert c_prime_full_1.shape == (N + 1, N + 1, N + 1, N + 1), "c prime 1 tensor is not the correct shape. " + \
                                                                 "Current shape is {}".format(c_prime_full_1.shape)
    # axes ret_2, req_2
    c_prime_2 = np.minimum(np.add.outer(np.maximum(maxrent_2 - poisson_vals, 0), poisson_vals), t[1])
    # axes ret_1, req_1, ret_2, req_2
    c_prime_full_2 = np.multiply.outer(np.ones((N + 1, N + 1), dtype=int), c_prime_2)
    assert c_prime_full_2.shape == (N + 1, N + 1, N + 1, N + 1), "c prime 2 tensor is not the correct shape. "  + \
                                                                 "Current shape is {}".format(c_prime_full_2.shape)
    
    # All the tensors being combined should agree on their dimensions
    V_s = np.sum(p_tensor * (r + gamma * V[c_prime_full_1, c_prime_full_2]))
    
    return V_s
    
def policy_evaluation(pi, V, gamma=gamma, P_mod=False, mode="sum", t=t, b=b, 
                      labmda_req=lambda_req, lambda_ret=lambda_ret, 
                      m=m , p=p, cost_park=cost_park, N=N, theta=theta):
    
    if mode == "sum":
        v_estimator = state_policy_evaluation_sum
    elif mode == "sum_vec":
        v_estimator = state_policy_evaluation_vec
    else:
        v_estimator = state_policy_evaluation
    
    while True:
        Delta = 0
        
        times = []
        
        for c_1 in range(t[0] + 1):
            for c_2 in range(t[1] + 1):
                t_i = time.time()
                v = V[c_1, c_2]
                V[c_1, c_2] = v_estimator(c_1, c_2, pi[c_1, c_2], V, gamma, P_mod, t, b, 
                                          labmda_req, lambda_ret, m , p, cost_park, N)
                Delta = max(Delta, abs(v - V[c_1, c_2]))
                times.append(time.time() - t_i)
                
        avg_t = round(sum(times)/len(times), 2)
        total_t = round(sum(times), 2)
                
        print("Current Delta value: {}. Avg. State Est. Time: {}s. Avg. Sweep Time: {}s".format(Delta,
                                                                                                avg_t,
                                                                                                total_t))
        
        if Delta < theta:
            break
                
                
def policy_improvement(pi, V, gamma=gamma, P_mod=False, mode="sum", t=t, b=b, 
                       labmda_req=lambda_req, lambda_ret=lambda_ret, 
                       m=m , p=p, cost_park=cost_park, N=N):
    
    if mode == "sum":
        v_estimator = state_policy_evaluation_sum
    elif mode == "sum_vec":
        v_estimator = state_policy_evaluation_vec
    else:
        v_estimator = state_policy_evaluation
    
    policy_stable = True
    
    for c_1 in range(t[0] + 1):
        for c_2 in range(t[1] + 1):
            
                old_action = pi[c_1, c_2]
                new_action = None
                new_value = - np.inf
                
                for a in action_set(c_1, c_2, t=t, b=b):
                    temp_value = v_estimator(c_1, c_2, a, V, gamma, P_mod, t, b, 
                                             labmda_req, lambda_ret, m, p, cost_park, N)
            
                    if temp_value > new_value:
                        new_value = temp_value
                        new_action = a
                
                pi[c_1, c_2] = new_action
                if new_action != old_action:
                    policy_stable = False
    
    print(pi)
                    
    return policy_stable

def policy_iteration(gamma=gamma, P_mod=False, mode="sum_vec", t=t, b=b, 
                     labmda_req=lambda_req, lambda_ret=lambda_ret, 
                     m=m , p=p, cost_park=cost_park, N=N, theta=theta):
    
    if mode not in ["sum", "mc", "sum_vec"]:
        raise ValueError("Parameter `mode` must be one of [`sum`, `mc`, `sum_vec`]")
    
    policy_stable = None
    V, pi = initialize(t, b)
    V_list = []
    pi_list = []
    
    print("Starting policy iteration with parameters \n" + \
          "\tN: {N}\n".format(N=N) + \
          "\tgamma: {gamma}\n".format(gamma=gamma) + \
          "\tP_mod: {P_mod}\n".format(P_mod=P_mod) + \
          "\tt: {t}\n".format(t=t) + \
          "\tb: {b}\n".format(b=b) + \
          "\tlambda_req: {lambda_req}\n".format(lambda_req=lambda_req) + \
          "\tlambda_ret: {lambda_ret}\n".format(lambda_ret=lambda_ret) + \
          "\tm: {m}\n".format(m=m) + \
          "\tp: {p}\n".format(p=p) + \
          "\tcost_park: {cost_park}\n".format(cost_park=cost_park) + \
          "\tt: {t}\n".format(t=t) + \
          "\ttheta: {theta}\n\n".format(theta=theta)
          )
    
    iteration = 0
    avg_t = 0
    
    while True:
        print("Iteration number {} start".format(iteration + 1))
        t_i = time.time()
        
        policy_evaluation(pi, V, gamma, P_mod, mode, t, b, labmda_req, 
                          lambda_ret, m, p, cost_park, N, theta)
        print("Policy Evaluation Complete")
        
        print(pi)
        policy_stable = policy_improvement(pi, V, gamma, P_mod, mode, t, b, labmda_req, 
                                           lambda_ret, m, p, cost_park, N)
        print("Policy Improvement Complete - Policy Stable: {}".format(policy_stable))
        
        V_list.append(V)
        pi_list.append(pi)
        
        if policy_stable:
            break
            
        delta_t = round(time.time() - t_i, 2)
        iteration += 1

        print("Iteration number {} complete: Latest iteration time {}s\n\n".format(iteration, delta_t))
            
    return V_list, pi_list<jupyter_output><empty_output><jupyter_text>Replication of Results: Summary of Attempts Fourth attempt (?): Cleaned vectorized sum trunctation approximation+ We used a trick to generate correct rank-4 tensors. Namely multiplying by all-1s tensors+ Could we perform better by directly broadcasting instead? Should be at least more memory efficient. It should also save on some opeerations. Could get a slight performance improvement+ Won't implement because current performance is satisfactory but wanted to summarize future attempt here Third attempt: Vectorized sum trunctation approximation+ Can we go even faster by vectorizing sum approach? YES!!!+ 3 seconds, N=10, sum approx. with vectorization for single sweep. ~0.01s for one state on average+ A greater than 10x improvement on performance for the same estimated O(0.01) error+ We finally achieved convergence!!! We believe the second approach would too but it just takes too long Second attempt: Unvectorized sum trunctation approximation+ 45 seconds, N = 10, sum approx. without vectorizing for single sweep. ~0.1s for one state. + 3x improvement for what we believe is much better error bounds O(0.01) First attempt: Vectorized MC approximation+ 127 seconds, N = 10^6, MC approximation fully vectorized. ~0.3s for one state+ We believe is O(1) error because of variance, based on a rough estimate and the fact we don't converge<jupyter_code>V_list, pi_list = policy_iteration(N=10, mode='sum_vec', theta=0.1)
fig, ax = plt.subplots(figsize=(10,10))         # Sample figsize in inches
sns.heatmap(pi_list[len(pi_list)-1].T, annot=True,  ax=ax)
ax.invert_yaxis()
plt.title("Policy Function - $\pi(C_1, C_2)$", size=20)
plt.xlabel("$C_1$", size=20)
plt.ylabel("$C_2$", size=20)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>New Environment<jupyter_code>V_list, pi_list = policy_iteration(N=10, P_mod=True, mode='sum_vec', theta=0.1)<jupyter_output><empty_output><jupyter_text>We can see in the plot below that: + Even at high $C_1$ and $C_2$ we still end up moving the one free car since it incurs no cost to us.+ We see lines at $C_1=10$ and $C_2=10$ where we'd rather move cars away from a lot that has fewer just to avoid the parking lot fee.+ Besides this we see essentialy the same features we saw in the original case, which makes sense!<jupyter_code>fig, ax = plt.subplots(figsize=(10,10))         # Sample figsize in inches
sns.heatmap(pi_list[len(pi_list)-1].T, annot=True,  ax=ax)
ax.invert_yaxis()
plt.title("Policy Function - $\pi(C_1, C_2)$", size=20)
plt.xlabel("$C_1$", size=20)
plt.ylabel("$C_2$", size=20)
plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10988142292490119, 'fraction_numerical': 0.023715415019762844, 'mean_word_length': 2.6129470964232286, 'pattern_counts': {'":': 4, '<': 26, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '316157', 'license': 'MIT', 'n_tokens_mistral': 9235, 'n_tokens_neox': 8732, 'n_words': 2732, 'path': 'S&B RL Ex. 2.5, 4.7.ipynb', 'repo_name': 'dbulhosa/suttonandbarto', 'size': '26565'}