hw3_AC_copy.ipynb

{"cells":[{"cell_type":"markdown","metadata":{"id":"0Nbq_oB2NI3E"},"source":["# 546 (Optimization for Learning and Control) hw3\n","\n","You are welcome (and encouraged) to work with others, but each individual must submit their own writeup.\n","\n","You are welcome to use your preferred computational toolkit -- upload the **commented** sourcecode alongside your writeup (e.g. the .ipynb file).\n","\n","You are welcome to consult research articles and other materials -- upload a .pdf of the article alongside your writeup and indicate which references were used where in the writeup."]},{"cell_type":"markdown","metadata":{"id":"1yz1n86PNI3H"},"source":["# Markov decision process (MDP) (14 pts)"]},{"cell_type":"markdown","metadata":{"id":"BoRzaEsNNI3I"},"source":["Let $(X,U,P,c)$ be an MDP:\n","\n","* $X$ is a finite set of $N$ *states*;\n","* $U$ is a finite set of $M$ *actions*;\n","* $P:X\\times U\\rightarrow \\Delta(X)$ specifies transition probabilities;\n","* $C:X\\times U\\rightarrow \\mathbb{R}$ specifies transition cost.\n","\n","Recall that $\\Delta(S) = \\{p : S \\rightarrow [0,1] : \\sum_{s\\in S} p(s) = 1\\}$ is the set of probability distributions over the finite set $S$."]},{"cell_type":"code","execution_count":314,"metadata":{},"outputs":[],"source":["import numpy as np\n","import pylab as plt\n","from scipy import linalg as la\n","import scipy.stats as stats\n","%matplotlib inline\n","import matplotlib.pylab as plt\n","import control as ctrl"]},{"cell_type":"markdown","metadata":{"id":"THQ7pkHANI3I"},"source":["a. *Implement an algorithm that generates a random MDP given finite sets $X$ and $U$ (i.e. let $N = |X|$, $M = |U|$ be parameters that are easy to vary, and generate $P$ and $c$ randomly).*"]},{"cell_type":"code","execution_count":315,"metadata":{},"outputs":[{"data":{"text/plain":["(array([[[0.10191865, 0.3195559 , 0.57852545],\n","         [0.61858119, 0.0608104 , 0.32060841]],\n"," \n","        [[0.07266062, 0.50025493, 0.42708445],\n","         [0.05598505, 0.82785363, 0.11616132]],\n"," \n","        [[0.40308553, 0.22569177, 0.3712227 ],\n","         [0.55298121, 0.00147458, 0.44554422]]]),\n"," array([[[0.59752625, 0.34340974, 0.60161398],\n","         [0.20718465, 0.28319941, 0.21725186]],\n"," \n","        [[0.42818152, 0.70044246, 1.90565948],\n","         [0.38738289, 0.29332724, 0.42693463]],\n"," \n","        [[0.43790514, 1.25689285, 0.47121546],\n","         [1.19326255, 0.41595356, 0.65965854]]]))"]},"execution_count":315,"metadata":{},"output_type":"execute_result"}],"source":["#generate random MDP (X,U,P,C)\n","\n","#X is a set of N sample space (finite set)\n","N = 3\n","X = list(range(N))\n","\n","#U is a set of M sample space (finite set)\n","M = 2\n","U = list(range(M))\n","\n","#P is transition probabilities\n","P_ = np.abs(np.random.randn(N,M,N))\n","P = P_ / np.sum(P_, axis=2)[:,:,np.newaxis] #sum to be 1 over X\n","P\n","\n","#C is a random cost \n","C = np.abs(np.random.randn(N,M,N))\n","\n","P,C"]},{"cell_type":"markdown","metadata":{"id":"PUvt8PuMNI3J"},"source":["b. *Implement an algorithm that simulates an MDP (i.e. write a function that inputs transition probabilities $P$, cost $c$, initial state distribution $p_0$, policy $\\pi:X\\rightarrow\\Delta(U)$, and time horizon $t$ and returns the state distribution $p_t$).*"]},{"cell_type":"code","execution_count":316,"metadata":{},"outputs":[],"source":["#p0 is a random initila state distribution\n","p0_ = np.abs(np.random.randn(N)) #1x3\n","p0 = p0_ / np.sum(p0_)\n","\n","#pi is policy with sum=1 over U \n","PI_ = np.abs(np.random.randn(N,M)) #3x2\n","PI = PI_ / np.sum(PI_,axis=1)[:,np.newaxis]\n"]},{"cell_type":"markdown","metadata":{},"source":["$$p(x^{+}) = \\sum_{x\\in X} p(x) (\\sum_{u \\in U} \\pi (u|x) P(x^+|x,u)) .$$"]},{"cell_type":"code","execution_count":317,"metadata":{},"outputs":[],"source":["def MP(X,U,P,p0,PI, T = 1): #markov process\n","    p_ = [p0]\n","    for t in range(T):\n","        p_prev = p_[-1]\n","        p_next = np.zeros((N))\n","        for xp in X:\n","            p_next[xp] = np.sum( [p_prev[x]* np.sum([PI[x,u]*P[x,u,xp] for u in U]) for x in X]) #x is x, xp is x+\n","        p_.append(p_next)\n","    return p_"]},{"cell_type":"markdown","metadata":{"id":"FZ48dOXuNI3J"},"source":["c. *Test your algorithm from (b.) (i.e. design one or more tests, explain why these tests are exhaustive, and provide summary statistics and/or visualizations that convince a skeptical reader that your algorithm is correct).*"]},{"cell_type":"code","execution_count":318,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["[0.32499672 0.24425737 0.43074591]\n","[0.32499672 0.24425737 0.43074591]\n","[0.32499672 0.24425737 0.43074591]\n"]}],"source":["T = 50\n","\n","for i in range(3): # 3 experiments with different initial cond p0\n","    p0_ = np.abs(np.random.randn(N)) #1x3\n","    p0 = p0_ / np.sum(p0_)\n","\n","    p_ = MP(X,U,P,p0,PI, T = T) #100 \n","    print(p_[-1]) #final p after T iterations\n","\n"]},{"cell_type":"markdown","metadata":{"id":"sDYop_ybNI3K"},"source":["Let us now consider the problem of minimizing the infinite-horizon discounted cost\n","$$c = \\sum_{t=0}^\\infty \\gamma^t C(x_t,u_t),$$\n","where $\\gamma\\in(0,1)$ is a ***discount factor***."]},{"cell_type":"markdown","metadata":{"id":"91q5DHAwNI3K"},"source":["Any policy $\\pi : X\\rightarrow\\Delta(U)$ has an associated ***value*** $v^\\pi : X\\rightarrow\\mathbb{R}$ defined by\n","$$\\forall x\\in X : v^\\pi(x) = E[c \\mid x = x_0]$$\n","that satisfies the *Bellman equation*\n","$$\\forall x\\in X : v^\\pi(x) = \\sum_{u\\in U}\\pi(u|x)\\sum_{x^+\\in X} P(x^+|x,u)\\left[C(x,u) + \\gamma \\cdot v^\\pi(x^+)\\right].$$"]},{"cell_type":"markdown","metadata":{"id":"_RoYRpTtNI3K"},"source":["d. *Noting that $v^\\pi$ appears linearly in this Bellman equation, implement a ***policy evaluation*** algorithm that computes $v^\\pi$ using linear algebra (i.e. determing $L$ and $b$ such that $L v^\\pi = b$ and use this equation to solve for $v^\\pi$).*"]},{"cell_type":"markdown","metadata":{},"source":["$$ v^{\\pi}(x^+) = L^{-1} b$$"]},{"cell_type":"markdown","metadata":{},"source":["where $$L(x,x^+) = \\delta(x,x^+) - \\sum_{u\\in U} \\gamma\\cdot \\pi(u|x)\\cdot P(x^+|x,u),$$\n","$$b(x) = \\sum_{x^+\\in X}\\sum_{u\\in U} \\pi(u|x)\\cdot P(x^+|x,u)\\cdot C(x,u)$$\n","\n","Note that the \"C\" here is the \"L\" in lecture note, where $L: X \\times U \\rightarrow \\Delta(X)$"]},{"cell_type":"code","execution_count":319,"metadata":{},"outputs":[],"source":["def policy_evaluation(gamma,X,U,P,C,PI): #policy evaluation: given the current pi, find the next v\n","    b = np.zeros((N))\n","    I = np.identity(N)\n","    L = I\n","\n","    for x in X:\n","        b[x] = np.sum( [[PI[x,u]*P[x,u,xp]*C[x,u,xp] for u in U] for xp in X] )  #x is x, xp is x+\n","        L = I - [np.sum([gamma*PI[x,u]*P[x,u,xp] for u in U]) for xp in X]\n","    v_pi = np.dot(la.inv(L),b)\n","\n","    return v_pi\n","    "]},{"cell_type":"markdown","metadata":{"id":"K0_oTUisNI3L"},"source":["e. *Test your policy evaluation algorithm from (d.) (i.e. design one or more tests, explain why these tests are exhaustive, and provide summary statistics and/or visualizations that convince a skeptical reader that your algorithm is correct).*"]},{"cell_type":"code","execution_count":320,"metadata":{},"outputs":[{"data":{"text/plain":["array([3.78836181, 4.30140692, 4.25000789])"]},"execution_count":320,"metadata":{},"output_type":"execute_result"}],"source":["gamma = np.random.rand()\n","v = policy_evaluation(gamma,X,U,P,C,PI) #value function: value of the policy for each state\n","v"]},{"cell_type":"markdown","metadata":{"id":"Tv56gQ-oNI3L"},"source":["f. *Using your policy evaluation algorithm from (d.), implement ***value iteration*** and ***policy iteration*** algorithms.*\n","\n","g. *Test your value iteration and policy iteration algorithms from (f.) (i.e. design one or more tests, explain why these tests are exhaustive, and provide summary statistics and/or visualizations that convince a skeptical reader that your algorithm is correct).*"]},{"cell_type":"markdown","metadata":{},"source":["$$u^*(x)= \\text{argmin}_{ u\\in U}\\left(\\sum_{x^+\\in X} P(x^+|x,u)\\left[C(x,u) + \\gamma \\cdot v^\\pi(x^+)\\right]\\right)$$\n","And then we update policy by $$\\pi(u|x)=\\delta(u,u^*(x))$$"]},{"cell_type":"code","execution_count":321,"metadata":{},"outputs":[],"source":["#policy iteration: \n","\n","#(1) policy evaluation: given any initial policy pi0, find next value of the policy v1\n","\n","#(2) greedy_update: improve policy with greedy update pi+ (i.e. make the direction of next policy \"1\")\n","\n","def greedy_update(gamma,X,U,P,C,v_pi):\n","    pi = np.zeros((N,M))\n","    for x in X:\n","        index = np.argmin( [np.sum([P[x,u,xp]*(C[x,u,xp] + gamma*v_pi[xp]) for xp in X]) for u in U])\n","        pi[x,index] = 1.\n","    return pi\n","\n","def policy_iteration(gamma,X,U,P,C,PI0): #finite number of iterations\n","    PI_ = [PI0]\n","    v1 = policy_evaluation(gamma,X,U,P,C,PI0)\n","    v_ = [v1]\n","\n","    while len(v_)<2 or not(np.allclose(v_[-1],v_[-2])): # loop until v_pi converges\n","        PI_.append( greedy_update(gamma,X,U,P,C,v_[-1]) )\n","        v_.append( policy_evaluation(gamma,X,U,P,C,PI_[-1]) )\n","        \n","    return PI_, v_\n","\n","# def policy_iteration(gamma,X,U,P,C,PI0): #finite number of iterations\n","#     PI_ = [PI0]\n","#     v_ = [0.]\n","\n","#     while len(v_)<2 or not(np.allclose(v_[-1],v_[-2])): # loop until v_pi converges\n","#         v_.append( policy_evaluation(gamma,X,U,P,C,PI_[-1]) )\n","#         PI_.append( greedy_update(gamma,X,U,P,C,v_[-1]) )\n","        \n","#     return PI_, v_\n"]},{"cell_type":"code","execution_count":322,"metadata":{},"outputs":[{"data":{"text/plain":["([array([[0.70088204, 0.29911796],\n","         [0.71123107, 0.28876893],\n","         [0.203178  , 0.796822  ]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]])],\n"," [array([3.78836181, 4.30140692, 4.25000789]),\n","  array([2.24640079, 2.34547899, 2.6664755 ]),\n","  array([2.24640079, 2.34547899, 2.6664755 ])])"]},"execution_count":322,"metadata":{},"output_type":"execute_result"}],"source":["PI0 = PI\n","PI_, v_ = policy_iteration(gamma,X,U,P,C,PI0)\n","PI_,v_"]},{"cell_type":"markdown","metadata":{},"source":["$$v^\\pi(x) = \\sum_{u\\in U}\\pi(u|x)\\sum_{x^+\\in X} P(x^+|x,u)\\left[C(x,u) + \\gamma \\cdot v^\\pi(x^+)\\right]$$"]},{"cell_type":"code","execution_count":323,"metadata":{},"outputs":[],"source":["#value iteration: v+ = Tv, iterate until v converges \n","\n","#(1) update_value : given any initial policy pi0 and initial value of the policy v0, find the next value v1\n","\n","#(2) greedy_update: from v1, improve policy with greedy update pi1 (i.e. make the direction of next policy \"1\")\n","    \n","def update_value(gamma,X,U,P,C,PI,v_pi):\n","    v_pi_next = np.zeros((N))\n","    for x in X:\n","        v_pi_next[x] = np.sum([PI[x,u] * np.sum([P[x,u,xp]*(C[x,u,xp] + gamma*v_pi[xp]) for xp in X]) for u in U])\n","    return v_pi_next\n","    \n","def value_iteration(gamma,X,U,P,C,PI0,maxiter=100): #infinite number of iteraitons\n","    PI_ = [PI0]\n","    v0 = np.zeros((N))\n","    v_ = [v0]\n","\n","    while len(v_)<2 or not(np.allclose(v_[-1],v_[-2])) and len(v_) < maxiter: # loop until v_pi converges or reaches max iteration\n","        v_.append( update_value(gamma,X,U,P,C,PI_[-1],v_[-1]))\n","        PI_.append( greedy_update(gamma,X,U,P,C,v_[-1]))\n","        \n","    return PI_, v_"]},{"cell_type":"code","execution_count":324,"metadata":{},"outputs":[{"data":{"text/plain":["([array([[0.70088204, 0.29911796],\n","         [0.71123107, 0.28876893],\n","         [0.203178  , 0.796822  ]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]]),\n","  array([[0., 1.],\n","         [0., 1.],\n","         [1., 0.]])],\n"," [array([0., 0., 0.]),\n","  array([0.42785894, 0.94090405, 0.88950501]),\n","  array([0.72362125, 1.07331741, 1.23414049]),\n","  array([0.97621014, 1.21256509, 1.4662377 ]),\n","  array([1.17654558, 1.34357585, 1.65004793]),\n","  array([1.33641201, 1.4617238 , 1.79963751]),\n","  array([1.46545905, 1.56572225, 1.92248604]),\n","  array([1.57063096, 1.6558593 , 2.02393503]),\n","  array([1.65697586, 1.73318028, 2.10804525]),\n","  array([1.72825387, 1.79904215, 2.17798198]),\n","  array([1.78733291, 1.85487045, 2.23625599]),\n","  array([1.83644627, 1.90203272, 2.28488617]),\n","  array([1.87736323, 1.94177866, 2.32551314]),\n","  array([1.91150507, 1.97521765, 2.35948092]),\n","  array([1.94002583, 2.00331657, 2.38789722]),\n","  array([1.96387039, 2.02690793, 2.41167907]),\n","  array([1.98381712, 2.04670269, 2.43158817]),\n","  array([2.00051022, 2.06330459, 2.44825869]),\n","  array([2.01448465, 2.07722428, 2.46221955]),\n","  array([2.0261857 , 2.08889247, 2.47391247]),\n","  array([2.03598474, 2.09867179, 2.48370663]),\n","  array([2.04419186, 2.10686707, 2.49191081]),\n","  array([2.05106622, 2.11373433, 2.49878341]),\n","  array([2.05682458, 2.11948843, 2.50454072]),\n","  array([2.06164832, 2.1243096 , 2.50936382]),\n","  array([2.06568924, 2.12834899, 2.51340437]),\n","  array([2.06907447, 2.1317333 , 2.51678936]),\n","  array([2.07191043, 2.1345687 , 2.51962519]),\n","  array([2.07428627, 2.13694422, 2.52200095]),\n","  array([2.07627667, 2.13893441, 2.52399129]),\n","  array([2.07794415, 2.14060178, 2.52565875]),\n","  array([2.07934113, 2.14199868, 2.52705571]),\n","  array([2.08051147, 2.14316898, 2.52822604]),\n","  array([2.08149196, 2.14414944, 2.52920652]),\n","  array([2.08231339, 2.14497086, 2.53002795]),\n","  array([2.08300156, 2.14565902, 2.53071612]),\n","  array([2.0835781 , 2.14623555, 2.53129265]),\n","  array([2.08406111, 2.14671856, 2.53177566]),\n","  array([2.08446576, 2.14712321, 2.53218032]),\n","  array([2.08480477, 2.14746222, 2.53251933]),\n","  array([2.08508879, 2.14774624, 2.53280335]),\n","  array([2.08532673, 2.14798418, 2.53304129]),\n","  array([2.08552608, 2.14818352, 2.53324063]),\n","  array([2.08569308, 2.14835053, 2.53340764]),\n","  array([2.085833  , 2.14849044, 2.53354755]),\n","  array([2.08595021, 2.14860766, 2.53366477]),\n","  array([2.08604842, 2.14870586, 2.53376297]),\n","  array([2.08613069, 2.14878813, 2.53384524]),\n","  array([2.08619961, 2.14885706, 2.53391417]),\n","  array([2.08625736, 2.1489148 , 2.53397191]),\n","  array([2.08630574, 2.14896318, 2.53402029]),\n","  array([2.08634626, 2.14900371, 2.53406082]),\n","  array([2.08638022, 2.14903767, 2.53409477]),\n","  array([2.08640867, 2.14906611, 2.53412322]),\n","  array([2.0864325 , 2.14908994, 2.53414705]),\n","  array([2.08645246, 2.14910991, 2.53416702])])"]},"execution_count":324,"metadata":{},"output_type":"execute_result"}],"source":["PI0 = PI\n","PI_, v_ = value_iteration(gamma,X,U,P,C,PI0,maxiter=100)\n","PI_,v_"]},{"cell_type":"markdown","metadata":{"id":"wPdI3b--NI3M"},"source":["## linear quadratic regulation (6 pts)"]},{"cell_type":"markdown","metadata":{"id":"n0Sz3znwNI3M"},"source":["Consider the linear-quadratic regulation (LQR) problem for the discrete-time linear time-invariant system\n","$$x^+ = A x + B u$$\n","with infinite-horizon cost\n","$$c(x,u) = \\frac{1}{2}\\sum_{t=0}^\\infty x^T Q x + u^T R u \\ dt$$\n","where  $Q = Q^T > 0$ and $R = R^T > 0$.  We know that the optimal policy is linear in state, \n","$$u = - K x,$$\n","where $K = (R + B^T Y B)^{-1} (B^T Y A)$ and $Y$ satisfies the discrete algebraic Riccati equation\n","$$Y = A^T Y A - (A^T Y B)(R + B^T Y B)^{-1}(B^T Y A) + Q.$$"]},{"cell_type":"markdown","metadata":{"id":"LmelWMGTNI3N"},"source":["For this problem, we'll consider a one-dimensional system, so $x,u\\in\\mathbb{R}$, with the values $A, B, Q, R = 1$."]},{"cell_type":"markdown","metadata":{"id":"vBV7m_UrNI3N"},"source":["If the state distribution starts out Gaussian, \n","$$x \\sim \\mathcal{N}(\\xi, \\Sigma),$$\n","and we apply zero input but add a state disturbance with zero mean and covariance $W$, then (so long as the disturbance is uncorrelated with the state) the state distribution after one step is Gaussian:\n","$$x^+ \\sim \\mathcal{N}(A\\xi, A\\Sigma A^T + W),$$\n","i.e. the covariance gets updated via\n","$$\\Sigma^+ = A\\Sigma A^T + W.$$\n","Since $A$ is stable, so long as $W$ is constant, this iteration will converge to a solution of the Lyapunov equation\n","$$S = A S A^T + W.$$\n","Also, since $A$ is stable, the steady-state mean is zero, thus we can determine the steady-state distribution in closed-form:\n","$$\\lim_{t\\rightarrow\\infty} x(t) \\sim \\mathcal{N}(0, S).$$"]},{"cell_type":"markdown","metadata":{"id":"iyhG0qxJNI3N"},"source":["In this problem, you'll make use of the facts above to build a finite MDP that approximates the LQR, solve the MDP using your algorithms from the previous problem, then make use of the facts above again to validate the solution you've obtained."]},{"cell_type":"markdown","metadata":{"id":"Ew--O5WQNI3O"},"source":["a. *Construct a finite MDP by discretizing the state and action spaces of the linear system and computing transition probabilities and transition costs to obtain  $(X,U,P,C)$ where:*\n","\n","* $X$ is a finite set of $N$ *states*;\n","* $U$ is a finite set of $M$ *actions*;\n","* $P:X\\times U\\rightarrow \\Delta(X)$ specifies transition probabilities;\n","* $C:X\\times U\\times X\\rightarrow \\mathbb{R}$ specifies transition costs."]},{"cell_type":"code","execution_count":410,"metadata":{},"outputs":[],"source":["n = 1; m = 1\n","A = np.array([[0.99]]) # to avois system becomes unstable\n","B = Q  = R = np.array([[1.0]])"]},{"cell_type":"code","execution_count":411,"metadata":{},"outputs":[{"data":{"text/plain":["(array([[0.60976745]]), array([[1.60366977]]))"]},"execution_count":411,"metadata":{},"output_type":"execute_result"}],"source":["Y = la.solve_discrete_are(A,B,Q,R)\n","K, Y, E = ctrl.dlqr(A, B, Q, R)\n","K, Y"]},{"cell_type":"code","execution_count":412,"metadata":{},"outputs":[{"data":{"text/plain":["array(100.50251256)"]},"execution_count":412,"metadata":{},"output_type":"execute_result"}],"source":["W = 2.\n","S = la.solve_discrete_lyapunov(A,W) #np.dot(A, np.dot( S, A.T )) + W\n","S "]},{"cell_type":"code","execution_count":413,"metadata":{},"outputs":[],"source":["def MDP_LQR(N,M,n,m): #x , u \\in R, N and M are number of sample space\n","    X_index = list(range(N))\n","    U_index = list(range(M))\n","    X = 15*2*(np.random.rand(N,n)-.5)\n","    U = 15*2*(np.random.rand(M,m)-.5)\n","    X[0],U[0] = 0,0\n","    P = np.zeros((N,M,N))\n","    C = np.zeros((N,M,N))\n","\n","    for xi in X_index:\n","        for ui in U_index:\n","            x_next = np.dot(X[xi],A.T)+np.dot(U[ui],B.T) #deterministic #np.dot(x,A.T)+np.dot(u,B.T)\n","            P[xi,ui] = stats.multivariate_normal.pdf(X, mean=x_next, cov=W)\n","            C[xi,ui] = 0.5* ( np.dot(np.dot(X[xi],Q),X[xi]) + np.dot(np.dot(U[ui],R),U[ui]) )\n","\n","    P = P / P.sum(axis=2)[...,np.newaxis]  #normalization to make sum prob = 1\n","\n","    return X,U,P,C\n"]},{"cell_type":"code","execution_count":415,"metadata":{},"outputs":[],"source":["N,M = 1000,5                         # number of discrete states, actions\n","X,U,P,C = MDP_LQR(N,M,n,m)   "]},{"cell_type":"markdown","metadata":{"id":"vC1d1EnyNI3P"},"source":["b. *Compare the steady-state distribution of the MDP with the steady-state distribution of the linear system subject to zero-mean disturbance with covariance $W$.  Discuss how the comparison varies with respect to $|X|$.*"]},{"cell_type":"markdown","metadata":{},"source":["1.  MDP method: the unique steady-state distribution is the left eigenvector associated with the unity eigenvalue.\n","$$ p^+(x^+) = \\sum_{x \\in X} (p(x) (\\sum_{u \\in U} \\pi (u|x) P(x^+|x,u)))$$\n","$$ p^+ = \\Gamma p $$\n","(i.e. $\\Gamma$ is the left stochastic) --> $\\forall \\lambda \\in spec\\Gamma : |\\lambda| \\leq 1$), so the system is stable.\n","\n","$$[\\Gamma]_{x^+,x} = \\sum_{u \\in U} \\pi (u|x) P(x^+|x,u) $$"]},{"cell_type":"code","execution_count":416,"metadata":{},"outputs":[],"source":["pi = np.zeros((N,M)); pi[:,0] = 1. #make sum be 1 --> zero input\n","G = np.sum(pi[...,np.newaxis]*P,axis=1) #Gamma, here is same as P\n","eigval,eigvec_left_T,eigvec_right = la.eig(G,left=True) #rows of eigvec_left_T is the left eigvectors of G\n","maxeigvec = eigvec_left_T.T[np.argmax(np.abs(eigval))] #eigvec corresponding to the largest eigvals\n","maxeigvec = maxeigvec/maxeigvec.sum() #normalization, we are finding the eigvec associated with eigval = 1\n","# maxeigvec.real"]},{"cell_type":"markdown","metadata":{},"source":["2. Lyapunov method:\n","Since $A$ is stable, so long as $W$ is constant, this iteration will converge to a solution of the Lyapunov equation\n","$$S = A S A^T + W.$$\n","Also, since $A$ is stable, the steady-state mean is zero, thus we can determine the steady-state distribution in closed-form:\n","$$\\lim_{t\\rightarrow\\infty} x(t) \\sim \\mathcal{N}(0, S).$$"]},{"cell_type":"code","execution_count":417,"metadata":{},"outputs":[],"source":["Lyapunov = stats.multivariate_normal.pdf(X,mean=0.,cov=S) #steady state distribution\n","Lyapunov = Lyapunov/Lyapunov.sum()\n","# Lyapunov"]},{"cell_type":"code","execution_count":418,"metadata":{},"outputs":[{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYcAAAEICAYAAAC0+DhzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU5dnw8d+VjV0IEEEMJESiAg0CCYtCUVAQ0AqIsqkIKIgVlWqriLSgj+WhtipaeVSsFvEVkYooFVyiFBEUQgIUxEjBSCACliUii5BMcr9/nMkwy5lkksxkslzfzyck5z7bfRJmrrl3McaglFJKuYsIdwaUUkpVPxoclFJK+dDgoJRSyocGB6WUUj40OCillPKhwUEppZQPDQ5KKaV8aHBQNYaI7BWRH0SkkVvanSKyNsj3iRGRt533MyJyldd+EZE/ichR59eTIiJu+xNF5F8iclpEvhGRa8px77XOe17mlf6uXV4qS0QWiUiBiJx0+4oM5j1UzaTBQdU0UcD9VXCf9cCtwCGbfVOA4cBlQBfgeuAut/1vAluBFsCjwNsiEleOe/8HGF+yISItgN7A4XJcozyeNMY0dvsqCtF9VA2iwUHVNH8GfisizUJ1A2NMgTFmvjFmPWD3Rnk78JQxJs8Y8z3wFDABQEQuBroDs40xPxtjlgM7gJHlyMIbwGi3T/BjgRVAQYUeSKkK0OCgappMYC3w20AOFpEfS/maUcE8dAb+7bb9b2dayb4cY8wJP/sDcQD4Ghjk3B4PLC7tBBGZUdqzlnG/X4vIMRHJEpHyBDFVi0WFOwNKVcAfgA0i8mxZBxpjQlHCaAwcd9s+DjR2tjt47yvZf2E577EYGC8iOUAzY8yXbs0aPowx84B55bwHwHPAg848DgLeEpFDxpgNFbiWqkW05KBqHGPMV8D7QEU/+VfWSeA8t+3zgJPGmsXSe1/J/hOUzzvAAOBe4PUK5rNMxpgtxpijxhiHMWY1VpXWjaG6n6o5NDiommo2MJkyPpF79cLx/ppZwXvvxGqMLnGZM61kX5KINPGzPyDGmNPAB8DdBBAcRGRmac9anlsD/osoqs7Q4KBqJGPMHuAt4L4yjmtcytdcf+eJSD0Rqe/cjBGR+m7dVRcDD4jIhSLSBqtaZpHzfv8BtgGzneeMwOrRtNx53URnl9TEAB5zJnClMWZvWQcaY+aW9qylPOdNItJYRCJEZBBWD62VAeRN1XLa5qBqsseB20J07V1AgvPnj5zf2wN7gZeAJKxeSAB/c6aVGIMVLPKBfcBNxpiSbqhtgVzg+7IyYIw5gNU4HUr3A69glRa+AyYbY9aG+J6qBhBd7EepqiMis4DDxpiXyjxYqTDS4KCUUsqHtjkopZTyocFBKaWUDw0OSimlfNSK3kotW7Y0iYmJ4c6GUkrVKFlZWUeMMbaTQtaK4JCYmEhmZma4s6GUUjWKiOT626fVSkoppXxocFBKKeVDg4NSSikftaLNwU5hYSF5eXmcOXMm3FlRVaB+/frEx8cTHR0d7qwoVSvU2uCQl5dHkyZNSExMpLR58FXNZ4zh6NGj5OXl0b59+3BnR6laodZWK505c4YWLVpoYKgDRIQWLVpoKVGpIKq1JQdAA0MdUpP+1uNf2cT6PUfAQL2oCBrUi2RUaltmDO0Y7qwp5VJrSw5KVTdLNu3jkkdXs273EYoNFAM/O4o5dqqQF9flkDhjFdOXbg13NpUCNDiElIhw223nlhtwOBzExcVx/fXXA7Bo0SLi4uLo2rUrnTp14uWXX/ZI79atG8nJyVx77bV88cUXIcnj/PnzOX36tGu7cWO/68IEpLLn11bjX9nEzBU7OFtU+izI7247QMdZH1RRrpTyT4NDCDVq1IivvvqKn3/+GYD09HQuvNBzVcvRo0ezbds21q5dy8yZM/nhhx9c6Vu3bmX37t3MmDGDG2+8kezs7KDn0Ts4qODKys0nZfaHrNt9JOBzfnYU02HmKrJy80OYM6VKp8HBTVZuPgv+tSeoL8ohQ4awatUqAN58803Gjh1re9z555/PRRddRG6u72j2/v37M2XKFBYuXOizb8KECdx9993079+fpKQkPvvsMyZNmkTHjh2ZMGGC67iPP/6Yyy+/nO7du3PzzTdz8uRJnnvuOQ4cOED//v3p37+/69hHH32Uyy67jN69e7uCVW5uLldffTVdunTh6quvZt++fQB89913XH755fTo0YPf//73Ff491UZZufmMfOELTpwtKve5jmK46YUvNECosNHg4JSVm88tf9vIUx/v4pa/bQzai3LMmDEsXbqUM2fOsH37dnr16mV7XE5ODjk5OXTo0MF2f/fu3fnmm29s9+Xn57NmzRqeeeYZfvWrX/Gb3/yGnTt3smPHDrZt28aRI0d44okn+OSTT9iyZQtpaWk8/fTT3HfffbRp04Z//etf/Otf/wLg1KlT9O7dm3//+9/069fPVdU1bdo0xo8fz/bt27nlllu47z5r6eb777+fu+++m82bN9O6devK/rpqldv+ttE2PQKIiYxgeNc2zB2RQsNo+5ehAaa9kRW6DCpVilrdW6k8NuYcpcBRTLGBQkcxG3OOkpoQW+nrdunShb179/Lmm28ydOhQn/1vvfUW69evp169erz00ks0b97c9jqlrdj3q1/9ChEhJSWFVq1akZKSAkDnzp3Zu3cveXl5fP311/Tp0weAgoICLr/8cttrxcTEuNpEUlNTSU9PB+DLL7/knXfeAeC2227joYceAmDDhg0sX77clf7www+X+TupCwY+tZbThcU+6fHN6rN+xtUeaeN6tSMrN59RL36Bd5PEwZ/OMvCptaQ/eFUIc6uULw0OTr2TWhATFUGho5joqAh6J7UI2rVvuOEGfvvb37J27VqOHj3qsW/06NE8//zzZV5j69atdOxo39WxXr16AERERLh+Ltl2OBxERkYycOBA3nzzzTLvEx0d7eoWGhkZicPhsD3OvetoTepGWhXGv7KJ3YdP+aQ3jon0CQwlUhNi+fZ/r2P8K5t82id2Hz7F9KVbmT+mW0jyq5QdrVZySk2I5Y07e/PAoEt4487eQSk1lJg0aRJ/+MMfXJ/oy+uzzz5j4cKFTJ48uULn9+7dmw0bNrBnzx4ATp8+zX/+8x8AmjRpwokTJ8q8xhVXXMHSpUsBeOONN+jbty8Affr08Uiv6+atzrZtfI5rHMNXjw8u8/zFd/Sia3xTn/R3tx3Qbq6qSmlwcJOaEMs9/TsENTAAxMfHc//995frnLfeeouuXbty8cUXM3fuXJYvX+635FCWuLg4Fi1axNixY+nSpQu9e/d2tV9MmTKFIUOGeDRI23nuuef4+9//TpcuXXj99dd59tlnAXj22WdZsGABPXr04Pjx4xXKX22xZNM+XlyX45PeIDqCzbMGBnydd6f1JTmukW+6BghVhaS0uuyaIi0tzXgv9pOdnV3hN1NVM4X7b37bK5v43KbUMHdECuN6tSv39YbMX0f2Id9SXUWvp5Q3EckyxqTZ7dOSg1JBYtfyMrxrmwq/kT8xIsX2BfpM+q4KXU+p8tDgoFQQeDckN2sQxdwRKZVqRE5NiOUfd19Bx9ZNPNIPnyxg+PPrK3xdpQKhwUGpSrLrYdSsYUxQqn5SE2L5YHo/GniNhdiWd5x5q4M/Yl6pEhoclKoEu8AAMLhzcAcEXmtzvRfX5egIahUyAQUHERksIrtEZI+IzLDZLyLynHP/dhHpXta5InKziOwUkWIR8WkQEZF2InJSRH5b0YdTKpSWbNpnGxgSmjcM+vTb88d0o/V59XzS532gpQcVGmUGBxGJBBYAQ4BOwFgR6eR12BAg2fk1BXghgHO/Am4E1vm59TOATk+pqq23Nu+zTX96dNeQ3G/BLak+jd6b91ojq7UEoYItkJJDT2CPMSbHGFMALAWGeR0zDFhsLBuBZiJyQWnnGmOyjTG23S5EZDiQA+ys0FNVEzp9tWXRokUcOHDAtZ2YmMiRI4HPUuqtsucHS6vz6vukTe2XFPRxMiVSE2K5q1+ST3rG3nydpE8FXSDB4UJgv9t2njMtkGMCOdeDiDQCHgYeK+O4KSKSKSKZhw8fLvUBVHh5B4fa4q4rLyIq8txn+an9kkK+mtuMoR2ZahMgDPAnrWJSQRRIcLDrvu09cs7fMYGc6+0x4BljzMnSDjLGLDTGpBlj0uLi4sq4ZID2Z8DnT1nfg6y4uJjk5GRKAllxcTEdOnTgyJEj/POf/6RXr15069aNa665xjVN9pw5c7jtttsYMGAAycnJrhlS165d65ocD6wZUxctWgRYn6pnz55N9+7dSUlJcY2EPnbsGMOHD3eNkN6+fTvFxcUkJiby448/uq7VoUMH1/1LzJkzh9tvv51BgwaRmJjIO++8w0MPPURKSgqDBw+msLAQgKysLK688kpSU1O59tprOXjwIG+//TaZmZnccsstdO3a1bW2xV//+teA8ghw9OhRBg0aRLdu3bjrrrtKnYSwKpRM7Q7w1pTL+d21l7D87iuqbJlPfwEiY28+41/ZVCV5ULVfIMEhD2jrth0PeH8M9HdMIOd66wU8KSJ7genATBGZFkA+K2d/Brx2A6z5o/U9yAEiIiKCW2+91TX/0CeffMJll11Gy5Yt6du3Lxs3bmTr1q2MGTOGJ5980nXe9u3bWbVqFV9++SWPP/54QJ/AW7ZsyZYtW7j77rv5y1/+AsDs2bPp1q0b27dvZ+7cuYwfP56IiAiGDRvGihUrANi0aROJiYm0atXK55rffvstq1at4r333uPWW2+lf//+7NixgwYNGrBq1SoKCwu59957efvtt8nKymLSpEk8+uij3HTTTaSlpfHGG2+wbds2GjRoUK48Ajz22GP07duXrVu3csMNN7jWkgiHJZv2MfqlL/nLR9bU7oA15UrEbvj7UJiXAHPj4bHmMKep59fjLa19Jdt/vrjC+ZgxtCP9klv6pK/bfYTkR1ezZFP4fkeqdghkVtbNQLKItAe+B8YA47yOWQlME5GlWG/ux40xB0XkcADnejDG/LLkZxGZA5w0xpQ9bWll7f0cigrAFFnf934ObXsG9RaTJk1i2LBhTJ8+nVdffZWJEycCkJeXx+jRozl48CAFBQW0b9/edc6wYcNo0KABDRo0oH///mRkZNCsWbNS73PjjTcC1pTbJdNsr1+/3jW19oABAzh69CjHjx9n9OjRPP7440ycOJGlS5cyevRo22sOGTKE6OhoUlJSKCoqYvBgaxK5lJQU9u7dy65du/jqq68YONCaQ6ioqIgLLrggKHlct26d65jrrruO2NjQ1OmXJSs3nx0r5/OPyLX8EBHLy47rMekfwME3rf8zZSkuhILCc9unfrCChLuoBtDrLhhYaq0qYE3SN291Ni+ty/EojhcWGWau2MFzn/yHBbemhqwNRNVuZQYHY4zD+cn9IyASeNUYs1NEpjr3vwisBoYCe4DTwMTSzgUQkRHAX4E4YJWIbDPGXBvsBwxY4i8hMsZ6kUfGWNtB1rZtW1q1asWaNWvYtGmTqxRx77338sADD3DDDTewdu1a5syZ4zrHezpsESEqKori4nNrBZw5c8bjmJJpu92n3LarihERLr/8cvbs2cPhw4d59913mTVrlm3e3acFd5/Wu2RacGMMnTt35ssvvwzod1GePLp/D4v9GfDJHH6Ru4nuUeemML8mMpPIvCDfy/EzbJgPG56DiAhI7AfjV/g9fMbQjnx98CfbLrWHTpxl5Atf6FxMqkICGudgjFltjLnYGHORMeaPzrQXnYEBZy+le5z7U4wxmaWd60xfYYyJN8bUM8a0sgsMxpg5xpi/VP4xA9C2J9y+EgY8an0PcqmhxJ133smtt97KqFGjiIyMBOD48eOutaVfe+01j+Pfe+89zpw5w9GjR1m7di09evQgISGBr7/+mrNnz3L8+HE+/fTTMu/br18/VzBau3YtLVu25LzzzkNEGDFiBA888AAdO3akRYuKrWNxySWXcPjwYVdwKCwsZOdOq7NZoNOC+8uje/oHH3xAfn6Ie+Xsz4Clt8CfEmFOLLwyEJO7gRhjBQYR6ytC7BvVgqMYih2Qs8bKw/wUyFxke+TiO3oR1zjG75VmrthB33mfMvz59XSYuYqLZq5myPx12rtJlUoX+3HXtmdQg8Lp06eJj493bT/wwAPce++9TJw40VWlBFaD780338yFF15I7969+e6771z7evbsyXXXXce+ffv4/e9/T5s2bQAYNWoUXbp0ITk5mW7dyp6/Z86cOUycOJEuXbrQsGFDjyA0evRoevTo4WrUroiYmBjefvtt7rvvPo4fP47D4WD69Ol07tyZCRMmMHXqVBo0aFBqycJfHmfPns3YsWPp3r07V155Je3ahehTcOYi+HQO/Oz5punqWeEWCYy/7hYAkfWsN3bjtXZ0RDRE1YeCsgOlp2L4cR+8fz+segAatoT+MyFtguuIzbMGMvCptbaLDAHk/XiGvB9LSpiG7EMnGPnCF0RHQFREBGccxa7nbNogiocGd9TSRh2nU3ZXsczMTH7zm9/w+eefl3nsnDlzaNy4Mb/9rQ4SD0SF/+YLB8CB0tdq9n6ZGIGI65+F/O9g00vgOAMNmsHVczzetP1aPhm+eR+KCq22iIpongQjXvL4QDNvdTavbviOAu/1RiuhecNohqRcwI3d47X9opYpbcpuLTlUoXnz5vHCCy/oimnVxfLJsONtwHet5xLG9Q8UmEjOEsXXxe1pNuyPXJp2jbUjgMZjHyNf9k1bPAK++8ytWFLke4y7YznwykCIrA+dboCRLzNjaEdmDO3Ikk37WLB2D9/n/1z+vHnf5nQhb2zaxxub9pHQvCFPj+6qQaIO0JKDqjXK9Tf/88VWb6FSGKCoGH6mHq8XDeTJIqujXVUMdgOsaq71T8OJg4H1hpIIaH+VRwP2kk37eHV9DojQOCaSbXnBXa2va3xT3p3WN6jXVFWntJJDrQ4Ol156aXh7uagqY4zhm2++KT04pM+GHcvg7Ck4W/qb5NnIxsz+eTRLi6/2SK+ywODNT3uIvQjoc5/fEs3w59e7gkRMpBAp4mpzqIhIgWVTr9DSRA1UJ4PDd999R5MmTWjRooUGiFrOGMPRo0c5ceKExxgRD8+kwPEABoZFNWRvz1kMXX8Rpws8q3U6nN+YTx64Mgg5roT9GfDJbDjwbyi0b3x2adgSxr5Z7k4W05duZdWOgxSWs92ieaMYRqXGhyd4qgqpk8GhsLCQvLw8nzEAqnaqX78+8fHxREdHe+5Inw1fPm/1HirN9c9C2gSycvO5+cUvKLZ5WVS78QKZi2DtXDh5mNLaTUqeraLGv7KJ9buPYCh77huA5LhGpD94VYXvp6pOnQwOSpXZCykiGpq19ejx0/uPn3DoxFmfQ4d3bVOpJT9Drqw2lGYJ0PeBSgUJsHpDLfpyL2cKSwlGWGNAbrismv/OlAYHVceUVT9frync+rZPdYu/cQLVPjCUCKBLbmVLEe5KG1dRolFMJB3Ob8zoHu2qV6lLARocVF2RPhs2vgRFpXTfbHkJTPOdVLHHE+kcPunbI6jGBIYSmYvg41n+B9pJJLTqBNc9HZQBn1m5+Szfksc/t33PibOld72tcb/LOkCDg6r9lk+2eiL5E9MEBj1h+6n5F3/4kJMFvm9sNbqbZuYiWP+UNbLaHz+BsqLce0H5U6N/p7VQacEhoLmVlKq29mfAi31LDwxJA2Bmnm1gGP78etvA0KJRdM1+E0ubANN3WNVILS+xxkB4O7ILHmsRtOnp353Wl7kjUmgY7f9tZVvecQY+tTYo91OhpcFB1Vz7M6wRwod22O9vEGu9OfqZ1XTe6my/n3QXju8RrFyGV9oEq3Twi5vs9xuH9Tv0M6lfeY3r1Y6v/2cIU/sl0bxhtO0xuw+foscT6UG5nwodrVZSNVNZ1UhtUmHKGr+7x7+yyXaaa6iGXVaD5fmeVmnBn6QBpU4PXhFZufnc80YWh37y7QHWrEEU22aHb5Z+pdVKqrZZPMJ/YIhpDH2mlxoY5q3OrnuBAawSRMoo//tz1sD/nB+0UgRAakIsG2deQ3JcI599P/7soO+8sqebV+GhwUHVHJmLrE+/OX7e+Bu1gpnflzkR3rLM/T5pMZHC8ruvqL2BocTIl+GOdGv0tJ2is9bU4EEMEADpD15lu+ZE3o9ntA2imtLgoGqGhQOsNy1/1SJN28Hv/lPmZbJy83HYDH+e1Kd93ZkbqG1PeOhbqz3Gn9UPBj1AbJ41kGYNfCeC3n34FMOfXx/Ue6nK0+Cgqr/lk/0P7ippdP6Nn0ZpN9OXbuWmF77gpzOeU2n0S25ZN+cDSpvgP0AUO6xgvHBAUG+5bfa1xDer75ued1wDRDWjwUFVb/sz4Kvl9vuufxYe3hvQiN++8z7l3W0HPOYGSmzRkLkjUlh8R69g5LRmSptQejXTgSxrao4gWj/jats2iG15x5m3Ojuo91IVp8FBVV+Zi+DVa32X2wSrYTXAaSCGP7/ebYlMiwBPjepa+9sYAlFSzZTkp5Rw6oeglyDSH7yK1ufV80l/cV1OUO+jKk6Dg6qeMhfBqulgbCZ4Sxpgv5KajSWb9tmOZRjWtU3daWMI1PgV/nszHciyOgME0YJbUm3TtXqpetDgoKqf5ZOt+m67MTgpowLui5+Vm8+jK3zbIuKb1dc5fvwZ+bJVXedvRPW8xKDdKjUh1m/1UlZuIIsaqVAKKDiIyGAR2SUie0Rkhs1+EZHnnPu3i0j3ss4VkZtFZKeIFItImlv6QBHJEpEdzu/BLc+q6u3PF9uPYWjY0nrTCrDEAPDAW9t81h9o3aQe62dcbXu8ckqbAJM+st93Jt8aZxIk/tZ9WL4lL2j3UBVTZnAQkUhgATAE6ASMFZFOXocNAZKdX1OAFwI49yvgRmCd17WOAL8yxqQAtwOvl/+xVI20cID9mgQRUdaKZuWYanre6mxyj532SV9wq31VhvLStqf/huqcNdZ8VkGak6lrfFOftD0/+JlVVlWZQEoOPYE9xpgcY0wBsBQY5nXMMGCxsWwEmonIBaWda4zJNsb4dFo3xmw1xhxwbu4E6ouIb8uVql3SZ9t3V22dAhM/KPf00nYD3ab2S9J2hvIoaaiub/M7O7TDmpMpCKWId6f19aleytibz/hXNlX62qriAgkOFwLur7Q8Z1ogxwRybmlGAluNMb4Ts6jaI302bJjvm96oFUxdX+7AMG91NsdOF3qk1dmxDMEwY6//nkw5a6w2okpKf/AqLvMqQazbfUS7toZRIMFBbNK8q3L9HRPIufY3FekM/Am4y8/+KSKSKSKZhw8fDuSSqjpaPhm+eM43vU1qQCOevWXl5vPS557dIeOb1a/bYxmCobSeTDuWBaWKaXQP327F7277vtLXVRUTSHDIA9q6bccDBwI8JpBzfYhIPLACGG+M+dbuGGPMQmNMmjEmLS4ursyHUNXQMynWG4t3d9UyJs4rzawVO3w6Of26f3IFM6g8jHzZWhvCzisDKx0gxvVqR0Lzhh5pBY7S16pWoRNIcNgMJItIexGJAcYAK72OWQmMd/Za6g0cN8YcDPBcDyLSDFgFPGKM2VDO51E1xcIBcNxrlTKJsAJDGRPn+TP8+fVkH/JsyOyZGKsD3YJpWoZVqrOzclqlA8TTo7t6bB87Xcj0pVsrdU1VMWUGB2OMA5gGfARkA8uMMTtFZKqITHUethrIAfYALwO/Lu1cABEZISJ5wOXAKhEp6Ts3DegA/F5Etjm/zg/O46pqwV/j8y9uqnBg8Ldwz8NDtJ0h6KassSY69HZ4F7wyqFIT9qUmxFIvyvNt6Z/by6xsUCGgi/2oquWv8blpu4Amz/Ony5yPfCbU0wXtQ6y0xYPuSC93R4ISaU+kc+RkgUfa1H5J2qEgBHSxH1U9PN/TPjC0Sa1UYMjKzefnQs/5l2IiRQNDqE3L8D+r6wrbfiQBeWCgb7vGkox9NkeqUNLgoKrG4hH2nzIr0fgMVmAY9dIXFBZ5loAn9Wlf4WuqckibYP0NvR3LqfBkfXYN0z+dcbBkkwaIqqTBQYVe5iL4bq1veiUan0v86YNsirw6tPRIjNUqiKo08DH7cRCVmKzPu2Ea4Jn0Uta/VkGnwUGFVuYi5yR6Xu/gLS+pdGDIys0nY6/nBG0CzNBG6Ko3foV9gDiyq0KD5FITYjmvvueqcYdPFuiEfFVIg4MKnf0Z8IVXnbREWG8i0yo/aGrya5t90u7SKTLCZ/wK+26uO5ZVKECM6+nbI2qWzSy7KjQ0OKjQyFwEfx8Cx77zTL/umYCn3C5NVm6+zxQZAlqdFG5T1tgPlNuxrNzzMM0Y2pGOrZt4pGUfOqFTalQRDQ4q+PZnwKoHrHWIS2ZRaZ5k9Wwpx8yqpbH7BOk9N48Kk2kZ9lNtVGAepidGpPikvaE9l6qEBgcVXPsz4O07PJf2lAgY8VLQAsP0pVt9RkK3aBjNu9P6BuX6KghGvmwfIHYss8a6BCg1IZZmDaI90k6ccWjbQxXQ4KCCJ322NceO97QYlwyu8IAob1m5+by7zXfE7MLbewTl+iqIRr5s30i9YX65RlH3bN/cJ+2lz2ynXFNBpMFBBUfmIvsBboh9P/gK2phz1CetY+sm2ghdXY1fAY1b+aa/f3/AVUx3XXmRT9qWfVpyCDUNDqry9mfAhz6rx1r63B+0UgNAbMMYIt0mgo8Q+3ppVY1cNdM+PcAqptSEWOKb1fdIO6LdWkNOg4OqnP0Z8OpgcPzsuy9lVKXHMribtzqbWe/uoNhAZIQwqFMr/jH1Ci01VHdpE/yvBbHpxYAuYTft+otatRRSGhxU5bz/gGfjc4k2qVadc5As2bSPF9flUGys/k/FxYbL2jbTwFBT+FsLwnEmoFHU43q1I65xjEfaNq1aCikNDqri9mfADzaDkpIGVGq+JDt/+tC3b3vvpBZBvYcKMX9rQRzZFdA8TN3aeX4QOHyyQOdbCiENDqpi9mfAynt901NGBWWQm7t5q7M5/rPndNwtG8doqaEmmrLG/zxMz5TedmTXMP3qhu9sjlTBoMFBld/yyVaX1cPfeKY3bhXUqqQSy7L2+6R1baeBoSEl0YkAACAASURBVMYav8K+iun4PpiX6Pe01IRYeiR6/t1z/ntSG6ZDRIODKp/FI6xeJnb89UqppIJC33WEp9p8ilQ1yLQMaNLGN/1MfqlVTDOGdPR40yoG3tmSF/TsKQ0OqjwyF1lTIHhr0iaoU2O4W7JpHycLPBu8eybGapVSbTDqNfv0A1l+16JOTYjliREpREac68/81uZ9WnoIAQ0OKnCfzvFNk0jrRR6CwADWC9+brgtdS7Ttad/+APDJHL+njevVjgGXnltW3lEMy7X0EHQaHFRgnkmBn70+nTWIhUkfBnWQm7us3Hy+PviTR9pUnZK7dhm/Aurb/D1zN5RavSRe20dOnA1uvpQGBxWAxSN850sCGLcsZIEBrPlz3Jf/HNSplU7JXRvN2GvfQH0gy2+AiGtSr9RtVXkaHFTplk+2b2dIGhDSwDBvdTYff/2DazsqUmy7MqpaYlqG/RxcB7Jsp9i4sXs8MVERCFYp4oOvDuo6D0GmwUH5lz7bvmdSZL2gj2Vwt2TTPl5al+OR1vmC87Q6qbYb+Jj9ILkN830aqFMTYpnzq86ANWL+2KlCXlyXowEiiAIKDiIyWER2icgeEfGZYU0szzn3bxeR7mWdKyI3i8hOESkWkTSv6z3iPH6XiFxbmQdUlfDFX+3ThzwZsltm5eYz690dGK/00T18l4xUtdCUNUCkb/r/u8knKf90gc//k3e3fR+SbNVFZQYHEYkEFgBDgE7AWBHp5HXYECDZ+TUFeCGAc78CbgTWed2vEzAG6AwMBv7PeR1VleYl2s+ZlDIqZD2TAO55I4tir1f81H5JjOulwaHOSBnpm3b2uM8cTDp9SmgFUnLoCewxxuQYYwqApcAwr2OGAYuNZSPQTEQuKO1cY0y2MWaXzf2GAUuNMWeNMd8Be5zXUVVl8QhrMJK3lFEhGQFdYvrSrRz6ybPXSYfzG2sjdF3jb5K+I7s81qFOTYilQ1wjj0MO/XRWxzwESSDB4ULAff6CPGdaIMcEcm5F7oeITBGRTBHJPHz4cBmXVAHzN9CtabuQBgaAj3Ye8kmb1Kd9SO+pqqlpGdb/OW85azx6ME3qm+RziK4SFxyBBAfvLsWAT1Wfv2MCObci98MYs9AYk2aMSYuLiyvjkiog6bOtFbq81Y+F39jMvhpE81Znc9bhOU1G43qRWp1Ul/1mh/8urs4eTON6tfNZCOjrA8erIne1XiDBIQ9o67YdD3gv4uvvmEDOrcj9VLD5W+YzaYDVDz2E5q3Odq3V4G7mUO+mLVXnTMuAqPq+6V/+n+tH74WA8n48o1N5B0EgwWEzkCwi7UUkBquxeKXXMSuB8c5eS72B48aYgwGe620lMEZE6olIe6xGbvuJVlTwrJ3rmxbVMKRdVsHqtrroy70eafWjI5g7IkVLDcrSa6pvWnGBq4F6XK92XOhVenh1fY7vOapcygwOxhgHMA34CMgGlhljdorIVBEp+autBnKwGo9fBn5d2rkAIjJCRPKAy4FVIvKR85ydwDLga+BD4B5j7LrNqKBZPhlO/uCb3vH6kN52yaZ9zFyxgzNes65OuDxRA4M6Z+Bj9nMwHdll/d8FOrdp6rHr2yOntGG6ksSYspoAqr+0tDSTmZkZ7mzUTOmz7auT2qQGfTU3b7/4w4ceM65GRQh39m2vvZOUvYUDrPYGb9c/S1bcMG5+8QuPqskeibH8Y+oVVZe/GkhEsowxaXb7dIR0XeavnaFViiswTF+6lYtnrab9jFUkzlhF8qOrKz0KNSs3n95zP/GZijuxRUMNDMq/KWvsG6jf/w2pEbu5pmMrj+TNe/O17aEStORQV+3PsFZz82Ikgpcu+j/ePNCa3GOn/Z4eIdCmaX1+3T+5XFVAJY3PdrSdQQVk+WTfaV3qN2Vv94fp/69E3N/SLotvynvT+lZt/mqQ0koOGhzqque6wbFzb9IGOEEjJpz9HVvMxeW6VHQkgBABOIoNKRc25V23F2RWbj4bc47y5qZc8n48Y3uNrvGe5yhVqgW9fJepBVYnzODXu7q4tiME/jH1Cp2Xyw8NDsqTV92tcf7zSOEdLC2+Oii3iBRo3bQ+FzZrwJZ9+Th8V/p00cCgyi1zkf2YnOiGXBWzhL1Hz5V6eybGskzbHmxpm4M6Z/EI28CwtiglaIEBoMjA9z+eIWOv/8AQ36w+c0ekaGBQ5Zc2wX6K78LTzD01xyMp22vBKBUYDQ51yf4M36kxDOQUt2ai4xG/p7VuUo/ld1/B8ruvoGdiLI1iIqkXKTSpV/H5EJPjGrF+xtXaxqAqzk8X18vZxkORS1zbJ84WacN0BUSFOwOqCq16wGPTGCgGfufwHGQUAdzQtQ3zx3TzuYR38Xze6mwWb9xLYZGhqMhQSu0RAkRHCr2TWrD4jl4Vewal3I1fYQ2GO3JuDk8B7op+n0+K01ztZ29t3qcfRMpJg0NdsXgEHNrhmvDKGKtK6dHCO1wvoGYNonlo8KXlehHNGNrRo/vpvNXZLMvcTzHQKCaSzm2actUl55N/uoDeSS20YVAF37QMnw4WEcBDkUsZ4/gDAPWitJKkvLRBui7warwrCQwz3Rqgk+Makf7gVWHJnlKVtj8DXr0WjFV2LWlL21KUxEjHE0RFwFt3aa8lb9ogXddtXeya1rbks8CKoj4aGFTt0bYnXPeMa1Oc/3SPzGF51CwcxbB8S17YslcTaXCozfZbxe3i77d4THqeU9yaBx33AFY3Ug0Mqlbw6sFUMvd/98gcxkR8ys7vj+t8S+WgwaG2KhkBfSwHcUYGY6CQSFcDdL/kltqNVNUuAx+z5gXDuaCMM0LcHfke/847zuiXvtAAESANDrXVJ3NcP5Z8gio0kYwp+D1bzMX0S26pPYZU7TRlDTRpc27VMIF2EUcYE/EpjmKYtSK0C1fVFhocaqP02ZC7waedYUNxJw0Mqm4Y9RolH4tKgsRDUW8xJuJTsg+d0HEPAdDgUNs4Z1p177IKVq+NO4seYe6IFA0MqvZr2xP6WD30Sj4kxcpJ/jf6FZ6KWqCLAQVAg0Ntsj8D3r/fZ/HuTUWXMtLxBBe3aqIDgVTdMfAxuP5ZiiIbAOfaH26M3MCNx17W0kMZNDjUJu/cBZwLDCXjGZ4sGgPAEyNSwpMvpcIlbQJRzRNcLwpXA3XU+3zycVkrFtdtGhxqi+WTIT/HdjzD7piOLL9bBwCpOqrX3cC510RJgPhLwR+151IpNDjUBumzYccyn3aGtUUpPOi4h0eGdtLAoOqutAlIn+kY8FgIKDbiNN+8/puwZau60+BQ03k1QJc4VtyIiY5H6BrfVNsZlBr4GD/GtAasAFFSehhXuILVf//fMGas+tLgUNN9PMsjMJR8Mvpz0RiS4xrpIDelnJo/uosCiXZtlwSIS3JfD1OOqjcNDjVZ+mxMwQlXX72SwPCC43r+GT1Ip8VQyku96/8C4lm9lGS+J+/Jy8OXqWoqoOAgIoNFZJeI7BGRGTb7RUSec+7fLiLdyzpXRJqLSLqI7HZ+j3WmR4vIayKyQ0SyRcT/KjR13fZl4Cwil/xnf6eoD08WjePRoZ3CmzelqqOS9gfxrF668NTX5D/pu35JXVZmcBCRSGABMAToBIwVEe93niFAsvNrCvBCAOfOAD41xiQDnzq3AW4G6hljUoBU4C4RSazg89VOmYvgmRTMiQOA50C3Bx33MLVfkrYzKOXPwMfYd8VcwDNANDudY3XuUEBgJYeewB5jTI4xpgBYCgzzOmYYsNhYNgLNROSCMs4dBrzm/Pk1YLjzZwM0EpEooAFQAOgisCWcazOY4/tcpQaATcXWQLep/ZI8Ft9RSvlKHHQP/5YOgFuAMFC8Yb41mFQFFBwuBPa7bec50wI5prRzWxljDgI4v5/vTH8bOAUcBPYBfzHGHPPOlIhMEZFMEck8fPhwAI9RS3w6x9UAXVKdVIzwpGMM/ZJbamBQKkAf9vx/7CpqA5wLEGLgp/dnhjln1UMgwUFs0ryXj/N3TCDneusJFAFtgPbAgyKS5HMRYxYaY9KMMWlxcXFlXLKWSJ+N+TnfpwH6Jcd1tLvsKp0zSalymDG0I3PbL2JT8aUe6U1+2GytS13HBRIc8oC2btvxwIEAjynt3B+cVU84v//XmT4O+NAYU2iM+S+wAbBdxq5OWTwCs2G+TwP02qIUfuozi/ljtDFNqfJafEcv/hk3hWLEo3rJHNlV5wNEIMFhM5AsIu1FJAYYA3hPSrISGO/stdQbOO6sKirt3JXA7c6fbwfec/68DxjgvFYjoDfwTQWfr3ZYPAKTs8Y2MDzSYI5WJSlVCTcOu5GXiq4D8A0QdbiBuszgYIxxANOAj4BsYJkxZqeITBWRqc7DVgM5wB7gZeDXpZ3rPGceMFBEdgMDndtg9W5qDHyFFVz+bozZXtkHrbH2Z3gEhhKbii5louMRFtyaGr68KVULpCbE8tMVs3zaH1wN1JmLwpq/cBFjymoCqP7S0tJMZmZmuLMREqee7U3DY9muwFAy0+pNBXO4afhI7bKqVJBc9thHpDvuIC7ixLnXGyASCZM+tNaIqGVEJMsYY1ttryOkq7G9Hy+gQX62a7skMMwsvEMDg1JB9vDgjkx1PIiDCOu1ZkomsiyCVQ+EO3tVToNDNbX34wW03TAT8WpnWFHUhy7DpmtgUCrIxvVqR7vLrmJ0wR/YXezWW9+AObQDnqlb66FocKiO9mcQ/8VMIjgXGAzW1Bif/+KPGhiUCpH5Y7rR85eDmeGY7Jri29VAfXwfLB4R7ixWGQ0O1dB3r91FpFcD9J7iC/n8F3/ULqtKhdiMoR3p+cvBfF2U4EpzNVDnrKkzI6g1OFQzez9eQGLhucXPS6qT9iTdpoFBqSoyY2hHnq43xWeBIDHAK4PqRIDQ4FCNbH9uNO02WEP33dsZsqK7MXSiTk6rVFW6ZtANfFZktTO4ptcQMBjY8GyYcxd6GhyqiS9enEbK0Q895kwC2NBoAGmz1oYza0rVSeN6teMeeZRDRU0BtxKEAfPN+9a67bWYBodqotPB5YBnO8MpU4++D60IU46UUq/f2ZvLHS/wfVFzV5qrgXrHslo9glqDQzVwbE4bmnIawNW/GmBb4z5hzJVSKjUhlqn9kriv6D4KiTw3/sEZIBwZr4Q7iyGjwSHMjs25kFhzCsBjFPSGRgO01KBUNVDSe2lMwe/Z6ezBVPIBLrLwRK2tXtLgEEZfP3kNseYk4NnOsLHNbRoYlKpGZgztyKNTJzA+5i/8aBoCtb96SYNDmGz6x1N0PLUZ8GxnOF0czRVTnw9TrpRS/qQmxPLbay/lzaIBgNcAuQ3zYeGA8GYwyDQ4hMkvdlqT0LpXJRUDG/q8Gr5MKaVKNa5XOzYm3cc7RVZ7oPsa1BzIqlUBQoNDGBx8LJmGpgA41wBdDLx9wYMMGnxDeDOnlCrV4jt68fkv/sgLjusBtwn6wAoQtWSKbw0OVWzna9NpXWwteudeavj08tcZPfUPYcyZUipQ88d046ML7mZLkbWCsUeAWPu/4cxa0GhwqEI/LvwVnXL+Dng2QH8bcaGWGJSqYUb3aMdIxxMeAQID5uShWjFBnwaHqpI+m6YH1gGegSG3qCWbh3wYxowppSpiXK92dI1vykjHE5wx0YBb+0POmhrfg0mDQxU5/eXffNaA/rsM44tf/Uun4Faqhnp3Wl+6xjfl2+I2rjRX9VINX2JUg0Oo7c+A96cTWXzGI/m0ieGyic9qYFCqhnt3Wl+ejJpMMec++Ll6p3/wUJhyVXkaHEJpfwa8Ohgy/04MDuDcf54lDCY1ITaMmVNKBctXkZfwaOEdwLnXuAEoOltj2x80OITSJ3PAFAHnipr5pjEvOK7naO+Z4cyZUiqIRqW2ZWnx1ax1n+IbZ4DIWVMjxz9ocAiV9NmQuwHXLL8GDMJDMTP5qc8sZgztGNbsKaWCZ8bQjvRLbslExyPkFrUEbMY/PN8znFkst4CCg4gMFpFdIrJHRGbY7BcRec65f7uIdC/rXBFpLiLpIrLb+T3WbV8XEflSRHaKyA4RqV/ZB61S6bNhw3wMzv8czgjxkuM6uvW5VgODUrXQ4jt6Ed+sPlc5nvMpQQBwZFeNqmIqMziISCSwABgCdALGikgnr8OGAMnOrynACwGcOwP41BiTDHzq3EZEooD/B0w1xnQGrgIKK/6IVWx/htVLAbf/FMAZE82fi8fRO6lFePKllAq59TOuJioCJjoe4fvi5r4H1KAuroGUHHoCe4wxOcaYAmApMMzrmGHAYmPZCDQTkQvKOHcY8Jrz59eA4c6fBwHbjTH/BjDGHDXGWXFfE7w51mOzpNSwurgnPRJitRFaqVruzr7WoLj7HPcBnmtQA7D5b1Wco4oJJDhcCOx3285zpgVyTGnntjLGHARwfj/fmX4xYETkIxHZIiI1py/Y/gw4fcS1WfKfYldRGx503MPDQ7Q6SanabsbQjvRMjGWLudhjgj6XgpM1onopkOAgNmnesdDfMYGc6y0K6Avc4vw+QkSu9smUyBQRyRSRzMOHD5dxySqQPhuWjHJtlvxnOFTUlMGOv9AgOlJLDUrVESUfBB903OPR/uCSs6baLxIUSHDIA9q6bccDBwI8prRzf3BWPeH8/l+3a31mjDlijDkNrAa648UYs9AYk2aMSYuLiwvgMUJo4QCrneHnfAzgMFYE3FXUhssdLwBwbedWYc2iUqrqpCbEktDcWhRoouMR+xLEjmXVegR1IMFhM5AsIu1FJAYYA6z0OmYlMN7Za6k3cNxZVVTauSuB250/3w685/z5I6CLiDR0Nk5fCXxdwecLveWTrW5qbn4yjbmpYA6DHX8BIDpSmD+mWzhyp5QKk6dHd3X97FGCcD9o1fRqGyDKDA7GGAcwDetNOxtYZozZKSJTRWSq87DVQA6wB3gZ+HVp5zrPmQcMFJHdwEDnNsaYfOBprMCyDdhijFkVhGcNvsxFVvR3Z2Bp0VVsMRe7ku7o075q86WUCrvUhFj6Jbd0bU90PGKtQe0eHYyB9++vlj2YxPg0pdc8aWlpJjMzs2pvuj8DXrkWa5mecxzF0KFgiWs7JlL4zx+HVm3elFLVRsrsDzlx1upw2V3+w9KY/yFainwbZFNGwciXqzRvIpJljEmz26cjpCtq1QN4BwYDrC9O8UibpKUGpeq0RZN6uX7eYi5mTMHvOS5NfQ/csaxalSA0OFTEwgFwaIdHkgH2RF/CRMcjrrQeibE6GlqpOi41IZZBnc51SNliLmbe2Zvsu21umG/VSlQDGhzKK322TwO0AWY57mTgidkYrF9q/egIZui4BqUUcNeVF3lsLy2+moyIrvYHO2dYCDcNDuWxfDJ88ZxP8rp6A3jDcW7WxaTzG/PGnb11XINSCrBKDz0TPd8PRp9+iO9bXOF78DerqsUkfRocAvV8T6tO0Hi2M5xo2ZXbj9/pkdYoRge8KaU8PTykIxFerdD9Dt7LwZSpvgcf2RX2ab41OARi8Qjrj+VBoM90HjzvKZ/DR/fQ1d2UUp5SE2K5uqPnYNiiYsPzEbfC9c/6nnAgK6xjIDQ4lCVzkTXU3VvjVjDwMdbt9py6IyZSdOlPpZStqVdeRJTXu+4/MveTFTcMGrb0PSGMYyA0OJRl7Vz79KseISs3nzOFntVMxTV/2IhSKkRSE2J5664ruCz+XFfWwiLDxpyjMPZN+5M2zA9LgNDg4E/mInjiAjj5g+++PtMhbQIvffatz64rLtL1GpRS/qUmxHK527ouBkjfeQja9oQ70u1LEBvmV/lEfRoc7KTPtopzjtO++65/FgY+RlZuPp9mewaO1ufVY/EdvXzPUUopNzsP/uSxvS3vOPNWZ1sBYsDv7U+q4on6NDh4c1vJzUfKKEibAMDGnKMUuVUhRUYIC25JDX3+lFI13pBfXOCT9vrGXOuHtAnWe42dz307wISKBgdvK+6yT295ice8Jyd+9ly5dHLf9tp9VSkVkHG92tG0YZRH2qmCIpZs2mdtjHwZkmy6sh7fB892q5JR1Boc3C0eAcdyfNPbpMK0c3+MrNx8Fn5+7jgBmjSIroIMKqVqi4ev9Z1B4a3N+85tjF9hHyDyc+CVQSEPEBocSiyfbN9lNWkATPFMf3DZNo9eSSLQO0kbopVSgRvXqx1d4z0n4Nvx/XGycvPPJYxfYXWA8WF81qsPNg0OYI1E9F6XAazAMH6FR9K81dnsPerZUJ3UspFWKSmlys27xqHYWO2ZHgY+Zl+COH0E/nyxb3qQaHBYOMBnIj3ANjAAfLjzkE/apL5JociZUqqWs2uYjm0Y43vg+BVW9ba3Uz+EbJqNuh0c9mfYB4aUUbaBAWBw59Ye2/2SW+qIaKVUhYzr1c5jOm8Bdh44bn/wlDVWxxhvB7JCMkiubgeHfy/xTWuTWupqTDOGdmRqvyQSWzRkar8kHdeglKqUu668iBjnnBoGWLJpnzXmwc60DGjUyjf9i+eC3kBdt4OD90J9Tdr4ND67m7c6m6v+/C8A1v6uvy7ko5SqtNSEWG5KjXdtG+DFdTnnurV6+91/fKuYjIG9nwc1X3U7OFw2FiLrAWJ9H/Wa30Pnrc7mxXU57D16mhfX5fiP7EopVU4ju8f7pD2T7j0TtJspa6xeTBIBCETVh8RfBjVPUWUfUou17QkT3rcibuIvrW0/3t32vc+2lhyUUsGQmhBLdKRQ6DbtwnGvgbY+Bj4Gl14X0PtXRdTt4ADWLzSAX+qpsw6P7XpRkaHKkVKqDuoQ15jsQydc2xfFNS77pADfvyqiblcrBWje6mxOnC3ySGvWUEdEK6WC54kRKa6V4gTonhDrOSCuigUUHERksIjsEpE9IjLDZr+IyHPO/dtFpHtZ54pIcxFJF5Hdzu+xXtdsJyInReS3lXnAYPCuUgJd7U0pFVypCbH8Y+oVjOvVjqhIYcmmfYx9eWPYAkSZwUFEIoEFwBCgEzBWRDp5HTYESHZ+TQFeCODcGcCnxphk4FPntrtngA8q8EzB57WAT0Lzhjq2QSkVdKkJsQjWAkAGKHAU886WvLDkJZCSQ09gjzEmxxhTACwFhnkdMwxYbCwbgWYickEZ5w4DSroHvQYML7mYiAwHcoCdFXyuoJm3OptDJ856pD09umuYcqOUqu12/3DCY3tLdS05ABcC+92285xpgRxT2rmtjDEHAZzfzwcQkUbAw8BjpWVKRKaISKaIZB4+fLi0QyvFu0qp9Xn1dB4lpVTInHV4Lj2cfehEWKqWAgkOYpPmvVKyv2MCOdfbY8AzxpiTpR1kjFlojEkzxqTFxcWVccmKycrN579epYbhXb3jolJKBY9de+asFTuqPB+BBIc8oK3bdjxwIMBjSjv3B2fVE87v/3Wm9wKeFJG9wHRgpohMCyCfQffOljyPqbk7tm6iYxuUUiE1rlc7Loxt4JGWfeiE/xHTIRJIcNgMJItIexGJAcYAK72OWQmMd/Za6g0cd1YVlXbuSuB258+3A+8BGGN+aYxJNMYkAvOBucaY5yv+iBXnXdcX37xhOLKhlKpj7rmqg0/a05+UMmI6BMoMDsYYBzAN+AjIBpYZY3aKyFQRmeo8bDVWA/Ie4GXg16Wd6zxnHjBQRHYDA53b1caSTfs8BqQAnN+kXphyo5SqS8b1auczlurIiYIqLT0ENELaGLMaKwC4p73o9rMB7gn0XGf6UeDqMu47J5D8hcJ8rygtwI02858opVQo9EhsTvrXP3ikffDVwSrrRq8jpG3MW53Nf08UeKRd2rqJ9lJSSlWZqVde5NOjp/MF51XZ/TU42LAbEd1dA4NSqgqlJsQysJPn2g05R05V2f01OHix676qVUpKqXBo6dXO+Wn2D1U25kGDg5eNOUc9uq8CjO3VTquUlFJVbmT3eCLd6paKDVU2nYYGBy8nvOZQj4ywX4hDKaVCLTUhlv8ZnkKkc7pWA7yVub9KSg8aHLx8ku3ZO6B9y8ZaalBKhc24Xu24+tLzXduOIsNLn30b8vtqcHCTlZvv0+CT1LJRmHKjlFKW/NOevSfTvw5924MGBzd/+iDbo71BBO668qLwZUgppbCm7nZnsNpHQ0mDg9OSTfvI2OsZiQd2bKVVSkqpsLObjC9956GQ3lODg9Nbmz2HpQtaalBKVQ/jerWjQ5xnFfe2vOPMW50dsntqcHCqF+X5q+iRGKulBqVUtTGpb5JPmt2A3WDR4IDVEL1t/4+u7cgIeHiITs2tlKo+xvVqR2uvQXGxDWNCdj8NDsDyLXkUFlkt0YJVv6elBqVUdXO113QapwuKQnavOh8csnLzeXPTPtfydNGRooPelFLV0o3d4z0m48s9dprxr2wKyb3qfHCYtWKHx7qlbZo10FKDUqpaSk2IpVG9SI+0dbuPhGTMQ50PDvvzT3tsHzl51s+RSikVftd0bOWTFooxD3U+OHRv51lKsPvFK6VUdTF/TDe6xjf1SPOeEy4Y6nRwyMrN94i4/ZJbMn9MtzDmSCmlyjawc2uPtoe/rf8u6FVLdTo4LN+SR0HRuRaH+OYNw5gbpZQKTO+kFq6ZWgGKi03Qq5bqdHDwXoLPe1sppaqj1IRYHh/2C6IihAggJjqC3kktgnqPqKBerYa5sXs8/8jKo9BRTHRUhK72ppSqMcb1asclrZuwMecovZNaBL2XpRhjyj6qmktLSzOZmZkVOrek3SEUv1yllKrORCTLGJNmty+gaiURGSwiu0Rkj4jMsNkvIvKcc/92Eele1rki0lxE0kVkt/N7rDN9oIhkicgO5/cB5X/kwKUmxHJP/w4aGJRSyk2ZwUFEIoEFwBCgEzBWRDp5HTYESHZ+TQFeCODcGcCnxphk4FPnNsAR4FfGmBTgduD1Cj+dUkqpCgmk5NAT2GOMyTHGFABLgWFexwwDFhvLRqCZiFxQxrnDgNecP78GDAcwxmw1xhxwpu8E6ouI52xTSimlQiqQ4HAhsN9tyHhM9gAABEtJREFUO8+ZFsgxpZ3byhhzEMD5/Xx8jQS2GmN8hi2LyBQRyRSRzMOHDwfwGEoppQIVSHCw6+Hp3Yrt75hAzrW/qUhn4E/AXXb7jTELjTFpxpi0uLi4QC6plFIqQIEEhzygrdt2PHAgwGNKO/cHZ9UTzu//LTlIROKBFcB4Y8y3AeRRKaVUEAUSHDYDySLSXkRigDHASq9jVgLjnb2WegPHnVVFpZ27EqvBGef39wBEpBmwCnjEGLOhEs+mlFKqggIa5yAiQ4H5QCTwqjHmjyIyFcAY86KICPA8MBg4DUw0xmT6O9eZ3gJYBrQD9gE3G2OOicgs4BFgt1sWBhlj/osfInIYyC3Xk5/TEquHVG2gz1I91ZZnqS3PAfosJRKMMbb18rViEFxliEimv0EgNY0+S/VUW56ltjwH6LMEok7PraSUUsqeBgellFI+NDjAwnBnIIj0Waqn2vIsteU5QJ+lTHW+zUEppZQvLTkopZTyocFBKaWUjzobHETkZhHZKSLFIpLmlp4oIj+LyDbn14vhzGdZ/D2Hc98jzqnSd4nIteHKY0WIyBwR+d7t7zA03Hkqr7Kmuq9JRGSvcxr9bSJSscVTwkREXhWR/4rIV25ptksGVHd+niUkr5U6GxyAr4AbgXU2+741xnR1fk2t4nyVl+1zOKdGHwN0xhqc+H/OKdRrkmfc/g6rw52Z8ghwqvuapr/zb1HTxgcswnoNuPO3ZEB1twjfZ4EQvFbqbHAwxmQbY3aFOx+VVcpzDAOWGmPOGmO+A/ZgTaGuqkYgU92rKmCMWQcc80q2XTKguvPzLCFRZ4NDGdqLyFYR+UxEfhnuzFRQIFOtV3fTnCsLvlpTiv1uasPv350BPnauzjgl3JkJgkCWDKhJgv5aqdXBQUQ+EZGvbL5K+wR3EGhnjOkGPAAsEZHzqibH9ir4HBWeLr2qlPFcLwAXAV2x/iZPhTWz5Vftf//l1McY0x2rmuweEekX7gwpl5C8VqKCcZHqyhhzTQXOOQucdf6cJSLfAhcDYWuEq8hzENhU62EV6HOJyMvA+yHOTrBV+99/eZSszmiM+a+IrMCqNrNrr6spfhCRC4wxB72XDKhpjDE/lPwczNdKrS45VISIxJU03IpIEta62DnhzVWFrATGiEg9EWmP9RwZYc5TwErW+nAagdXwXpMEMtV9jSAijUSkScnPwCBq3t/Dm+2SATVRqF4rtbrkUBoRGQH8FYgDVonINmPMtUA/4HERcQBFwFRjTJU0AFWEv+cwxuwUkWXA14ADuMcYUxTOvJbTkyLSFasqZi9+VgSsrowxDhGZBnzEuenqd4Y5WxXVClhhzcxPFLDEGPNheLMUOBF5E7gKaCkiecBsYB6wTETuwLlkQPhyGDg/z3JVKF4rOn2GUkopH1qtpJRSyocGB6WUUj40OCillPKhwUEppZQPDQ5KKaV8aHBQSinlQ4ODUkopH/8fozPMk6dgP7UAAAAASUVORK5CYII=","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["plt.plot(X.flatten(),maxeigvec.real,'.',label='MPD method')\n","plt.plot(X.flatten(),Lyapunov,'.',label='Lyapunov method')\n","plt.title('N = 100, M = 5')\n","plt.legend()\n","plt.show()"]},{"cell_type":"markdown","metadata":{"id":"sJmIXjP8NI3P"},"source":["Now consider the problem of minimizing the infinite-horizon discounted cost\n","$$c = \\sum_{t=0}^\\infty \\gamma^t C(x_t,u_t),$$\n","where $\\gamma\\in(0,1]$ is a *discount factor*.\n","\n","The un-discounted case $\\gamma = 1$ corresponds exactly to LQR; unfortunately, policy iteration won't converge when $\\gamma = 1$!"]},{"cell_type":"markdown","metadata":{"id":"BbyWp2KYNI3Q"},"source":["c. *Compare the optimal value and control for the MDP obtained using value or policy iteration with the optimal value and control for the original LQR problem.  Discuss how the comparison varies with respect to $|X|$, $|U|$, and $\\gamma$.*"]},{"cell_type":"code","execution_count":427,"metadata":{},"outputs":[],"source":["gamma = 0.9\n","n = m = 1\n","N = 100\n","M = 5\n","X,U,P,C = MDP_LQR(N,M,n,m)   "]},{"cell_type":"code","execution_count":428,"metadata":{},"outputs":[],"source":["#(1) find the optimal value and optimal policy from LQR\n","Y = la.solve_discrete_are(A,B,Q,R)\n","K, _, _ = ctrl.dlqr(A, B, Q, R)\n","val_Opt = 0.5*np.diag(np.dot( np.dot(X,Y), X.T )) #optimal value = 1/2*X.T P X, P (i.e. Y here) satisfies DARE\n","pi_Opt = -np.dot(K,X.T).flatten() #optimal policy u = -Kx\n"]},{"cell_type":"code","execution_count":429,"metadata":{},"outputs":[],"source":["#(2) policy iteration and value iteration\n","PI0 = np.abs(np.random.randn(N,M)) #initial policy\n","PI0 = PI0 / np.sum(PI0,axis=1)[:,np.newaxis]\n","X_index = list(range(N))\n","U_index = list(range(M))\n","PI_policy, v_policy = policy_iteration(gamma,X_index,U_index,P,C,PI0)\n","PI_value, v_value = value_iteration(gamma,X_index,U_index,P,C,PI0,maxiter=100)"]},{"cell_type":"code","execution_count":433,"metadata":{},"outputs":[{"data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAegAAAHSCAYAAAAnsVjHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZyT9bn//9eVzAy4ICAIoiwDLqiAjDIFrKVqtWBbFTcqalu1oqf9as851XPU2lNrF/1ZbXvs6TmtR7SCrRviAlitiIriUZYMRgWUisgyQhlAQBCYTJLr90eSIZNJZmG2zMz7+XjMI7nv3MvnvpPJlc9u7o6IiIjkl0BbJ0BERERqU4AWERHJQwrQIiIieUgBWkREJA8pQIuIiOQhBWgREZE8VNDWCUjXu3dvLy4ubutkiIiItJqysrIt7n5Y5vq8CtDFxcWEQqG2ToaIiEirMbO12dariFtERCQPKUCLiIjkIQVoERGRPKQALSIikocUoEVERPKQArSIiEgDhCvCPPDeA4Qrwq1yvrzqZiUiIpKPnlz5JHcuupO4xykKFjF1/FRK+pS06DmVgxYREalDuCLMnYvuJOpR4sSJxCKENrX8mB0K0CIiInUIbQoR81j1csAClPYtbfHzKkCLiIjUobRvKV2CXQgQoMAKuHXMrS1evA2qgxYREakWrggT2hSitG9pdRAu6VPC1PFTa61vaQrQIiIiJILzNXOvIRKL1GoIVtKnpNUCc4qKuEVEREjUNUdiEeLEqYpXtUpDsLooBy0iIp1aqli7e1F3ioJFVMWrKAwUtkpDsLooQIuISKeVWax90xduYkdkR6vWNeeiAC0iIp1WZrH2jsgOpoyY0tbJAhpZB21mPcxsppl9YGbvm9kpZnaomb1kZh8mH3umbf8jM1tlZivNbELzJ19ERGT/lfYtpShYRNCCeVGsna6xOejfAX9z94vNrAg4ELgVeNnd7zKzW4BbgJvN7ARgMjAMOAKYZ2bHuqf19hYREWlDbdWFqiEaHKDN7BDgy8CVAO4eASJmNhE4PbnZdGA+cDMwEXjc3SuBj81sFTAaeKuZ0i4iItJkbdGFqiEaU8Q9BNgMPGRmb5vZA2Z2ENDX3TcCJB/7JLc/Eliftn95cl0NZnatmYXMLLR58+b9uggREZG6tPZMVM2hMUXcBcDJwA/cfZGZ/Y5EcXYulmWd11rhfj9wP0BpaWmt10VERJqirgFI8lljctDlQLm7L0ouzyQRsDeZWT+A5GNF2vYD0vbvD2xoWnJFREQaZ85Hc6iMVebNACQN1eAA7e7/ANab2dDkqjOBFcBs4IrkuiuAWcnns4HJZtbFzAYDxwCLmyXVIiIiDRCuCPPsqmfxZAFu0IJ51VK7Lo1txf0D4JFkC+7VwFUkgvwMM7saWAdMAnD35WY2g0QQjwLXqQW3iIi0ptCmENF4FADDmHj0xHZRvA2NDNDuHgay/fQ4M8f2dwB37Ee6REREmizVzzk1fOd5R53X1klqMI0kJiIiHVY+93OujwK0iIh0aPnaz7k+mm5SREQ6hPbY17kuykGLiEi71177OtdFOWgREWn3Mmelai99neuiAC0iIu1ePs9Ktb9UxC0iIu1ee26tnYsCtIiIdAjttbV2LiriFhGRdqejtdjORjloERFpVzpii+1slIMWEZF2pSO22M5GAVpERNqVjthiOxsVcYuISLvSEVtsZ6MALSIi7U5Ha7GdjYq4RURE8pACtIiISB5SgBYREclDCtAiIiJ5qFEB2szWmNl7ZhY2s1By3aFm9pKZfZh87Jm2/Y/MbJWZrTSzCc2deBER6Zg6w0hh9dmfVtxnuPuWtOVbgJfd/S4zuyW5fLOZnQBMBoYBRwDzzOxYd481OdUiItJhdZaRwurTHEXcE4HpyefTgfPT1j/u7pXu/jGwChjdDOcTEZEOrLOMFFafxgZoB+aaWZmZXZtc19fdNwIkH/sk1x8JrE/btzy5TkREJKfOMlJYfRpbxH2qu28wsz7AS2b2QR3bWpZ1XmujRKC/FmDgwIGNTI6IiHQ0nWWksPo0KkC7+4bkY4WZPUOiyHqTmfVz941m1g+oSG5eDgxI270/sCHLMe8H7gcoLS2tFcBFRKRzCFeEawTlzhqYUxpcxG1mB5lZt9RzYDywDJgNXJHc7ApgVvL5bGCymXUxs8HAMcDi5kq4iIh0HKmGYb9f+nuumXtNp269ndKYHHRf4BkzS+33qLv/zcyWADPM7GpgHTAJwN2Xm9kMYAUQBa5TC24REckmW8Owzp6DbnCAdvfVwMgs67cCZ+bY5w7gjv1OnYiIdAqphmFV8apO3TAsnWazEhGRNqeGYbUpQIuISF5Qw7CaNBa3iIi0Og3lWT/loEVEpFVpKM+GUQ5aRERa1ZyP5lAZq+z0Q3nWRwFaRERaTbgizLOrnsWTA0sGLagW2zkoQIuISKsJbQoRjUcBMIyJR09U8XYOCtAiItJq0ifC6BLswnlHndfWScpbaiQmIiKtRv2dG04BWkREWpX6OzeMirhFRETykAK0iIhIHlKAFhGRZqeRwppOddAiItKsNFJY81AOWkREmlW2uZ2l8ZSDFhGRZhGuCBPaFKJ7UXfN7dwMFKBFRKTJMou1b/rCTeyI7FBf5yZQgBYRkSbLLNbeEdnBlBFT2jpZ7ZrqoEVEpMnSh/BUsXbzaHQO2syCQAj4xN3PMbNDgSeAYmAN8E1335bc9kfA1UAM+Gd3f7GZ0i0iInlEQ3g2v/0p4v4X4H3gkOTyLcDL7n6Xmd2SXL7ZzE4AJgPDgCOAeWZ2rLvHmiHdIiLSxlKNwlIBWUN4Nq9GFXGbWX/gG8ADaasnAtOTz6cD56etf9zdK939Y2AVMLppyRURkXyQahT2+6W/55q512hAkhbQ2Droe4GbgHjaur7uvhEg+dgnuf5IYH3aduXJdSIi0s6pr3PLa3CANrNzgAp3L2voLlnWeZbjXmtmITMLbd68uaHJERGRNqRGYS2vMXXQpwLnmdnXga7AIWb2F2CTmfVz941m1g+oSG5fDgxI278/sCHzoO5+P3A/QGlpaa0ALiIi+UeNwlqeuTc+JprZ6cC/JVtx3wNsTWskdqi732Rmw4BHSdQ7HwG8DBxTVyOx0tJSD4VUTCIiIp2HmZW5e60iiOYYqOQuYIaZXQ2sAyYBuPtyM5sBrACiwHVqwS0iItIw+5WDbinKQYuISGeTKwetkcRERETykAK0iIhIHlKAFhGRGsIVYR547wENPtLGNJuViIhUe3Llk9y56E5iHqNLsAtTx09VF6o2ohy0iIgAiZzznYvuJOpRHCcSi2iEsDakHLSISCeXmvRiw64NxH3fSM4BC2iEsDakAC0i0kmFK8LM+WgOz656lmg8SkGggIJAAdF4lIAFuHXMrSrebkMK0CIinVCqrjnq0ep1MY9x0TEX0e/gfhq+Mw8oQIuIdDLpdc0phlEYKOTco85VYM4TCtAiIp1IuCLMH9/5I7G0kZeDBLnw2As576jzFJzziAK0iEgnkd6FynEMI2hBbh1zK5OGTmrr5EkGBWgRkU4gs1jbME454hS+P/L7yjXnKQVoEZEO7MmVTzJv3Ty6BrvW6EIVtKCCc55TgBYR6aCeXPkkP1/48+rlgkABxFEXqnZCAVpEpIOat25ejeXjDz2erwz8irpQtRMK0CIiHdRZA8/izQ1vVi9fcPQFagzWjihAi4h0UKlgPG/dPM4aeJaCczujAC0i0oFNGjpJgbmd0mxWIiIieajBAdrMuprZYjN7x8yWm9nPkusPNbOXzOzD5GPPtH1+ZGarzGylmU1oiQsQERHpiBqTg64EvuLuI4ES4GwzGwvcArzs7scALyeXMbMTgMnAMOBs4A9mFmzOxIuIdFbhijAPvPcA4YpwWydFWkiD66Dd3YFdycXC5J8DE4HTk+unA/OBm5PrH3f3SuBjM1sFjAbeao6Ei4h0VuGKMNfMvYZILEJRsIip46eq21QH1Kg6aDMLmlkYqABecvdFQF933wiQfOyT3PxIYH3a7uXJdSIisp9Sk11EYhHixKmKVxHaFGrrZEkLaFQrbnePASVm1gN4xsyG17G5ZTtErY3MrgWuBRg4cGBjkiMi0qmkcs6VsUocJ0CAwkAhpX1L2zpp0gL2qxW3u28nUZR9NrDJzPoBJB8rkpuVAwPSdusPbMhyrPvdvdTdSw877LD9SY6ISKcw+6PZNYLz2CPGqni7A2tMK+7DkjlnzOwA4CzgA2A2cEVysyuAWcnns4HJZtbFzAYDxwCLmyvhIiKdSbgizKxVs/BkQWRBoECTXXRwjSni7gdMT7bEDgAz3P05M3sLmGFmVwPrgEkA7r7czGYAK4AocF2yiFxERBoptClENL5vqsjzjz5fwbmDa0wr7neBk7Ks3wqcmWOfO4A79jt1IiKdWLgiTGhTiNK+pZT2LaUoWERVvIrCQCHnHnVuWydPWpiG+hQRyUPZulJNHT+1OmAr99zxKUCLiOSh0KZQra5UU0ZMUWDuRDQWt4hIHkoVaQctqK5UnZRy0CIieaikT4mKtDs5BWgRkTxV0qdEgbkTUxG3iIhIHlKAFhERyUMK0CIiInlIAVpERCQPKUCLiIjkIQVoERGRPKQALSIikocUoEVEWkG4IswD7z1AuCLc1kmRdkIDlYiItKBwRZjZH81m1qpZROPR6okvNACJ1EcBWkSkhaRmpKqMVeI4QPXEFwrQUh8FaBGRFhCuCHP3krvZG9tbvc4wTXwhDaYALSLSzMIVYa5+8Woi8Uj1uoJAARcefSHnHnWucs/SIArQIiLNLLQpRFW8qsa6C4++kJ+c8pM2SpG0R2rFLSLSzEr7llIYKKxeLgoUce5R57ZhiqQ9anAO2swGAA8DhwNx4H53/52ZHQo8ARQDa4Bvuvu25D4/Aq4GYsA/u/uLzZp6EZE8VNKnhAcnPMicj+bgOOcddZ6KtaXRzN0btqFZP6Cfuy81s25AGXA+cCXwqbvfZWa3AD3d/WYzOwF4DBgNHAHMA45191iuc5SWlnooFGrSBYmItLZwRZjQphClfUsViKXRzKzM3Wu1HGxwDtrdNwIbk893mtn7wJHAROD05GbTgfnAzcn1j7t7JfCxma0iEazf2v/LEBHJL6muVJFYRH2cpVntVx20mRUDJwGLgL7J4J0K4n2Smx0JrE/brTy5LvNY15pZyMxCmzdv3p/kiIi0mdCmEJFYhDjx6j7OIs2h0QHazA4GngL+1d0/q2vTLOtqlae7+/3uXurupYcddlhjkyMi0qZK+5ZSFCwiaEH1cZZm1ahuVmZWSCI4P+LuTydXbzKzfu6+MVlPXZFcXw4MSNu9P7ChqQkWEcknJX1KmDp+quqgpdk1phW3AQ8C77v7b9Nemg1cAdyVfJyVtv5RM/stiUZixwCLmyPRDZVquNG9qDs7Ijv0zyMiTZKrMVhJnxJ9t0iza0wO+lTg28B7ZpaajuVWEoF5hpldDawDJgG4+3IzmwGsAKLAdXW14G5umWPgBghQFCzipi/cpGAtIo2mxmDS2hrTivsNstcrA5yZY587gDv2I11Nlmq4kRqgPk6cSCzCnYvuJOYxghbk1jG3MmnopLZInoi0M9kagylAdy5la7excPVWxg7pxahBPVv8fB12qM9Uw43UP1SAAAELEPMYjhP1KHcuupNjeh6jfzIRySq9SDv1nVIVr1JjsE7o0UXruG3WMuLuFBUEeGTK2BYP0h02QJf0KeHfTvwtc1e/yQl9D6f7wVV0L+rOnYvuJOpRAOIe169gEcnqyZVPVpe4dQl2Yer4qWoM1kmVrd3GbbOWEY0nSmQj0TgLV29VgN5fZWu38dMndxGJDmdBwJhUOoALT+7PrWPgzkV3Evc4RcEiSvuWahQgEakhXBGu8WM+EosQ2hRiyogp+o7ohBau3kosvq+XcMCMsUN6tfh5O2yAXrh6K5FonLhDJOY8umgdTy0t57ZzxnBOn19QcODHTDxuHEB1w4+CQAETj56ocXNFOrnQphBxj1cvByygIu1OIls989ghvehSGCBSFScQMH4+cbjqoJti7JBeFBUEqKyK4yRGSIlUxblt1jJicSMYOIpjuhzK7gPnVtdTR+IRZv59JrNWzVKgFunE0tuwBCzArWNu1XdBB5YKyj0PLOLnzy0nEo3XqGceNagnj0wZ26oNxKARk2W0huaeLKNs7TaeWlrOzLJyYrE4ZkYs7tXDmRUEjDsmd+PX795Q3R0rXYEVqKW3SCelqq/OoWztNi5/YCGRaJxAWowIGtwwfijXnXF0i6ehyZNltEepXz4Xndy/+tdRekV/3J0tW/vxbyf+lqf+Pou/736FmEdJjUiqlt4inZcGH+kc0qtDwQkGDHensCDQKvXMdenQATolFahT0pvK9zywiJ8+uYvKqjOwA4ZQdMhSCnouwSxR/xT3OHM+msPsj2ZjGOceda7+aUVE2plcfZhT1aFV0TiFBQFuO2cY23ZHWrUoO5cOXcSdS/obtXD1Vn4zdyVpDfQo6LGIrn1ngcUpsAIsANF4ojVnUaCIByc8qCAtItJOpPowx+JOl8LafZhbewCSTJ2yiJv1i2HNAjigF+zZCsXjYMDoWjnqooJE67w4iaHSotvHsLvycAoOXE2scDsFPRdVb1sVr+KhZQ+xN7aXswaepfppkXZA9cmdV60+zFW1+zBnxoR80XED9PrFMP08iFYCcbAABLvA2XfVCtap1nk9Dyxi2YYdPLFkPbE9g4jsGUTwgLUEu5eBJXLQZsYr618B4M0Nb/LGJ29w1fCr9E8vkqc0hnbnkS0nvHD1VuJpJcWBQOv0YW4OHTdAr1kAsQiQ7Mvo8USwfv5GiMchEICv/wZKr6z162n4Ed0T9dRxxyqL2bPuGgq6LwWH4AGfEDygvHrbV9a/whufvKFuWSJ5SmNodw7prbHTu0il6phTrbRbqw9zc+i4Abp4HASLauagLQDxGOCJIP38jdD3hMT2aUXhlxWPY+g/nZLW8hsq/zEISNZPd90XoM0gEo/w5N+fZNaq2Tw44QH984vkgfTpZjWGdseTmVtOb41dlTYUZ1v1YW4OHTdADxgNV8yuWQd9QK9kDjpRXI3H4Z3HIPxYzUAeKGDUMeMZdfBhcPRlMHF4dY46vmMMe4HCHksIdtmAJ3PoZlAZizDrgwUK0CJtLHMcbU0z27Fkyy1ntsZOL8bO1zrm+nTcAA2JID1gdO31z9+YCM7BLoDXLgqPReCD5xLLbz/CZV+7h7NGLGbTZ3t5u8fZ3Pb2GKLbxxA4YC2FhyylsEcZTgy8gOjuwa11dSKSxZMrn+SXC39JPPk/HYlF2BHZwZQRU9o4ZdJcsuWWrzvj6HabU86lYwfobEqvTBRrr1mQKAYHCD++LweNQfqIYrEq+OsN9PEYfYARFXMYdfzlVJa/w1N7TuKxTRdQ9dnJFB64mkDkaM7/yqmtfkkikpCa5CIVnEHjaLdXdXV9ypVbbq855Vw6X4CG2jnr9KLwf4Rh6V8gXpV4zQKJXHVKLMKwjx8C4KSCpXy370dsOfzLxHcX0fOEozmujg/HE+8uYO7qNxk/5ItccuK4lrgykU4ttClEzGPVywE0jnZ7lKvBV0p7rldujM4ZoDNlBuyRlyXqpnE4vAReuAlilckX9+WwDThq63yO2jo/EcjX/wkOn521WP2Jdxfwi7J/AYuyqOwJyj//Fz5nrUYnE2lGpX1L6RLsUmOSC41V0H6kcs2fbN+TtcFXuo6WW86mwQHazP4EnANUuPvw5LpDgSeAYmAN8E1335Z87UfA1UAM+Gd3f7FZU96SMgN23xP2Bewu3eH/7q29T6rues2CrAF67uo3waKYOU6U6X//bXUDs6c/fIaHzv6TgrRII2UOQFLSp4Sp46dqUJJ2KD3XXBAMUBBITFyRD2NiA4mxNd55FDAYeWn29k3NrDE56GnAfwMPp627BXjZ3e8ys1uSyzeb2QnAZGAYcAQwz8yOdU8re2pPMgN2z8Hw9sOw8Z1Ed63qgVCK9tVrZxg/5IssKXuMxO8VcI8nsuBANF6l1t8ijZRrABJNcpH/ytZu4+ml5Thw0cn9a3WTisXiTB49kCN6HJAfRdjrF8O0byQbFANvPwJXPtfiQbrBAdrdXzez4ozVE4HTk8+nA/OBm5PrH3f3SuBjM1sFjAbealpy80TplYm/HEOJZnNJzy4cW7GZUFGQbnHnnkO7k3yrKXTnkIp41v1EJDsNQNI+la3dxqVTEzllgJmh9Tx27Sm1Gn5dmAzcrS71vZ7+fb5mQaLBcEodpaXNqal10H3dfSOAu280sz7J9UcCC9O2K0+u61hydePKZs0CTtq7h5P2xHCMoZV7mdPtIAz4xs7d9O29M+euv3lrOvPWzuOsQWdx4ylXNE/aRdq50r6lGoAkj+Vqhb1w9VaqovsyJFUxb/tuUumZrb/dkgjAwaJEA+IBo5MDXxXuy0HXUVranFqqkZhlWZd12iwzuxa4FmDgwIEtlJw8kBrZLBbBAgWMrIpRsnUbAB4oIlAyPutuv3lrOtNW/hqAaSvDLFnxMKf2P50DjuinOjbp1FTfnH9SQbnngUX8/LnlWVthjx3Si8Lk0JsAhUFr225SqXkbYpHEiFOpasv0XPKA0XDlX/O6DjqbTWbWL5l77gdUJNeXAwPStusPbMh2AHe/H7gfEtNNNjE9+St9ZLPicQSguuGZjbws55u94P2/QMATHxx3lgc3sXzjE7ARAhgXHFRC/y4BSodMoGT4Za15RSJtTvXN+SO9kVfAEg28nNqtsEcN6slj14ytVQfdotYv3tfQN/P7NjVvg8eAQGKeBrfaueTGlJg2k6YG6NnAFcBdycdZaesfNbPfkmgkdgywuInnav8y3+AGvNlf/fxzPuoGpGZjSQZqzIi789Tnb2OfQ9HWMh4ABWkRaTXpxdjpjbzACQYM9+ytsFs1p7x+MUw7Z19X2bcfSeSGU9+/aaWbBItqzXjYlhrTzeoxEg3CeptZOfBTEoF5hpldDawDJgG4+3IzmwGsAKLAde22BXcb+8bQy+nz7h083e1gPuhSRMx9X11BMli7GVGcl1c8lz1AZ2v0ICLSBJmDidx2zrAajbxuO2cY23ZHWrc+OWcDr8i+bWJVNRt4ZZRu5tN3ZGNacV+a46Uzc2x/B3DH/iRK9ikefx0Ao9+fzdKuR/Dq5+/z+oG7iKXPb+pOgUMgMKb2AdYvJvzYhYSKApS+9Z+UXPp0m34AwxVh5nw0hy17ttD7gN4apEWkncocD3vb7kjrNvLKDMbpdcm1GngV7ctBBwtrN/Bqg+LrhtBIYu1A8fjrYPx1DAIuAB59/GZCW14kYl3Y9vkggsHP2b53JF/+1ndq7Rv+4CmuOaw7ETOK3Jn6wVOU5PogtnBOO1wR5uoXryYS3/dr9tlVz/LghAcVpEXySF3jYKdkGw+71YquswXj9LrkWg28nstdB53HFKDbocsm/4rL+BVQ/z/SC1EjYkbcjKrkctZQ2MCcdubITdnWpS8D1c9Dm0JUxatqHK+u/qvZziUizS/9e2TlP3Zy26xlxOJOl8La42CntPh42HVlGLIF48y65DZu4NUcFKDbufp+se62L+H+N4wYcS9gt30p63YNyWmHK8Jc8+J3icSrKAoUctMxk/lgwyKe3fURUY8l1o3+EXcvuZtILEJBIPHxisajFAWLuOkLN1EYKKyRg87VfzXXKFEi0jxSo3lV7Kzktb9vJhqLJ4bXdIglWnoRqco+DnZKs+eY6+uPnJItGOdxXfL+UoDu4M4//lRmln2PeNGqOqfDDHXtWiOnHeratVZOO7TymcTITWZEYhHu/OBhYiQ7uJtRFYswb+VTRGKVxPFkbjnZ1SI5J++DEx5sUB20RokSaX6poLx5ZyWvrKwgGqvZs7Uq5jUGrAgErGXGwc7W7akh/ZFTcgXjdppTzkUBuoMbNagnf/nWJfUWQ3Xrfhbus6pz2t26n1Vrm9K9eylypwowd+JmeLIlublT6M5ZlTGWxuNUGQTdwRKjjxd6jFI7kJLKCCVbt8KuLXCwQf9IrfOARokSaS6PLlrHC8s2MqzfIfzpzTXVA4RkMhKDhmBGNJboy/zzicP3P4ecq4g6V7enhvZHTulgwTgbBehOoCHFUFu29mP3umsJHLga3z2ELcX9am1TctxF3L/0L4SKgnSPwz2HHkIVTtCdibt2c96eKkoOO5JjPlpAqEshpXsT/4Chrl0orayiZM0SePtfanZ5yOyTmDqXRokSabJHF63j1mfeA2DBh1tyblcQgEu+MJALT+4PsP/1yg0pos7V7SmP+yO3FQVoARItMgteGUzVp4NyTu9WFj+Ge3bdxChfzss2jGuOCVK+7TW+2Ps4zj66R/Wv3JLw45RUfg6BIGCUfPZ54h8OrzngPNTuk5impDJCyfYd0CN7LltE6vbCso01llPjHAEUBI2vDO1D725dao3m1eDAnJ5LhoYVUefq9tQB65CbSgFagIa1yFy4eiuLo0ez0I8mAIQWGnEfzOMBY1LpAC4cmPwnT/8ng5rP336k5q/nbH0SIXefxkZ4cuWTzFs3j7MGnsWkoZMata9IPipbu42nlpZjwLAjutc7EMjXhverkXP+p3FD2FkZbZ4hNjP/R0smN3zIzFzdnjpBsXVjKEBLtfqKwtP7PVpyqNG4QyTmPLpoHU8tLeeRKWOBY1gYPZSx8UQuPPV81KCe+wac37UZDj4sd5/EXH0aG+jJlU/y84U/B+DNDW8CKEhLu5Pe/Qng0vvfIpLWsCtg1JqMIt1lYxITEL2wbCNfG96verlOdY1bnS7zfxRreBG1AnGDKEBLg6XnslOz1VRWxXGoHhT/6aXlPLW0nEg00WUj1eAk9SUCx7DwoOsZO6Ke+q26+jQ2wLx182otK0C3LfVrb5zMoTQvOrk/VRmtrlOjeNXVFeqyMQOzB+b0QHx4SSKYHtALXrgp97jV6TL/R0demvhTEXWzUYCWRknPZQ89vBtPLS1nZlk5sVhiNCGHfcP/xRKhOzN4V1bFCRiceXxf/um0o7J/sTSxPuqsgWdV55xTy9J21K+98TKH0nQSrawzc9C52oxUCyXbw/8AACAASURBVE2D92fB8ROh7wmJoLyrAv7+IqQPHGSBxF88bdqEOtqI1NnVSZqFArTst1Swvujk/jWK4Z5eWk5VNE4wmYNOD96pHHfMYe6KTcxfWcFj156SO0jv5z97KrdcVx20cnQtL3WPN+za0On7tTdk+Mx0mUNpXnRyfy46uX+j6qAJTYPn/iXx/KNXIFAA8Wj2bT2eKAqzQLIemdxtRFJUVN2iFKClyTLrrtMbmwE1ns9Ysp5ofF8OoCrmOYvnGvuFlmnS0Ek5i7WVo2t56fe4IFBAQaCAmMc6Tb/2VIOuLTsTxcXzV1YQjXuddcbpcjXcHBX4MJFrDfaCgq0QGAfkCJLvz6q5nCs4QyIwB7sk6o7/8Q7tbdzqjkgBWppdZsBOf/7zicP5SXKcX0gU2WXt0pVR/9aQL7TG0EhlLS/9Hsc8xkXHXES/g/t1yBKLzB+TZWu31WrQlVJfnXG6UYEPGVWwYF8QDk2D529MFkP7vqCaq5fD8RMTOeeU9Bx0oBCOnZBorJmqg1bdcV5RgJZWddmYgQw9vBtPLy2vs6tHZv1bQ7/QGqqhI5WpGHz/Zd7jjjC1aLZSnWw/Jheu3lqrQRckR+uqr844JbMb09l3JYNzWi7YcwyFmVJ6ZeIxsw5aueN2QQFaWl1DRjbLNpVdY6WGOMzWvaQhI5VlKwYHFLAbqL2PBpcatzr1QxLIWqqT7cfk2CG9+ELBh5zH6/RmB1vozvsM5sv9gwwaNZ7jGvJjM7Mb0/uzkoN/pLFA/b0cSq/cF6hBQbkdUYCWvNTUqeyyDXGYGaRjewZRueVgYodkD/6ZxeCzP5rNnI/mqN46TX0lDCV9StrlPSpbu41Lpy6sHrd6Zmg9F5cOyFqqM3ZIL0YXrGKUL6fMhjF2yBcZFfiQx4vuwOI1R8GzfwTgxYfh8AYMvJPZjen4ibD2LYhWJgYBOeV66HqIiqU7MAVoyVtNmcruhWUbmRx4ma8FF/NCbDQvLOtdI0A3pI67tG8pBYFCquJVBK0Aw2rVW0PnzVHne0O7zEE+GvNjb+HqrQyPfcCFBQsAeDo2DmNArUAMiXriR4vuTAbSWQQCp8CaBQQy5j5PqKdIOl22bkx9T1A/405EAVo6pOsPeYPRhQ8C8OXAeyw+5EhgTPXrC1dvZVjsA8YE3mdx7HgWrj6m1hd3bM8gdq+bQrxoFVWRozl25DAKArOqA3b3ou5c/eKU6jrWByc8QEmfkupcZfei7uyI7OiwwTsfGtrlqhO+77WPeOWDCuJxrzFDU0MbHJ558BqmFP2SIhL1vZOCr7HxoCIGZgZiRqcF43iiX3H1xA+FNYe1BaABRdLpMrsxqVtTp9LiAdrMzgZ+BwSBB9z9rpY+p8iYvW/glmiU48lluLH69TMPXsN3C++kkChVFLD24BHA0TWOsXD1Vk74fA9jdn/GYt/Dux/14Mh14+nR9R227x3Ja73WEolFwJxILMKsDxK5ratfnEIkHkme2SgKFHHZ0T9gxaZ/MH7IF7nkxMaNipav2npK0LK127jngYcZ5cu555Vh/PuU7wDwq6nTOY/XuSiwgy2B7jwdS9zvsXX8GMt03N53cIthyeUiYhRvmpcc2CMtEA8YnX3UuwGjM4a17QOHj1RLaWmUFg3QZhYE/gf4KlAOLDGz2e6+oiXPK8LxE7Fk9xJLLqc7bu87eCCKeZygxThu7ztAzdHGMoP4C1t28x+RqRRGolSxmF+v/iZdDooTBQrcOaQiziwWVAftRHx2IrFKpq38DeAsKXsM+H2tID3rtQcoW/Mio4onMPTAYrateIWeJ3yF476wL01N7RfeWA2pX66vEVh6mg+qKMt6Xfvr47df5aHAL5PvzzP89e0BAPw5+IvqnC/AN4PzcYwg8Zw/xmopHoelzbhkwcJ9dcCZw8/WNaKWArE0QUvnoEcDq9x9NYCZPQ5MBBSgpWVldi9Jb8UKyS/gLhCLJL6IsxQ5ZgbxM/0tColSYHHwKFezmq9v3MzbBxRx0p4Ih/feyYyKbnTxOMl8FgEHcyduhpsRIMqbK56qEaBnvfYAv1h9L1UGz69ewY1bd/BZMEDJS9P4gIc57gtnUbZ2Gz/7yz306PoOby45iCMPjdH7wF7sqtrGqOIJTDxtSo20f7BkHtvfepiivVuIdO3N6j79+L/PFmEYXx16efX26T8M0o8Rrghz9d++S5VXUWiFPHj2n7IG4K5rtzBsxWq6nlAMfWq+lp7DffrlbvxH8GGOJkrV6ql8wGNNDtKnBFfUeD9OCSa+VgqIYbZvuwJiGIbhOX+M1ZJrxqVcdcAKxtICWjpAHwmsT1suJ70iUKQlZXYvSdeQsb4zgnj3ky8mvmkJ8VgVgYJCjhh7CYe/sJSTdnwOwUICJeP5ZnguX/57Imh3j8XYEQxySCzO3b16VOe0zymqeZqyNS9SZRA3owq4u3cP4kBhD+dfV8zkuC+cxetLHmbDEbNZa4nA/x7gkUTpwPOr3weoDrAfLJnHzpe+xdxuXdl6UACA1/ccQLQoEbXeXH1v9bn3/TCoeYwXyh4nGo/gZsTiEV4oe5ySr9UM0B8smceg5y7NGXTTc7iOESBO0Bw8yrYVr0ATA/SRJeOJh39f/X4cWTIegPjbv8PTWk8HAgVgQYhHc/4Yyypb0FUgllbU0gHasqyr0XvfzK4FrgUYOLABU6GJNJf6vmyzBPFARg4qc/lIoF/495y0YxdGHCcAgQKOiW4j1KWA0kiMk866rMZpRhVP4PnV7xPFwZ1YMrcdBdb36QFAPL6oOojjjqc9RnHK1rxYHVxDK2by2369qLKMf7/kcmp7oPqYmccYULGdQvfqHxUDKrbXuj3bVrzC0Wk52Mygm57DjboRJ4C7U0UBPU/4SiPeqNzvT+DKObV+ZAWuyjKlKaj1s7Q7LR2gy4EBacv9gQ3pG7j7/cD9AKWlpbWH3hFpS/W1os2yXB00DuiFJRsFnQSclCNApIJi2ZoXOazfiUzb9AxRryIYKOJroyYDcOYJ5/DYknepwomTKDZ3IOBOgcOowROqj7e+Tw+im4wa5byeCP5Aje1TPwwyj1F6wsX8z0uzCR9QQMmeKN2/enGtW9PzhK9QtXoqeDRr0M3M4a4bcxsbN37SbHXQQONyuQrM0s6Ye8vFRDMrAP4OnAl8AiwBLnP35dm2Ly0t9VAo1GLpEWkPcjXOCi97lNDqF4nGDmbDZxUcXNgzax10uCLMVS9cSTQ1IxEQNGNEpBvdvahBddCQKMKur1FXvdusX6ycq0g9zKzM3Wt1g2jRAJ088deBe0l0s/qTu9+Ra1sFaJHmEa4IM+ejOWzZs4XeB/TuEONgi3RUuQJ0i/eDdvfngedb+jwisk97HWJTRPYJtHUCREREpDYFaBERkTykAC0iIpKHFKBFRETykAK0iIhIHmrxblaNYWabgbX7sWtvYEszJ6cj0H2pTfekNt2T2nRPatM9ya457ssgdz8sc2VeBej9ZWahbH3IOjvdl9p0T2rTPalN96Q23ZPsWvK+qIhbREQkDylAi4iI5KGOEqDvb+sE5Cndl9p0T2rTPalN96Q23ZPsWuy+dIg6aBERkY6mo+SgRUREOpR2HaDNbJKZLTezuJmVpq0vNrM9ZhZO/t3XlulsTbnuSfK1H5nZKjNbaWYTch2jIzOz283sk7TPxtfbOk1txczOTn4WVpnZLW2dnnxhZmvM7L3k56NTTq9nZn8yswozW5a27lAze8nMPkw+9mzLNLa2HPekRb9P2nWABpYBFwKvZ3ntI3cvSf59r5XT1Zay3hMzOwGYDAwDzgb+YGbB1k9eXvjPtM9Gp5xpLfne/w/wNeAE4NLkZ0QSzkh+Pjprt6JpJL4n0t0CvOzuxwAvJ5c7k2nUvifQgt8n7TpAu/v77r6yrdORT+q4JxOBx9290t0/BlYBo1s3dZJHRgOr3H21u0eAx0l8RkRw99eBTzNWTwSmJ59PB85v1US1sRz3pEW16wBdj8Fm9raZvWZm49o6MXngSGB92nJ5cl1ndL2ZvZsssupUxXRp9HnIzYG5ZlZmZte2dWLySF933wiQfOzTxunJFy32fZL3AdrM5pnZsix/df3a3wgMdPeTgBuAR83skNZJccvbz3tiWdZ1yCb89dyfPwJHASUkPie/adPEtp1O83nYD6e6+8kkiv+vM7Mvt3WCJG+16PdJQXMerCW4+1n7sU8lUJl8XmZmHwHHAh2iwcf+3BMSOaQBacv9gQ3Nk6L80tD7Y2ZTgedaODn5qtN8HhrL3TckHyvM7BkS1QHZ2rl0NpvMrJ+7bzSzfkBFWyeorbn7ptTzlvg+yfsc9P4ws8NSDaDMbAhwDLC6bVPV5mYDk82si5kNJnFPFrdxmlpd8osl5QISjeo6oyXAMWY22MyKSDQgnN3GaWpzZnaQmXVLPQfG03k/I5lmA1ckn18BzGrDtOSFlv4+yfscdF3M7ALg98BhwF/NLOzuE4AvAz83sygQA77n7q1aud9Wct0Td19uZjOAFUAUuM7dY22Z1jZyt5mVkCjOXQP8U9smp224e9TMrgdeBILAn9x9eRsnKx/0BZ4xM0h8Pz7q7n9r2yS1PjN7DDgd6G1m5cBPgbuAGWZ2NbAOmNR2KWx9Oe7J6S35faKRxERERPJQhyziFhERae8UoEVERPKQArSIiEgeUoAWERHJQwrQIiIieUgBWkREJA8pQIuIiOQhBWgREZE8pAAtksfMrNjMXjWz3Wb2gZnlHGfczHqY2fTkpPIVZnZ7I85zpZm5mf17xvpyMzt9/68g67m+aWZvJq9pfpbXS5IzSe1OPpZkvP5DM/uHme1IziDUpYHnPT15jU9nrB+ZXF8rLU2RfO/czHal/f2kOc8hHZsCtEh+ewx4G+gF/BiYaWaH5dj2P4EDgWISEzx828yuasS5PgVuboWZ3z4F7iUxdGQNyXHBZwF/AXqSmHd4VnI9ZjYBuAU4k8R1DgF+1ohzbwa+aGa90tZdAfy90VfRcD3c/eDk3y9a8DzSwShAS14zs5OT83rvNLMnzewJM/tl8rWeZvacmW02s23J5/3T9p1vZr9M5tZ2mdkcM+tlZo+Y2WdmtsTMitO2dzP7f2b2YfJ8vzCzo8zsreT2M9ICRZ3nbqZrPxY4Gfipu+9x96eA94CLcuxyLnC3u+929zXAg8B3G3HK94G3gB/uf6rr5+7z3H0G2WfPOp3EGNj3unulu/8Xiakxv5J8/QrgQXdf7u7bgF8AVzbi9BHgWRKTg5CcVOebwCP7cSkiLUoBWvJWMhg+A0wDDiWRm7wgbZMA8BAwCBgI7AH+O+Mwk4FvA0eSmLf1reQ+h5IISD/N2P5sYBQwFrgJuB+4nMTUjMOBSxtx7vRrec7Mtuf4yzVF3TBgtbvvTFv3TnJ9zlNlPB9ex7bZ/AT4oZkdWt+GZnZLHde0vZHnTRkGvOs1Jwl4l33XPIzEPUh5B+ibkSOuz8PAd5LPJwDLqWeqzbqu08xuqed8a5NVBQ+ZWe9GpFM6OQVoyWdjSeSm/svdq9z9adKmyHT3re7+VDLHuBO4Azgt4xgPuftH7r4DeAH4KJmDiwJPAidlbP8rd/8sObPTMmCuu69O2/+kRpy7mruf4+49cvydk2O3g4EdGet2AN1ybP834BYz62ZmR5PIPR+YK0050hkG5gI3N2Dbu+q4ph6NOW+a+q458/XU81z3JFu63wQONbOhJAL1ww3YJ+d1unutovqkLcAXSPyIG5VMo3Lq0mAK0JLPjgA+ychNrU89MbMDzex/zWytmX0GvA70SBZbpmxKe74ny/LBGeds0PYNPHdT7QIy64MPAXZm2Rbgn5Np/JBEPe5jQPl+nPc24Ptmdvh+7NtU9V1z5uup57nuSS5/Bq4HziBRStPs3H2Xu4fcPerum5LnG98KdfzSQShASz7bCBxpZunFtgPSnt8IDAXGuPshJOYBh5rFvC2lUec2sxcyWvOm/72Q4xzLgSFmlp47HJlcX4u7f+rul7v74e4+jMT/9+Js29bF3T8AngZurWs7M7u1jmva1djzJi0HTsx4z09k3zUvJ3EPUkYCm9x9ayPP82fg/wHPu/vu+jau6zrNrM77lCb1Q7M1Pp/SARS0dQJE6vAWEAOuN7M/At8g0Tp5fvL1biRyjNuTdaaZ9cktqVHndvevNfYE7v53MwsDPzWz/wC+RiJYZW0kZmZHAduTf+OBa0krdrdEN6L57n57A07/MxJ1vzmDibvfCdzZkGvJSGcQKCTx/RMws65AzN2rSLy3MeCfzew+4Jrkbq8kHx8GppnZIyR+wP0HiTYKqWNPS6btyrrS4O4fm9lpwOqGpNndM0ta6mVmY0i8Fx+SaJH+XyTuf2YRvkhWykFL3nL3CHAhcDWJL7pvAc8BlclN7gUOIFHXt5BEHWxraa1zTwZKgW0kuiVd7O6bAcxsXEZOdRSJVt47gf8PuDxZl54yAPi/hpzU3T8mkcs8qMlXUNu3Sfy4+SMwLvl8avK8EeB8EnXD20nUo5+fXI+7/w24G3gVWJv8S/9x1JhrfMPd62wc1kRDSHwudpJoz1DJvkaGIvWymtV7IvnNzBYB97n7Q22dlvYk2QXsSXc/pa3T0lKSrf7fAU5M5sZF2jUFaMlryWLIlSRyqpcD9wFD3H1jmyZMRKSFqQ5a8t1QYAaJ1tMfkSjiVXAWkQ5POWgREZE8pEZiIiIieUgBWkREJA/lVR107969vbi4uK2TISIi0mrKysq2uHutWeryKkAXFxcTCoXaOhkiIiKtxszWZluvIm4REZE8pAAtIiKShxSgRURE8lBe1UGLiEj9qqqqKC8vZ+/evW2dFGmErl270r9/fwoLCxu0vQK0iEg7U15eTrdu3SguLqbmzJySr9ydrVu3Ul5ezuDBgxu0T4cO0OGKMKFNIboXdWdHZAelfUsBCG0KUdq3lJI+JW2cQhGRxtu7d6+CcztjZvTq1YvNmzc3eJ8OG6DDFWGumXsNlbFKHCdAgIJA4nKj8ShFwSKmjp+qIC0i7ZKCc/vT2PeswzYSC20KEYlFcBJjjceJUxWvoipeVf08tEl9rkVEWtq9997L7t27q5e//vWvs3379iYfd/78+ZxzzjlNPk6+6rABurRvKUXBIgLJSwwQoDBQSGGgkKAFKQwUVhd5i4hIy8kM0M8//zw9evRowxS1Dx22iLukTwlTx09VHbSISAv47W9/y5/+9CcApkyZwvnnn8/ZZ5/NmDFjePvttzn22GN5+OGHeeCBB9iwYQNnnHEGvXv35tVXX60eNXLXrl2cffbZfOlLX2LhwoWMHDmSq666ip/+9KdUVFTwyCOPMHr0aBYvXsy//uu/smfPHg444AAeeughhg4d2sZ3oBW4e978jRo1ykVEpG4rVqxo9D6hNZ/6f7/yoYfWfNrk84dCIR8+fLjv2rXLd+7c6SeccIIvXbrUAX/jjTfc3f2qq67ye+65x93dBw0a5Js3b67eP7X88ccfezAY9HfffddjsZiffPLJftVVV3k8Hvdnn33WJ06c6O7uO3bs8KqqKnd3f+mll/zCCy90d/dXX33Vv/GNbzT5elpTtvcOCHmWmNhhc9AiIpJQtnYblz+wkEg0TlFBgEemjGXUoJ77fbw33niDCy64gIMOOgiACy+8kAULFjBgwABOPfVUAL71rW/xX//1X/zbv/1bnccaPHgwI0aMAGDYsGGceeaZmBkjRoxgzZo1AOzYsYMrrriCDz/8EDOjqqpqv9PennTYOmgREUlYuHorkWicuENVNM7C1VubdLxEpq+2zFbKDWm13KVLl+rngUCgejkQCBCNRgH4yU9+whlnnMGyZcuYM2dOpxmgRQFaRKSDGzukF0UFAYIGhQUBxg7p1aTjffnLX+bZZ59l9+7dfP755zzzzDOMGzeOdevW8dZbbwHw2GOP8aUvfQmAbt26sXPnzv0+344dOzjyyCMBmDZtWpPS3p4oQIuIdHCjBvXkkSljuWH80CYXbwOcfPLJXHnllYwePZoxY8YwZcoUevbsyfHHH8/06dM58cQT+fTTT/n+978PwLXXXsvXvvY1zjjjjP0630033cSPfvQjTj31VGKxWJPS3p5YrqKKtlBaWuqaD1pEpG7vv/8+xx9/fFsno4Y1a9ZwzjnnsGzZsrZOSl7L9t6ZWZm71+r3qxy0iIhIHlKAFhGRJisuLlbuuZkpQIuIiOQhBWgREZE8pAAtIiKShxSgRURE8pACtIiItLjTTz+dVDfa5phu8r777uPhhx8GEoOXbNiwoclpTJk/fz5vvvlm1nO1Jo3FLSIirer5559v8jG+973vVT+fNm0aw4cP54gjjmjw/tFolIKC7CFw/vz5HHzwwXzxi1+sda7WpBy0iIg0ypo1azjuuOO44oorOPHEE7n44our53t++eWXOemkkxgxYgTf/e53qaysrLV/cXExW7ZsAeDhhx/mxBNPZOTIkXz7299m586dDB48uHpCjM8++4zi4uJaE2Tcfvvt/PrXv2bmzJmEQiEuv/xySkpK2LNnD2VlZZx22mmMGjWKCRMmsHHjRiCRi7/11ls57bTT+N3vfsecOXMYM2YMJ510EmeddRabNm1izZo13Hffffznf/4nJSUlLFiwoPpcAOFwmLFjx3LiiSdywQUXsG3btupj33zzzYwePZpjjz2WBQsWNPk+N0uANrM/mVmFmS1LW3eomb1kZh8mH5s2tpyIiOy/9YthwW8Sj81g5cqVXHvttbz77rsccsgh/OEPf2Dv3r1ceeWVPPHEE7z33ntEo1H++Mc/5jzG8uXLueOOO3jllVd45513+N3vfke3bt04/fTT+etf/wrA448/zkUXXURhYWHWY1x88cWUlpbyyCOPEA6HKSgo4Ac/+AEzZ86krKyM7373u/z4xz+u3n779u289tpr3HjjjdXzUL/99ttMnjyZu+++m+LiYr73ve/xwx/+kHA4zLhx42qc7zvf+Q6/+tWvePfddxkxYgQ/+9nPql+LRqMsXryYe++9t8b6/dVcOehpwNkZ624BXnb3Y4CXk8siItLa1i+G6efBK3ckHpshSGdOLfnGG2+wcuVKBg8ezLHHHgvAFVdcweuvv57zGK+88goXX3wxvXv3BuDQQw8FYMqUKTz00EMAPPTQQ1x11VUNTtfKlStZtmwZX/3qVykpKeGXv/wl5eXl1a9fcskl1c/Ly8uZMGECI0aM4J577mH58uV1HnvHjh1s376d0047Lev1XXjhhQCMGjWqeqrMpmiWAO3urwOfZqyeCExPPp8OnN8c5xIRkUZaswBiEfBY4nFN04tfs00t2di5Hdw965SUp556KmvWrOG1114jFosxfPjwRh1z2LBhhMNhwuEw7733HnPnzq1+PTWHNcAPfvADrr/+et577z3+93//t8nTWKamygwGg9VTZTZFS9ZB93X3jQDJxz4teC4REcmleBwEi8CCicficfXvU49sU0sed9xxrFmzhlWrVgHw5z//uTq3mc2ZZ57JjBkz2Lo1MT/1p5/uy+d95zvf4dJLL21Q7jl9OsuhQ4eyefPm6rRVVVXlzBmnT2M5ffr06vW5psfs3r07PXv2rK5fru/6mqrNG4mZ2bVmFjKz0ObNm9s6OSIiHc+A0XDFbPjKjxOPA0Y3+ZDZppbs2rUrDz30EJMmTWLEiBEEAoE6W0APGzaMH//4x5x22mmMHDmSG264ofq1yy+/nG3btnHppZfWm5Yrr7yS733ve5SUlBCLxZg5cyY333wzI0eOpKSkpEaXqXS33347kyZNYty4cdXF7ADnnnsuzzzzTHUjsXTTp0/n3//93znxxBMJh8Pcdttt9aZvfzXbdJNmVgw85+7Dk8srgdPdfaOZ9QPmu/vQuo6h6SZFROrX1tNNtsbUkjNnzmTWrFn8+c9/brFztIXGTDfZkv2gZwNXAHclH2e14LlERKSD+MEPfsALL7zQLP2l27NmCdBm9hhwOtDbzMqBn5IIzDPM7GpgHTCpOc4lIiJtq6Wnlvz973/fYsduT5olQLt7rkqCM5vj+CIiIp1NmzcSExERkdoUoEVERPJQhw7QZWu38T+vrqJs7ba2ToqIiEijdNjZrMrWbuPyBxYSicYpCBiTSgdw4cn9GTVIQ4KLiLS2gw8+mF27djX5OPfddx8HHngg3/nOd5g2bRrjx49v1CxWdZk/fz5FRUXVs1iln6stdNgAvXD1ViLROHGHSMx5dNE6nlpazm3nDGPb7ghjh/RSsBYRaWc6wzSTKR22iHvskF4UFQRIjfLqQKQqzm2zlvHrF1dyyf++xaOL1rVlEkVE2qWbb76ZP/zhD9XLt99+O7/5zW/YtWsXZ555JieffDIjRoxg1qzaw1/Mnz+fc845p3r5+uuvZ9q0aQA5p4lM1xmmmUzpsAF61KCePDJlLJeOGUhRQYCgQSBgxOKOA9G4c9usZZSt3aa6ahHp8MIVYR547wHCFeEmH2vy5Mk88cQT1cszZsxg0qRJdO3alWeeeYalS5fy6quvcuONNzZ4Ao2qqqo6p4nM1JGnmUzpsEXcAKMCHzKq9wKu/XJXNm78hJ19xzJ1wWpKD3iRdQfs4ePPT+HppQN4b8Vj9Oj6Dm8uOYgjD40xpNfhVFR+hAHnDk/UPcxe9nD18ocbQszb+AZDD+rPIUXd6N61Fzv2bqV0yARKhl+WMz3hZY8SWv0ipUMm1DpmXfs9OfeGGuf7LLKT0GerOazgIL50+Jh6z51+3oamryHHAgitfrHW9WfbJv21OcseZktkJ72LutV97esXJ2bdKR7XLGMHi3RW4Yow18y9hkgsQlGwiKnjp1LSp2S/j3fSSSdRUVHBhg0b2Lx5Mz179mTgwIFUVVVx66238vrrrxMIBPjkk0/YtGkThx9+eL3HTJ8mEiAWi9GvX78Gp6m+/TOnmbzkkkvYuHEjkUiEwYMHXYfShAAAGwNJREFU13nsbNNMTpq0b+yt5p5mMqXjBujU/KfRSoqJU2wBWDeV3t2C/FPfXkTMKDx0Nv7pP9hwxBLWWKIY/D3Ad6ysPsxTS+7AzIimLceS06O9uevv1dsFgKKtZUyFrAEnvOxRrllyJxGDgk/LiDtEk+Xvzy65kwdz7Pfk3Bv4+Ya5tc5HAIjv4JUNczGgS45zp5+3oenLtV2Na9haBgZVJO5b6vpv2hDi7k/m1tgmmvbaXcnXANi7I/e1p96/WCQx+04zDfAv0hmFNoWIxCLEiVMVryK0KdSkAA2JHOzMmTP5xz/+weTJkwF45JFH2Lx5M2VlZRQWFlJcXFxrCseCggLi8Xj1cur11DSRqVmoGqu+/TOnmbzhhhs477zzmD9/Prfffvt+nTOluaeZTOmwRdzV85+S/CB4HGJVLO1SSMSMuBlRMzbZMqoMPBl0U4+YgRmxVHBOW65+Pe0xbkaVJXKL2YRWv0jEktuRDM7JY9a137yNb2Q9X2rfVJpzHaPGeRuavgYeqyrtfqXWzdv4Rq1t0l+rSrvuOq+9BeavFemsSvuWUhQsImhBCgOFlPatNS9Do02ePJnHH3+cmTNncvHFFwOJnGafPn0oLCzk1VdfZe3atbX2GzRoECtWrKCyspIdO3bw8ssvA42bJjKlo04zmdJxA3Rq/tPUJVoAgoWUVkYpcifoTqHD2f3HUegQSNaTWKq+xB2S2xWkLQfSX097/P/bu/8guev6juOv9+7dEUMDPQMkB8cFgzRIUg3c1YSRooyUiINNgWolVPuHEp2Rv+jUsXaMNB0dx/qr06EiUac67cUfEBqkYCItCkyblL1rWi+xSDxJcuSahOMKcdTc3e67f+x+N7vf++6P+/Hdn8/HzE129777/Xz2u3v7zufX+5PInS/o1g0bWL1JXa5suZI6/Ow5yz3vxp7rIssLnlup7KJyq61flefq1NnrFtThxp7rZh2TLPhdZ8HrLvvaY9i/FmhX6y9arx037dDdV9+94O7twNq1a3X69Gldcskl+a7kO++8U6lUKj82fOWVV8563qWXXqr3vOc9euMb36g777xTV199tSSpq6ur6m0iA626zWRg0babXAyLvt1kMIb5muXSrybyX/IHUl9WamZSA294T9GY6Uz6NzQ68b968eWkEue8KEl6+ZXrlXHpgt98SpI0+er16u89oeenDqgnuULndv6GLj7vInUkf8EYNGPQQE3Ue7tJzN9ctpts7QA9T0NHJrVreEwuad3F52vb7hHNZEpfp64kiVAA1A4Bunk1yn7QTat/VfesQLtt94gyGZdMCsfqwkQo//jBjQRpAMCCEaCrsGVDn9asXKZ9oxM6/atp3f/U6KxjXNL0TEb7RicI0ACABSNAV6mwVd23/Fw9PjKutT3n6dUzM3pwaEzpdEadHQltXL28zjUF0A7cXRas6kBTmOuQMgF6HrZs6NOWDX35+7df06t9oxMV83sPHZms6jgAKGfJkiWamJjQ8uXLCdJNwt01MTGhJUuWVP0cAvQiiBqzDivcXaurI6Ftt6zVweOvyJUN8ARsANXq7e3V2NiYTp06Ve+qYA6WLFmi3t7eqo8nQNdI0e5aMxl9YveI0rnZZg+mjmnn1msJ0gCq0tnZWTE9JZpf6yYqaTDB7lpJkxJm+eAsSdNp177RiTrWDgDQaGhB10iwu9a+0Ql1L+3Svd87qKmZbBrSzqQxuQwAUIQAXUOFY9VrVi7LJ0OpNAY9uP+oHh8Z183reoompwEAWhcBuk6qmVgmZYPzxx/+sSTp6edf0g+fO6nVF5yrg+OvErABoIURoBvc4yPjRff3HjqRv/308y/pP34+oStWLGPpFgC0GAJ0g7t5XY+efv6lkr//pwPHs/tBdyZIMwoALYRZ3A1uy4Y+ffrW39abes9XZ9IUlZKgMM0oAKA10IJuAkHmsiAT2fMnTuuR/zqe3VJZ2f9lkWYUAFoLAbqJFE4se9+1l+WXbE3+cqrsGDQpRgGg+RCgm1S1s8DDKUYZpwaA5sAYdIsrTDHKODUANA8CdIsrTDFaaZx66Mik7nvysIaOTNawhgCAKHRxt7jCFKOVxqnpCgeAxkGAbgPVjFdHdYUToAGgfmIP0Gb2gqTTktKSZtx9IO4yMXdBV/j0TIYlWwDQAGrVgr7B3Uunw0LdzaUrnCVbABA/uriRV6krPDxOve2WtRXXYAMA5qcWAdol7TUzl/QVd3+gBmUiBoXj1FPTGW3bPaKMuzoSpncPXKrbKmybCQCoXi0C9Fvc/biZXSTpB2b2P+7+VPBLM9sqaask9fWxdWIjKxynNjNl3LPBOu0a3H9UDw2P6R8/uFGS6AYHgAUyd69dYWb3SvqFu38u6vcDAwOeSqVqVh/MXTAG3b20S9sfPagz0xkFn6CkSe99c58eGh6jGxwAqmRmQ1ETqGNtQZvZuZIS7n46d/smSdvjLBPxKhynXrNymR4aHtODQ2NKp7Ozv12K7AZnbTUAzE3cXdwrJD1sZkFZg+7+/ZjLRI0Ewfr2a3rzXdqStGt4bFY3eGGa0YeGx/TS6TO6YNk5up1xawCIFGuAdvdRSW+KswzUX3j2d7BcK+gGD9ZWdy/t0h0P/Lum0meHVR5MHdPOrdcSpAEghGVWWHThbvCgdb1vdELT6eI5D9NpJ2sZAEQgQCNW4dZ1Z9KKWtCdSSNrGQBEIECjZvpXdWvn1murHoMmaxmAdkaARk1Vs3GHxO5aAMB+0GhIUbtrAUA7IUCjIQVZy5ImdtcC0Jbo4kZDqnZ3LQBoVQRoNKxqx6tLGdx/VI+PjOvmdT3asoE87wCaCwEaLWlw/1F9/OEfS5Kefj67FTlBGkAzYQwaLenxkfGy9wGg0RGg0ZJuXtdT9r6UXcp135OHNXRkslbVAoCq0cWNlhR0Z5cag2adNYBGR4BGy9qyoa/kuHPUOmsCNIBGQhc32lK166zpBgdQL7Sg0ZaqWWdNNziAeiJAo21VWmddKt0oyVMA1AIBGigh6AafnsmosyOh7qVdtKgB1AwBGigh3A1OixpALRGggTLC3eDVtKjZxxrAYiBAA1WqtkUdBO2OZEJ/2N+r26/pJVADmDMCNDAH5VrU4aA9NZPRzv1HtWt4TNtuWavJX07RqgZQNQI0ME+llmp1dSR0Zjojl+SSpqYz2rZ7RBl3JpcBqBqJSoAF6F/VrY/c8Pp8wA2C9pYNfepKmpImJRKmjPusrnAAKIcWNLDIgm7w267p1b7RCXUv7dL2Rw8WdYWHMbEMQBgBGohJ4Xj1mpXLSgbgcMYyxqsBSARooCbKZS0rmljGeDWAHMaggTor3Lij3Hg1G3cA7YUWNFBnhbPBS41XF3aDs74aaA8EaKABVBqvjlpf/eDQGIEaaGEEaKDBRI1XB93gReurCxKhMFYNtB7GoIEmEF5fbbnHXdFrqxmvBpofLWigSRSur941PKbvpo4pnfFZa6uD8eoz0xklTHr7G1boQ2+9nBY20GRiD9Bm9g5JfyMpKemr7v6ZuMsEWlk4EUp4vfS+0Yl8V3japb2HTuhfnzupG9ZcpIuWnaPbGLMGmkKsAdrMkpLuk/R7ksYkPWtmj7j7oTjLBdpBqbXVG1cvVzJhmsl4/rGZtOsHh05Ikr47NKaddzFmDTS6uFvQb5Z02N1HJcnMviVpsyQCNCo6cPKAUidSGlgxoPUXrZ/z87+79x49Mf6Mbuy5Tu++6Quzzz8yqNToHg2s3qT167ZE1yF0TKX7hc85f8lyvfLrCQ2s3iRJZcsqPE+pY8PnLTx/4XH9q7r10Y2jevrnOzWd/LU600uU/tUlSr7mRUnSS/93vfaNXqH+Vd3a/aOvauiFPeq/bJM2v/WDC74+83ldC1Hq+n9v5Jt6aeq0Luhapnete/+8yj4wMqhHRr4pk/Sude8v+X5XrserWt51nt7w2qsi36+5vL7gnC5Vfb7Cv4MrLh7Iv6YrX3uVfvLyofztws9q+HXPpX7VvN9R1zZ8vnLPL3fNS70/5V5TtX9/wbX//QrXZbHEHaAvkXSs4P6YpA0xl4kWcODkAd219y5NpafUlezSjpt2zClIf3fvPdp+fK8k6d+O75X23lMUpA+MDOquZz+tKZO6Joa0Q4r8oy085qPHU/rsi3tL3t+Re17wnIwkk9Q5MSSZNKPosgrL6ShxbOExwXk99+85oXMeGBnUVybu19T5QQlnpPNfyZfXed731Ju4XLt/9EP91eiXNG3SY6M/kaR8kJ7P9ZnP61qIqPIl6QO5xyRJv35Fu579lBJmcyr7wMigPpD6tKZy9//p2U/rYxHvd/i9KV2PV6XjY5HvV7Wv76PHU/pMrnxJ0vExJSq8pvDfQfL4XqV19vkquJ1Q9n3K5N6n4HV/rcS5y/19dLw8JLk0E/HZiLq2X8vdzn9eyjw/quzw3174/Ykqr7A+lco9MDJY9H7uTpW+Losp7lncFvGYFx1gttXMUmaWOnXqVMzVQbNInUhpKj2ljDKazkwrdSI1p+c/Mf5M9oZZ8f3g/KN7ssHOTNOWvT+rDqFjnhh/puz91OieoudIkud+N63SZYXLiTo26ryF5y88Z2p0j6Yt99qDH529P2PSqVee1tAL2eMyuceGXthTsk7VXJ/5vK6FiCo/6rXPmM257NToHk3r7HlKvd9zqYcU/X5V+/qeGH+m+JxVvKbw30FaivxMFJ5rJvS6S5277N+HlP9sRX4+I8ooOl+Z55e75qXen6jyIs9Votzw+zmtxfkMVxJ3gB6TdGnB/V5JxwsPcPcH3H3A3QcuvPDCmKuDZjGwYkBdyS4lLanORKcGVgzM6fk39lyXveFefD84/+pN6nIp6a5OV75rq9wxN/ZcV/b+wOpN+eckcuUmcr/rVOmywuVEHRs+r4XOX3jOgdWb1Om51x786Oz9ztwx/Zdlj0u6q8Ol/stmnyP4XTXXZz6vayGiyo967R251zyXsgdWb1Knzp6n1Ps9l3pI0e9Xta/vxp7ris9ZxfnCfwdJKfIzUXiujtDrLnXusn8fUv7zE/n5jCij6Hxlnl/umpd6f6LKizxXiXLD72fwNxQ38+CNiuPkZh2Sfirp7ZJelPSspC3ufjDq+IGBAU+l5tZSQutiDHp+Y9DBsYXjsFe+9ir9z8uHZo2flRqDHjoyqb/8h7/Wby75L/3fr9+kT/7xn0VOKmMMupp6MAbNGHR5Zjbk7rNaIbEG6FzB75T0JWX/8/Z1d/9UqWMJ0EBjuO/Jw/r83ueUcSlp0j03rdFHbnh9vasFtKRSATr2ddDu/pikx+IuB8DiCVKLhjftCBs6Mllyn2sAC0MmMQCzFO6wVSr4Fu6wxd7VwOIjQAOIVCoRSqBwh60gHzgBGlg8bJYBYF6CbvCkqWQ3OJt2APNHCxrAvFTqBg93gW+7Za0mfznFeDVQJQI0gHkr1w1e2AU+NZ3Rtt0jyrgzXg1UiS5uALEo7AJPJEwZ96LxagDl0YIGEIvCLvDupV3a/ujBisu2AJxFgAYQm8Iu8DUrl7FmGpgDAjSAmqi0bAtAMcagAQBoQARoAE2LddZoZXRxA2hKpBpFq6MFDaApRaUalWhVo3XQggbQlKJ23Apa1WemM0qY9PY3rNCH3no5LWs0JQI0gKYUlWr0vicP68x0Ri4p7dLeQyf0w+dO6t7fX0eaUTQdAjSAphVeurVx9XIlE6aZjOcfm047aUbRlBiDBtAy+ld1a/vmdUomLP9YwkSaUTQlWtAAWsqWDX1as3KZdg2PySWtu/h80oyiKRGgAbSccNc3aUbRjAjQAFrefNKMDh2ZzLfCb7+ml8COmiNAA0DI0JFJ3bEjmwRFkh5MHdPOrdcSpFFTBGgACNk3OqHpXHCWsjPBg8lldJWjVgjQABCycfVydXYk8i3ozqSpe2kXqUVRUwRoAAjpX9WtnXdtLBqDjkotGgTooSOTtKyx6AjQABAhamJZOLWoxKYdiA8BGgCqEJVaVIretKN/VbeGjkzqoeExmaTbmAWOeSBAA0CVolrVpTbtuOOBf9dUOpty9FvPHtNd171Oy17TSTc4qkaABoAFKLVpx3T6bD7wdMZ1/1OjMknndNINjuoQoAFggaI27ehMWr4FHXDNnmAGlMJmGQCwyPpXdWvn1mv1e1etUDJhCrbuSEjkA0fVaEEDQAz6V3Vrx/sH8kuwupd2sSc15oQADQAxmk8e8ADrq9sbARoAGhDrqxFbgDazeyXdJelU7qGPu/tjcZUHAK0kan21JNZWt5G4W9BfdPfPxVwGALSc8Prq7qVdRWurv506pj8auJRA3cLo4gaABhReX71vdKJobfVM2jW4/6geGh6j+7tFxR2g7zaz90tKSfpTd5+MuTwAaBnhCWbhtdWsq25tC1oHbWZPmNlIxM9mSV+WdLmk9ZLGJX2+xDm2mlnKzFKnTp2KOgQA2l6wtnrLhj7ddNUKdXUklLTK66oH9x/V+762X4P7j9awtlgM5u6Vj1poIWaXSXrU3deVO25gYMBTqVTs9QGAZlfNEqzB/Uf18Yd/nL//4etXkw+8AZnZkLsPhB+PcxZ3j7uP5+7eKmkkrrIAoN1Us7768ZHxovsPPD0qSUqYafvmddqyoS+2+mHh4hyD/qyZrVd2mOQFSR+KsSwAQMjN63r09PMv5e+7Z7+QM+76xO4R/fC5k7pg2Tm6nZngDSm2AO3u74vr3ACAyoIW8uMj41rbc56++szPNZPJDmumM669h05Ikr7z7FH90e/0sWSrwdRkDLpajEEDQHwG9x/Vtt0j+SBdKNgKc9sta3Xw+CtyiZZ1jZQagyZAA0AbGToyqV3DY/p26phmQtthJiRZwpTOBfCupGnn1msJ0jGr+SQxAEDjCSaX3XZNr3YNj+nk6TP60U9PKZ3OyMyKWtfTaWeNdR0RoAGgDRXOAi/cEvPe7x3U1ExGUjYxCntX1w8BGgDaXGGwXrNymXYNj1U1Bs12mPEiQAMA8qrdvzpqO0xJBOxFRIAGAMxZeDvMXcNjemh4LB+wt92yVpO/nCJYLwABGgAwZ+HtMF3KB+ypmYy27R5Rxj3fuiZIzx0BGgAwZ+HtMCVp1/CYpmeys8HTGWe3rQUiQAMA5iU8Xh0E7O6lXdr+6MF86zo8EzxYi00ylPII0ACARRGeDR41YWzoyKTu2LEvv5TrwdQxkqGUQIAGACy6UrPB941OaDoXnKXiZCgs2ypGgAYA1MzG1cvV2ZGYlQwlvGyLWeAEaABADfWv6tbOuzbOGoO+78nDzAIPIUADAGoqqvu7cNkWs8CzCNAAgLorXLZVaRa41B5pRgnQAICGUM0scCk6zWgrBmkCNACg4ZTLCR5OM9qqs8AJ0ACAphJOMxo1C7wVWtUEaABAUwmnGQ3PAg9PLGvWljUBGgDQdMJd4FGtaql4vDphpu2b12nLhr56VXtOCNAAgKYX1aqWiserM+7atntEa1Yua4qWNAEaANASSq2vTpgp4y5JymR81rrqRu0CJ0ADAFpW/6pubd+8LpuZLOPq6ixeV93IKUYJ0ACAlrZlQ1/JddWFXeCNlmKUAA0AaHml1lVXm2K0Ht3gBGgAQNuqJsVo0A1+ZjqjZKJ2M8EJ0ACAtlYpxei+0Qmdmc7IJc1kajcTnAANAEBOqZngyYRpJpObCe6zZ4LHIRHr2QEAaHLBTPCOhClhUleJHbYWGy1oAAAqKDcTPC4EaAAAqlBuh604LKiL28zebWYHzSxjZgOh3/25mR02s+fMbNPCqgkAQHtZaAt6RNJtkr5S+KCZXSXpvZLWSrpY0hNm9lvunl5geQAAtIUFtaDd/Sfu/lzErzZL+pa7n3H3n0s6LOnNCykLAIB2Etcs7kskHSu4P5Z7DAAAVKFiF7eZPSFpZcSv/sLdd5d6WsRjXuL8WyVtlaS+vubYoxMAgLhVDNDufuM8zjsm6dKC+72Sjpc4/wOSHpAkMztlZkfmUd4Fkl6ax/NaHddlNq7JbFyT2bgms3FNoi3GdVkV9WBcy6wekTRoZl9QdpLYFZL+o9KT3P3C+RRmZil3H6h8ZHvhuszGNZmNazIb12Q2rkm0OK/LQpdZ3WpmY5KulfTPZrZHktz9oKTvSDok6fuSPsIMbgAAqregFrS7Pyzp4RK/+5SkTy3k/AAAtKtWycX9QL0r0KC4LrNxTWbjmszGNZmNaxIttuti7pGTqwEAQB21SgsaAICW0tQBulQucDO7zMx+ZWYHcj/317OetUR+9PLM7F4ze7Hgs/HOetepXszsHbnPwmEz+1i969MozOwFM/tx7vORqnd96sHMvm5mJ81spOCx15rZD8zs+dy/tds1ogGUuCaxfp80dYDW2VzgT0X87mfuvj738+Ea16ueIq9JKD/6OyT9nZkla1+9hvDFgs/GY/WuTD3k3vv7JN0s6SpJd+Q+I8i6Iff5aNdlRX+v7PdEoY9J+hd3v0LSv+Tut5O/1+xrIsX4fdLUAbpMLvC2RX50VOnNkg67+6i7T0n6lrKfEUDu/pSkl0MPb5b0jdztb0j6g5pWqs5KXJNYNXWAruB1ZvafZvYjM/vdelemAZAf/ay7zey/c11WbdVNV4DPQ2kuaa+ZDeVSESNrhbuPS1Lu34vqXJ9GEdv3ScMHaDN7wsxGIn7K/W9/XFKfu18t6R5ls5qdV5sax2+e16Tq/OjNrsL1+bKkyyWtV/Zz8vm6VrZ+2ubzMA9vcfdrlO3+/4iZXV/vCqFhxfp9Eleqz0Uzn1zg7n5G0pnc7SEz+5mk35LUEhM+4s6P3uyqvT5mtkPSozFXp1G1zedhrtz9eO7fk2b2sLLDAVHzXNrNCTPrcfdxM+uRdLLeFao3dz8R3I7j+6ThW9DzYWYXBhOgzGy1srnAR+tbq7p7RNJ7zewcM3udqsyP3mpyXyyBW5WdVNeOnpV0hZm9zsy6lJ1A+Eid61R3ZnaumS0Lbku6Se37GQl7RNKf5G7/iaRSuxm2jbi/Txq+BV2Omd0q6W8lXahsLvAD7r5J0vWStpvZjKS0pA+7e00H9+ul1DVx94NmFuRHn1H75kf/rJmtV7Y79wVJH6pvderD3WfM7G5JeyQlJX09l0O/3a2Q9LCZSdnvx0F3/359q1R7ZrZT0tskXZDbb+GTkj4j6Ttm9gFJRyW9u341rL0S1+RtcX6fkEkMAIAG1JJd3AAANDsCNAAADYgADQBAAyJAAwDQgAjQAAA0IAI0AAANiAANAEADIkADANCA/h8HFL31G5CnygAAAABJRU5ErkJggg==","text/plain":["<Figure size 576x576 with 2 Axes>"]},"metadata":{"needs_background":"light"},"output_type":"display_data"}],"source":["fig, axs = plt.subplots(2, 1,figsize = (8,8))\n","axs[0].plot(X.flatten(),val_Opt,'.',label = 'optimal')\n","axs[0].plot(X.flatten(),v_policy[-1] - v_policy[-1][0],'.',label = 'policy iteration')\n","axs[0].plot(X.flatten(),v_value[-1] - v_value[-1][0],'.',label = 'value iteration')\n","\n","axs[1].plot(X.flatten(),pi_Opt,'.',label = 'optimal')\n","axs[1].plot(X.flatten(),U[np.argmax(PI_policy[-1],axis=1)],'.',label = 'policy iteration')\n","axs[1].plot(X.flatten(),U[np.argmax(PI_value[-1],axis=1)],'.',label = 'value iteration')\n","\n","plt.legend()\n","plt.title('gamma = 0.9, N = 100, M = 5')\n","plt.show()"]},{"cell_type":"markdown","metadata":{"id":"FXYHuKHmNI3R"},"source":["# paper / project presentation (4 pts)\n","\n","You will give a combined paper / project presentation during the final two weeks of class meetings.\n","\n","The paper presentation constitutes 30% of the course grade.  \n","\n","The project constitutes 50% of the course grade, with the project presentation constituting 1/5 of that portion, so 10% of the course grade.\n","\n","The presentations will be evaluated according to a simple rubric:  ***what did we learn, and how do we know we learned it?***  Your goal should be to (i) teach us something new (i.e. not covered in the course) and (ii) provide simple examples or exercises we can use to assess our understanding.  You should endeavor to adopt the notation and terminology from lecture (this may require translation from source material), and all new concepts and techniques must be defined and explained.\n","\n","Each presenter will be allotted 20min of lecture time, including answering questions that arise during or after the presentations -- ***there will be a hard stop after 20min***.  To satisfy this constraint while maximizing your expected grade, you should prepare < 10min of material for the paper portion of the presentation and < 5min of material for the project presentation.  (*Hint:* you probably need to prepare less material than you expect -- to calibrate, I encourage you to practice your presentations with a colleague who actively asks questions.)\n","\n","a. *Update your availability in the shared document accessible via Canvas -> Collaborations.*\n","\n","b. Specify whether you are planning a solo or group project -- if group, specify your group members."]},{"cell_type":"markdown","metadata":{},"source":["a. done\n","b. solo"]},{"cell_type":"markdown","metadata":{},"source":[]}],"metadata":{"anaconda-cloud":{},"colab":{"collapsed_sections":[],"name":"hw3.ipynb","provenance":[]},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.8.3"}},"nbformat":4,"nbformat_minor":0}