supply_chain.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "import scipy.sparse as spa\n",
    "from cvxpylayers.torch import CvxpyLayer\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import latexify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "# generate problem data\n",
    "n = 4   # nodes\n",
    "k = 2   # suppliers (with prices p)\n",
    "c = 2   # retail (with demand d)\n",
    "m = 8   # links\n",
    "\n",
    "supply_links = [0, 1]\n",
    "retail_links = [6, 7]\n",
    "internode_links = [2, 3, 4, 5]\n",
    "\n",
    "# Incidence matrices (nodes x links)\n",
    "A_in = np.array([[1, 0, 0, 0, 0, 0, 0, 0],   # 1 (supply)\n",
    "                 [0, 1, 0, 0, 0, 0, 0, 0],   # 2 (supply)\n",
    "                 [0, 0, 1, 0, 0, 1, 0, 0],   # 3 (retail)\n",
    "                 [0, 0, 0, 1, 1, 0, 0, 0],   # 4 (retail)\n",
    "                 ])\n",
    "\n",
    "A_out = np.array([[0, 0, 1, 1, 0, 0, 0, 0],   # 1 (supply)\n",
    "                  [0, 0, 0, 0, 1, 0, 0, 0],   # 2 (supply)\n",
    "                  [0, 0, 0, 0, 0, 0, 1, 0],   # 3 (retail)\n",
    "                  [0, 0, 0, 0, 0, 1, 0, 1],   # 4 (retail)\n",
    "                  ])\n",
    "\n",
    "# Prices\n",
    "mu_p = torch.tensor([0, 0.1]).double()\n",
    "sigma_p = torch.tensor([0.2, 0.2]).double()\n",
    "mean_p = torch.exp(mu_p + sigma_p ** 2 /2).double().view(k, 1)\n",
    "var_p = (torch.exp(sigma_p ** 2) - 1) * torch.exp(2 * mean_p + sigma_p ** 2)\n",
    "\n",
    "# Demands\n",
    "mu_d = torch.tensor([0.0, 0.4]).double()\n",
    "sigma_d = torch.tensor([0.2, 0.2]).double()\n",
    "mean_d = torch.exp(mu_d + sigma_d ** 2 /2).double().view(c, 1)\n",
    "var_d = (torch.exp(sigma_d ** 2) - 1) * torch.exp(2 * mean_d + sigma_d ** 2)\n",
    "\n",
    "# Uncertainty distribution (prices and demands)\n",
    "w_dist = torch.distributions.log_normal.LogNormal(torch.cat([mu_p, mu_d], 0), \n",
    "                                                  torch.cat([sigma_p, sigma_d], 0))\n",
    "\n",
    "# Capacities\n",
    "h_max = 3. # Maximum capacity in every node\n",
    "u_max = 2. # Link flow capacity\n",
    "\n",
    "# Storage cost parameters, W(x) = alpha'x + beta'x^2 + gamma\n",
    "alpha = 0.01\n",
    "beta = 0.01\n",
    "\n",
    "# Transportation cost parameters\n",
    "tau = 0.05 * np.ones((m - k - c, 1))\n",
    "tau_th = torch.tensor(tau, dtype=torch.double)\n",
    "r = 1.3 * np.ones((k, 1))\n",
    "r_th = torch.tensor(r, dtype=torch.double)\n",
    "\n",
    "\n",
    "# Plotting settings\n",
    "fig_width_pt = 469.75   # Get this from LaTeX using \\the\\textwidth\n",
    "inches_per_pt = 1.0/72.27  # Convert pt to inch\n",
    "fig_width = fig_width_pt * inches_per_pt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.13\n",
      "1.02\n",
      "1.52\n"
     ]
    }
   ],
   "source": [
    "print(np.round(np.exp(0.1 + 0.04/2), decimals=2))\n",
    "print(np.round(np.exp(0.0 +  0.04/2), decimals=2))\n",
    "print(np.round(np.exp(0.4 +  0.04/2), decimals=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([0.8140, 1.1559, 0.9484, 1.2243], dtype=torch.float64)\n",
      "tensor([[1.0202],\n",
      "        [1.1275]], dtype=torch.float64)\n",
      "tensor([[0.3268, 0.3268],\n",
      "        [0.4050, 0.4050]], dtype=torch.float64)\n",
      "tensor([[1.0202],\n",
      "        [1.5220]], dtype=torch.float64)\n",
      "tensor([[0.3268, 0.3268],\n",
      "        [0.8915, 0.8915]], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "print(w_dist.sample())\n",
    "print(mean_p)\n",
    "print(var_p)\n",
    "print(mean_d)\n",
    "print(var_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define linear dynamics\n",
    "# x = (h, p^{wh}, d) \n",
    "# u = u\n",
    "# w = (p^{wh}, d)\n",
    "# x_{t+1} = Ax_{t} + Bu_{t} + w\n",
    "A_d = np.bmat([[np.eye(n), np.zeros((n, k+c))],\n",
    "              [np.zeros((k+c, n)), np.zeros((k+c, k+c))]])\n",
    "A_d_th = torch.tensor(A_d, dtype=torch.double)\n",
    "B_d = np.vstack([A_in - A_out,\n",
    "                 np.zeros((k+c, m))])\n",
    "B_d_th = torch.tensor(B_d, dtype=torch.double)\n",
    "n_x, n_u = B_d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup policy\n",
    "# Parameters\n",
    "P_sqrt = cp.Parameter((n, n))\n",
    "q = cp.Parameter((n, 1))\n",
    "x = cp.Parameter((n_x, 1))\n",
    "h, p, d = x[:n], x[n:n+k], x[n+k:]\n",
    "\n",
    "# Variables\n",
    "u = cp.Variable((n_u, 1))\n",
    "h_next = cp.Variable((n, 1))\n",
    "\n",
    "# Cvxpy Layer\n",
    "stage_cost = cp.vstack([p, tau, -r]).T @ u\n",
    "next_stage_cost = cp.sum_squares(P_sqrt @ h_next) + q.T @ h_next\n",
    "constraints = [h_next == h + (A_in - A_out) @ u, \n",
    "               h_next <= h_max,  \n",
    "               0 <= u, u <= u_max,\n",
    "               A_out @ u <= h, u[retail_links] <= d,\n",
    "              ]\n",
    "prob = cp.Problem(cp.Minimize(stage_cost + next_stage_cost), constraints)\n",
    "policy = CvxpyLayer(prob, [x, P_sqrt, q], [u])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stage_cost(x, u):\n",
    "    assert x.shape[0] == u.shape[0]\n",
    "    batch_size = x.shape[0]\n",
    "    r_batch = r_th.repeat(batch_size, 1, 1)\n",
    "    tau_batch = tau_th.repeat(batch_size, 1, 1)\n",
    "    \n",
    "    h, p, dh = x[:,:n], x[:, n:n+k], x[:, n+k:]\n",
    "    \n",
    "    m = len(u)    \n",
    "    \n",
    "    # Selling + buying + shipping cost\n",
    "    s_vec = torch.cat([p, tau_batch, -r_batch], 1).double()\n",
    "    S = torch.bmm(s_vec.transpose(1, 2), u)\n",
    "    H = alpha * h + beta * (h ** 2)  # Storage cost\n",
    "        \n",
    "    return torch.sum(S, 1) + torch.sum(H, 1)\n",
    "\n",
    "def simulate(x, u):\n",
    "    assert x.shape[0] == u.shape[0]\n",
    "    batch_size = x.shape[0]\n",
    "\n",
    "    A_batch = A_d_th.repeat(batch_size, 1, 1)\n",
    "    B_batch = B_d_th.repeat(batch_size, 1, 1)\n",
    "    \n",
    "    zer = torch.zeros(batch_size, n, 1).double()\n",
    "    w = w_dist.sample((batch_size,)).double().view((batch_size, k + c, 1))\n",
    "    w_batch = torch.cat([zer, w], 1).double()\n",
    "    \n",
    "    return torch.bmm(A_batch, x) + torch.bmm(B_batch, u) + w_batch\n",
    "\n",
    "    \n",
    "def loss(policy, params, time_horizon, batch_size=1, seed=None):\n",
    "    P_sqrt, q = params\n",
    "    if seed is not None:\n",
    "        torch.manual_seed(seed)\n",
    "        \n",
    "    # Batchify input\n",
    "    x_b_0 = h_max * torch.rand(batch_size, n, 1).double()\n",
    "    w_0 = w_dist.sample((batch_size,)).double().view((batch_size, k+c, 1))\n",
    "    x_batch = torch.cat([x_b_0, w_0], 1).double()\n",
    "\n",
    "    # Repeat parameter values\n",
    "    P_sqrt_batch = P_sqrt.repeat(batch_size, 1, 1)\n",
    "    q_batch = q.repeat(batch_size, 1, 1)\n",
    "\n",
    "    cost = 0.0\n",
    "    x_t = x_batch\n",
    "    x_hist = [x_batch]\n",
    "    u_hist = []\n",
    "    for t in range(time_horizon):\n",
    "        u_t = policy(x_t, P_sqrt_batch, q_batch, solver_args={\"acceleration_lookback\": 0})[0]\n",
    "        x_t = simulate(x_t, u_t) \n",
    "        cost += stage_cost(x_t, u_t).mean() / time_horizon\n",
    "        x_hist.append(x_t)\n",
    "        u_hist.append(u_t)\n",
    "        \n",
    "    return cost, x_hist, u_hist\n",
    "\n",
    "\n",
    "def monte_carlo(policy, params, time_horizon, batch_size=1, trials=10, seed=None):\n",
    "    if seed is not None:\n",
    "        torch.manual_seed(seed)\n",
    "    results = []\n",
    "    x = []\n",
    "    u = []\n",
    "    \n",
    "    for i in range(trials):\n",
    "        cost, x_hist, u_hist = loss(policy, params, time_horizon, batch_size=batch_size, seed=seed)\n",
    "        results.append(cost.item())\n",
    "        x.append(x_hist)\n",
    "        u.append(u_hist)\n",
    "    return results, x, u\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(policy, params, time_horizon, lr, epochs, batch_size):\n",
    "    opt = torch.optim.SGD(params, lr=lr)\n",
    "    val_costs = []\n",
    "    best_params = []\n",
    "            \n",
    "    for epoch in range(epochs):\n",
    "        with torch.no_grad():\n",
    "            val_cost_mc, x_behav, u_behav = monte_carlo(policy, params, time_horizon, 1, trials=10, seed=0)\n",
    "            val_cost = np.mean(val_cost_mc)\n",
    "            val_costs.append(val_cost)\n",
    "\n",
    "        torch.manual_seed(epoch)\n",
    "        opt.zero_grad()\n",
    "        cost, _, _,  = loss(policy, params, time_horizon, batch_size, seed=epoch+1)\n",
    "        cost.backward()\n",
    "        print(\"epoch %d, valid %.4e\" % (epoch, val_cost))         \n",
    "        # TODO: Print gradient norm (possibly clip it)\n",
    "        # torch.nn.utils.clip_grad_norm_(params, 1)\n",
    "#         for p in params:\n",
    "#             print(p.grad.data.norm(2).item())\n",
    "        opt.step()\n",
    "#         scheduler.step(val_cost)        \n",
    "    return val_costs, [np.array(p.detach().numpy()) for p in params], x_behav, u_behav"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Baseline cost:  -0.2786436537604188\n",
      "Perform training\n",
      "epoch 0, valid -2.7864e-01\n",
      "epoch 1, valid -2.8086e-01\n",
      "epoch 2, valid -2.8249e-01\n",
      "epoch 3, valid -2.8470e-01\n",
      "epoch 4, valid -2.8614e-01\n",
      "epoch 5, valid -2.8761e-01\n",
      "epoch 6, valid -2.8670e-01\n",
      "epoch 7, valid -2.9183e-01\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sidereus/miniconda3/envs/python37/lib/python3.7/site-packages/diffcp/cone_program.py:259: UserWarning: Solved/Inaccurate.\n",
      "  warnings.warn(\"Solved/Inaccurate.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch 8, valid -2.9450e-01\n",
      "epoch 9, valid -2.9483e-01\n",
      "epoch 10, valid -2.9606e-01\n",
      "epoch 11, valid -3.0273e-01\n",
      "epoch 12, valid -3.0218e-01\n",
      "epoch 13, valid -3.0280e-01\n",
      "epoch 14, valid -3.0113e-01\n",
      "epoch 15, valid -3.0609e-01\n",
      "epoch 16, valid -3.0719e-01\n",
      "epoch 17, valid -3.0744e-01\n",
      "epoch 18, valid -3.0588e-01\n",
      "epoch 19, valid -3.0885e-01\n",
      "epoch 20, valid -3.1254e-01\n",
      "epoch 21, valid -3.1304e-01\n",
      "epoch 22, valid -3.1392e-01\n",
      "epoch 23, valid -3.1485e-01\n",
      "epoch 24, valid -3.1601e-01\n",
      "epoch 25, valid -3.1466e-01\n",
      "epoch 26, valid -3.1662e-01\n",
      "epoch 27, valid -3.1747e-01\n",
      "epoch 28, valid -3.1957e-01\n",
      "epoch 29, valid -3.2224e-01\n",
      "epoch 30, valid -3.2307e-01\n",
      "epoch 31, valid -3.2327e-01\n",
      "epoch 32, valid -3.2457e-01\n",
      "epoch 33, valid -3.2553e-01\n",
      "epoch 34, valid -3.2729e-01\n",
      "epoch 35, valid -3.2815e-01\n",
      "epoch 36, valid -3.2954e-01\n",
      "epoch 37, valid -3.2963e-01\n",
      "epoch 38, valid -3.3178e-01\n",
      "epoch 39, valid -3.3258e-01\n",
      "epoch 40, valid -3.3211e-01\n",
      "epoch 41, valid -3.3208e-01\n",
      "epoch 42, valid -3.3302e-01\n",
      "epoch 43, valid -3.3432e-01\n",
      "epoch 44, valid -3.3491e-01\n",
      "epoch 45, valid -3.3517e-01\n",
      "epoch 46, valid -3.3426e-01\n",
      "epoch 47, valid -3.3696e-01\n",
      "epoch 48, valid -3.3725e-01\n",
      "epoch 49, valid -3.3794e-01\n",
      "epoch 50, valid -3.3704e-01\n",
      "epoch 51, valid -3.3555e-01\n",
      "epoch 52, valid -3.3560e-01\n",
      "epoch 53, valid -3.3458e-01\n",
      "epoch 54, valid -3.3537e-01\n",
      "epoch 55, valid -3.3447e-01\n",
      "epoch 56, valid -3.3766e-01\n",
      "epoch 57, valid -3.3952e-01\n",
      "epoch 58, valid -3.4020e-01\n",
      "epoch 59, valid -3.3983e-01\n",
      "epoch 60, valid -3.4014e-01\n",
      "epoch 61, valid -3.3866e-01\n",
      "epoch 62, valid -3.3906e-01\n",
      "epoch 63, valid -3.3942e-01\n",
      "epoch 64, valid -3.4063e-01\n",
      "epoch 65, valid -3.4106e-01\n",
      "epoch 66, valid -3.4153e-01\n",
      "epoch 67, valid -3.4051e-01\n",
      "epoch 68, valid -3.4164e-01\n",
      "epoch 69, valid -3.4203e-01\n",
      "epoch 70, valid -3.4197e-01\n",
      "epoch 71, valid -3.4171e-01\n",
      "epoch 72, valid -3.4249e-01\n",
      "epoch 73, valid -3.4331e-01\n",
      "epoch 74, valid -3.4436e-01\n",
      "epoch 75, valid -3.4467e-01\n",
      "epoch 76, valid -3.4475e-01\n",
      "epoch 77, valid -3.4489e-01\n",
      "epoch 78, valid -3.4526e-01\n",
      "epoch 79, valid -3.4398e-01\n",
      "epoch 80, valid -3.4311e-01\n",
      "epoch 81, valid -3.4236e-01\n",
      "epoch 82, valid -3.4182e-01\n",
      "epoch 83, valid -3.4144e-01\n",
      "epoch 84, valid -3.4275e-01\n",
      "epoch 85, valid -3.4372e-01\n",
      "epoch 86, valid -3.4382e-01\n",
      "epoch 87, valid -3.4510e-01\n",
      "epoch 88, valid -3.4504e-01\n",
      "epoch 89, valid -3.4644e-01\n",
      "epoch 90, valid -3.4684e-01\n",
      "epoch 91, valid -3.4667e-01\n",
      "epoch 92, valid -3.4905e-01\n",
      "epoch 93, valid -3.4823e-01\n",
      "epoch 94, valid -3.4606e-01\n",
      "epoch 95, valid -3.4528e-01\n",
      "epoch 96, valid -3.4534e-01\n",
      "epoch 97, valid -3.4530e-01\n",
      "epoch 98, valid -3.4405e-01\n",
      "epoch 99, valid -3.4339e-01\n",
      "epoch 100, valid -3.4331e-01\n",
      "epoch 101, valid -3.4378e-01\n",
      "epoch 102, valid -3.4320e-01\n",
      "epoch 103, valid -3.4293e-01\n",
      "epoch 104, valid -3.4187e-01\n",
      "epoch 105, valid -3.4212e-01\n",
      "epoch 106, valid -3.4064e-01\n",
      "epoch 107, valid -3.4147e-01\n",
      "epoch 108, valid -3.4102e-01\n",
      "epoch 109, valid -3.4078e-01\n",
      "epoch 110, valid -3.4131e-01\n",
      "epoch 111, valid -3.4212e-01\n",
      "epoch 112, valid -3.4124e-01\n",
      "epoch 113, valid -3.4238e-01\n",
      "epoch 114, valid -3.4228e-01\n",
      "epoch 115, valid -3.4068e-01\n",
      "epoch 116, valid -3.4148e-01\n",
      "epoch 117, valid -3.4015e-01\n",
      "epoch 118, valid -3.3930e-01\n",
      "epoch 119, valid -3.3993e-01\n",
      "epoch 120, valid -3.3927e-01\n",
      "epoch 121, valid -3.3844e-01\n",
      "epoch 122, valid -3.4094e-01\n",
      "epoch 123, valid -3.4171e-01\n",
      "epoch 124, valid -3.4244e-01\n",
      "epoch 125, valid -3.4239e-01\n",
      "epoch 126, valid -3.4109e-01\n",
      "epoch 127, valid -3.4152e-01\n",
      "epoch 128, valid -3.4136e-01\n",
      "epoch 129, valid -3.4079e-01\n",
      "epoch 130, valid -3.4047e-01\n",
      "epoch 131, valid -3.4128e-01\n",
      "epoch 132, valid -3.4085e-01\n",
      "epoch 133, valid -3.4101e-01\n",
      "epoch 134, valid -3.4121e-01\n",
      "epoch 135, valid -3.4030e-01\n",
      "epoch 136, valid -3.4004e-01\n",
      "epoch 137, valid -3.3978e-01\n",
      "epoch 138, valid -3.4031e-01\n",
      "epoch 139, valid -3.4120e-01\n",
      "epoch 140, valid -3.4175e-01\n",
      "epoch 141, valid -3.4133e-01\n",
      "epoch 142, valid -3.4062e-01\n",
      "epoch 143, valid -3.4073e-01\n",
      "epoch 144, valid -3.4138e-01\n",
      "epoch 145, valid -3.4159e-01\n",
      "epoch 146, valid -3.4226e-01\n",
      "epoch 147, valid -3.4223e-01\n",
      "epoch 148, valid -3.4246e-01\n",
      "epoch 149, valid -3.4305e-01\n",
      "epoch 150, valid -3.4303e-01\n",
      "epoch 151, valid -3.4271e-01\n",
      "epoch 152, valid -3.4380e-01\n",
      "epoch 153, valid -3.4423e-01\n",
      "epoch 154, valid -3.4427e-01\n",
      "epoch 155, valid -3.4534e-01\n",
      "epoch 156, valid -3.4542e-01\n",
      "epoch 157, valid -3.4569e-01\n",
      "epoch 158, valid -3.4623e-01\n",
      "epoch 159, valid -3.4612e-01\n",
      "epoch 160, valid -3.4461e-01\n",
      "epoch 161, valid -3.4408e-01\n",
      "epoch 162, valid -3.4426e-01\n",
      "epoch 163, valid -3.4395e-01\n",
      "epoch 164, valid -3.4337e-01\n",
      "epoch 165, valid -3.4391e-01\n",
      "epoch 166, valid -3.4291e-01\n",
      "epoch 167, valid -3.4321e-01\n",
      "epoch 168, valid -3.4245e-01\n",
      "epoch 169, valid -3.4171e-01\n",
      "epoch 170, valid -3.4274e-01\n",
      "epoch 171, valid -3.4171e-01\n",
      "epoch 172, valid -3.4113e-01\n",
      "epoch 173, valid -3.4085e-01\n",
      "epoch 174, valid -3.3846e-01\n",
      "epoch 175, valid -3.3782e-01\n",
      "epoch 176, valid -3.3990e-01\n",
      "epoch 177, valid -3.4007e-01\n",
      "epoch 178, valid -3.3958e-01\n",
      "epoch 179, valid -3.3945e-01\n",
      "epoch 180, valid -3.3893e-01\n",
      "epoch 181, valid -3.3936e-01\n",
      "epoch 182, valid -3.3951e-01\n",
      "epoch 183, valid -3.3997e-01\n",
      "epoch 184, valid -3.3942e-01\n",
      "epoch 185, valid -3.4041e-01\n",
      "epoch 186, valid -3.3969e-01\n",
      "epoch 187, valid -3.3975e-01\n",
      "epoch 188, valid -3.4118e-01\n",
      "epoch 189, valid -3.4133e-01\n",
      "epoch 190, valid -3.4220e-01\n",
      "epoch 191, valid -3.4241e-01\n",
      "epoch 192, valid -3.4192e-01\n",
      "epoch 193, valid -3.4229e-01\n",
      "epoch 194, valid -3.4166e-01\n",
      "epoch 195, valid -3.4147e-01\n",
      "epoch 196, valid -3.4198e-01\n",
      "epoch 197, valid -3.4197e-01\n",
      "epoch 198, valid -3.4156e-01\n",
      "epoch 199, valid -3.4093e-01\n",
      "Final cost:  -0.3409255369772052\n",
      "Performance improvement: 22.35 % over baseline cost\n"
     ]
    }
   ],
   "source": [
    "# Perform training\n",
    "time_horizon = 20\n",
    "epochs = 200\n",
    "batch_size = 5\n",
    "lr = 0.05\n",
    "\n",
    "# Initialize value function V(x) = x'Px + q'x\n",
    "# centered at h_max/2 (between 0 and h_max) of each node\n",
    "torch.manual_seed(0)\n",
    "P_sqrt = torch.eye(n, requires_grad=True, dtype=torch.double)\n",
    "q = -h_max * torch.ones(n, 1, dtype=torch.double)\n",
    "# P_sqrt = torch.zeros(n, n).double()\n",
    "# q = torch.zeros((n, 1)).double()\n",
    "q.requires_grad_(True)\n",
    "params = [P_sqrt, q]\n",
    "\n",
    "# Baseline\n",
    "P_sqrt_baseline = torch.eye(n, dtype=torch.double)\n",
    "q_baseline = -h_max * torch.ones(n, 1, dtype=torch.double)\n",
    "# P_sqrt_baseline = torch.zeros(n, n).double()\n",
    "# q_baseline = torch.zeros((n, 1)).double()\n",
    "baseline_params = [P_sqrt_baseline, q_baseline]\n",
    "baseline_costs, x_behav_bl, u_behav_bl = monte_carlo(policy, baseline_params, time_horizon, batch_size=1, trials=10, seed=0)\n",
    "baseline_cost = np.mean(baseline_costs)\n",
    "print(\"Baseline cost: \", baseline_cost)\n",
    "\n",
    "print(\"Perform training\")\n",
    "val_cost, params, x_behav, u_behav = train(policy, params, time_horizon, lr, epochs, batch_size)\n",
    "print(\"Final cost: \", val_cost[-1])\n",
    "\n",
    "improvement = 100 * np.abs(baseline_cost - val_cost[-1])/np.abs(baseline_cost)\n",
    "print(\"Performance improvement: %.2f %% over baseline cost\" %  improvement)\n",
    "\n",
    "# Store final values\n",
    "P_sqrt_train = P_sqrt.detach().numpy()\n",
    "q_train = q.detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAACQCAYAAACMC/vEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1xUZf4H8M8zgQpKIt5CzRBLV61QLm5e8yVgXtbygtXupuUlrF1JLC9Z2/62y0vT3CzzBibuq7yLad5XAcmtzAsoIr6oFMEUUxJRIQxwvr8/OEMDzCDIDGcGP+/Xa17NnPOcc77z9Mh3nuec8xwlIiAiIqLyDHoHQERE5IiYIImIiCxggiQiIrKACZKIiMgCJkgiIiILmCCJiIgscNE7gLrUokUL8fHx0TsMIiJyEElJSb+ISEtL6+6qBOnj44OjR4/qHQYRETkIpVSWtXUcYiUiIrKACbIG8vLyEB4ejrS0NL1DISIiO9N9iFUpFQYgD4AngAwRSbZSxgtAAIBNIhKnLQ/RtgOAPNNyeykuLsYXX3yBkydP4uuvv4bBwN8XRET1la5/4ZVSngBCRSRORGIBzLNQxh+liTNaRCYD2GS2raeIxGrb+ts73pYtW2LhwoU4ePAgli1bZu/DERGRjvTuAj0N4IzZ5zwtIZrzAjDZ7HOuWZnZSilf7X1zSwdQSoUrpY4qpY7m5OTUOuDnnnsOoaGheP3113H69Ola74+IiByT3gnSE6XDqya5KE2IZbTepXmC9BKRZBHJAzALQJJSapOIzLJ0AK3nGSgigS1bWrySt0aUUoiOjkaDBg0wbNgw5Obm1nqfRETkePROkJZ4WluhlJoHINhskS9Kz0tCKbXPznGV8fHxwZdffonMzEyMGjUKt27dqqtDExFRHbHrRTpKqZlVrI7G7xfnmHgByLCyrzAAG0wX8Wifj4pIBoAxSql5SqkQe1+oY9K3b19ERUVh/PjxWLx4MaZOnVoXhyUiojqi9HxgsnahzTzTEKpSap+IhFoo54/Sq1QzzM45hgCI0xKk6YrWDNNnSwIDA8WWEwWICIYNG4YDBw7g5MmT4Cw9RETORSmVJCKBFtfpmSCBcrd5ACg956gtT0LpcKovgHiUnp8ESs9BNtPKzDTbNuN2vUdbJ0gAyMrKQrdu3RAQEIA9e/bAzc3NpvsnIiL7cegEWZfskSABYO3atXjuuecwePBgbNmyBQ0bNrT5MYiIyPaqSpCOeJGO0/nLX/6C5cuXY/fu3XjmmWdQXFyM3NxcxMfH6x0aERHdId1n0qkvwsPDUVRUhIiICAwZMgSpqam4fPkykpOT0aNHD73DIyKiGmIP0oamTJmCBQsWID4+Hm3atAEA7NmzR+eoiIjoTjBB2thrr72G9PR0HDlyBN27d8fevXv1DomIiO4AE6QddO7cGS4uLhg0aBC++eYb5Ofn6x0SERHVEBOkHQ0aNAjFxcVITEzUOxQiIqohJkg76tu3L9zc3PDf//5X71CIiKiGmCDtqGHDhhgwYAB27tzJ+VqJiJwME6SdjR8/HmfPnsWaNWv0DoWIiGqAM+nYmdFoRFBQEHJzc5Gens5ZdoiIHAhn0tGRwWDA3LlzkZmZicjISBw5cgR3048SIiJnxQRZB0JDQ/HMM89g+fLl6NmzJ0aMGIGCggK9wyIioiowQdYBpRTWr1+P7OxszJ8/Hzt27MCAAQNw48YNvUMjIiIrmCDrkLe3N2bMmIG1a9fi6NGj2L59u94hERGRFUyQOhg1ahTc3Nxw5MgRvUMhIiIrmCB14OrqCn9/fxw+fFjvUIiIyArdH3ellAoDkAfAE0CGiCRbKBMCIBdAIACISHR1t3VUPXv2xLJly1BcXAxXV1e9wyEiogp07UEqpTwBhIpInIjEAphnpcw8EUnWEuO86m7ryHr27ImbN2/i5MmTeodCREQW6D3E+jSAM2af85RS/uYFRCRPRAIAQCnlC+BodbfVtglXSh1VSh3NycmxbfS18Mc//hEALA6zXr9+va7DISKiCvROkJ4oHSI1yQXgZamgNsw6WXtVe1sRiRaRQBEJbNmypU2CtgUfHx+0aNEChw4dKrd8586daN68OdauXatTZEREBOifIC3xtLRQROIA7MPvCbLa2zoipRR69uyJuLg4jBo1CkOHDkV0dDSeffZZlJSUYOPGjXqHSER0V7PrRTpKqZlVrI7G7xfYmHgByLC2gYjEKaU2KaU21HRbR9S7d2/s2rULRUVFcHV1xe7du9GmTRuEhIQgLi4ORUVFaNCggd5hEhHdleyaIEVkflXrlVIbUf7iGs+KV6IqpcIBdBSRWdqiXO2/t93W0UVGRqJfv37o1asXAGDXrl3o1q0bTp06ha1bt+Lrr7/GwIEDdY6SiOjupPvTPMxu1QBQNpQKpVQSgGBtsWmmdX+tzPyqtrVGj6d53IkbN26gefPmmDp1Kj744AO9wyEiqreqepqH7gmyLjlLggSA4OBgXLp0qdJtIHl5ebj33nthMDji6WMiIufCx105oSFDhiAtLQ07duwoW3bw4EG0bdsWM2bM0DEyIqK7Q7V7kEqp7iJyXHvfAYAvABGRBDvGZ1PO1IP85ZdfEBISghMnTmDKlCkICgpCZGQkcnNz0ahRI2RmZqJ169Z6h0lE5NRs1YP0Nb0RkbMiEm++jGyrRYsWOHjwIMaPH4/Fixdj3LhxZVe6/vbbb/j444/1DpGIqF6rsgeplGqK0hlrJgNoiso35keJyBd2jdCGnKkHaS4/Px+nTp1Cx44d0bx5c4wZMwZ79+7FuXPn0LRpU73DIyJyWja5SEcpFaz1Gp2WsybIipKSkhAYGIiPPvoIU6dO1TscIiKnZZMh1orJUSnlU7uw6E4FBAQgICAAMTExuJuuQiYiqkvVTpBKqfeVUqOUUj2UUj8CCFNKTbJjbFSFCRMm4MSJEzh27JjeoRAR1Us1uUhng3a+MRzAfBFZAOCsfcKi2/nzn/+Mhg0bYtWqVXqHQkRUL9UkQXZQSt0LIBTABm0ZrxDRSbNmzTBq1Ch8/vnnWL58OXJzc2+/ERERVVtNEuQxAG8ACAOglFLvg7d56OrNN99E27Zt8fLLL8PX1xcrVqyA0WjUOywionqhJhfpnAUwF6VJcQyAOdowK+mkW7duOHnyJJKSktCjRw+Eh4dj4MCB+PHHH8uVKykp0SlCIiLnVZOLdIJR+vSMjgCaAVihlOKjJnSmlIK/vz8SEhKwYsUKHD9+HI888giioqJQXFyMcePGoU2bNsjMzNQ7VCIip1KT+yBHi8jm2y1zZPXlPsiqXLx4ERMmTMCePXvw0EMP4ccff0TDhg3h7++PAwcOwMXFrk84IyJyKraaau5qNZeRjry9vbFjxw7MmjULGRkZWLRoEVatWoWDBw9i2rRpKCoq0jtEIiKnUKMeJIAOAEwPJfYHkCcin9opNpu7G3qQ5goLC+Hm5gYAiIiIwOLFi9G5c2esWbMGAQEBOkdHRKQ/W82ksxlAPErnZh0DINmZkuPdyJQcAeCTTz7Bzp07UVBQgNGjR+PGjRs6RkZE5PhqcpFODwDhIvKSiLwM4KotpptTSoUppUK0//pbKROilPJXSoUrpcIrbBuulIpSSoXUNpb6bujQodiwYQPOnTvHZ0oSEd1GTR+YPN/0RkSOoZb3QSqlPAGEikiciMSi9CpZS2XmiUiyiESbymjJNENEokVkMoBNtYnlbtG7d2+89tpriIqKQmJiot7hEBE5rBrNpKPdC2nOs5bHfxrAGbPPeRV7kSKSJyIBAKCU8gVgOonohdLHcJnkWuqBaj3Mo0qpozk5ObUMt35455130L59e8ycOZOTnRMRWVGTBKmUUsu0CctHKaWWoTRJ1YYnKj9j0uI+tSHUydoLWq/TPEF6iUhyxe20HmagiAS2bNmyluHWD25ubvjXv/6FI0eOYMuWLXqHQ0TkkGp6kU40gJ7aK9pOF+lY7JWKSByAfSjfawQAKKXmAQi2Qyz11tixY9GlSxe88cYbnMeViMiCGt01rp13rPbzlZRSM6tYHY3S3qN5QvQCkFHF8eOUUpuUUhtMvUWlVBhKnzRSqfdI1rm4uGDBggUYPnw4OnXqhNdeew1+fn7o06cPmjblHPRERHadVkVE5le1Xim1EeUvzPGsmOi0q1Y7isgsbVGu2Tp/lN5ukqGdn4SIWE2wVN7QoUNx7NgxvPLKK3jjjTcAAO3bt8e+ffvQqVMnnaMjItJXtScKsFsApT3AsvOQ2lAqlFJJ+H3Y1HQTp79WZr6WHOPxe8L0EpFmVR3rbpsooCYuX76MpKQkPP/88wCALVu2oE+fPjpHRURkX1VNFKB7gqxLTJC398MPP2DIkCHIzMzEjBkz8O6778LV1VXvsIiI7MJWc7HSXaBTp044duwYJk6ciHnz5uHVV1/VOyQiIl3w0Q5Uyb333ovo6Gh4eHjgww8/REBAAF544QW9wyIiqlPsQZJV8+bNQ3BwMCZNmoQZM2YgPz9f75CIiOoMEyRZ5eLigs2bN2P8+PFYsGABAgICcOHCBb3DIiKqE0yQVKWmTZtixYoV2L9/Py5evIjHH38ce/fuRVpaGqepI6J6jQmSqmXAgAHYu3cvcnJy8MQTT+Dhhx/Gu+++q3dYRER2wwRJ1fbYY4/hxx9/xP79+zF69Gi8++67SElJ0TssIiK74H2QdEdyc3PRtWtXtGnTBgkJCfD09ITRaMT169dRXFyMFi1aQCmld5hERFXifZBkc15eXoiOjsbx48fx4IMPYty4cfD29kazZs3QqlUr9O/fH6dPn9Y7TCKiO8YESXfsySefRHJyMvz8/LB161YMGDAA//73v/Hee+8hNTUVfn5+2Lx5s95hEhHdEQ6xkk2ISLkh1QsXLmDMmDE4dOgQFi9ejJdeeqlsvYjg0KFD8PPzg5ub2x0dz2g0QinFYVwiqhUOsZLdVUxUbdu2RVxcHAYPHoy//e1v6NKlC+bMmYNt27Zh+PDh6NWrF6ZNm1atfefn5yM8PBzdu3fHl19+iWXLlsHLywsLFy60x1chIgLAHiTZWUlJCdatW4clS5bg0KFDAAB3d3c88sgjSEpKQnp6Ojp27FhpuytXrmDGjBnIy8tDSkoKzp49i/bt2yMrKwsA0LhxYzRu3BhZWVlo1KhRnX4nIqo/+DQPDROkvq5evYqTJ0/C19cXBoMBHTt2xOjRo/H555+XK3f9+nUEBwfjxIkT6NSpE9zc3DB//nz06dMH69atg7u7Ozw9PREaGoqYmBiMHz9ep29ERM6OCVLDBOlYZs2ahQ8++ADBwcHo378/XnrpJeTl5WHs2LFISkrCF198geHDh1vcVkTQvXt3iAhSUlJ4LpKI7gjPQZJDevPNNzFlyhRcuXIF//d//wcfHx88+uijSE9Px4YNG6wmR6D0nGdkZCRSU1OxcePGOoya7Olu+sHuSIqLi7FmzRoMHDgQW7Zs0Tsch6F7D1IpFQYgD4AngAwRSbZQJgRALoBAABCRaEv7EJG4qo7FHqTjSk9PxwcffAAAeO+99+Dt7X3bbYqKitCvXz+kp6cjOTnZ4rnMY8eO4b333kPnzp0xaNAg9OrVCw0bNrR5/FQ7ly5dwuzZsxEbG4vExET4+/vrHdJdw2g0YsCAAfjf//6HBg0awM3NDSdPnkS7du0qlc3JycHu3bvRuXNnBAUFwWBw/j6Www6xKqU8AcwTkcna530iEmqhTLyIBGifr4pIswrrN2n7YYK8y2RmZsLf3x/t2rXDypUr0b17dyQnJ6Nhw4YwGAwYOHAgiouLUVBQgFu3bsHd3R1jx47Fxx9/zESpk1OnTuHVV1/FhQsXMHnyZJw7dw7Lly/HzZs30ahRI3Tv3h1fffWVXYfNRQTZ2dlISUlBixYtEBQUdNcO08fExGDixIn46KOPMGTIEPTo0QN9+/bFqlWrUFxcjE2bNuHo0aPIzs7GwYMHUVJSAgBo0aIFOnTogE6dOuGpp55CQUEBYmNj0bdvX0yfPh1nzpzB9u3b4erqisDAQPTp00fnb2pZVQkSIqLbC0A4gJlmnzcB8K+ivC+AfRWWhQCYCSDkdscLCAgQqn/27Nkjnp6eAkDc3d0FQNnL29tbTp8+LdeuXZMvv/xSJk6cKACkX79+kpWVVW4/p06dkt27d+v0LRzXyZMnJT8/v9yywsJCuXnzZrW2z8rKkl9//VVERJYuXSouLi7i6ekpAQEBAkAMBoM8++yzkp6eLsuXLxcAEhsbW6MYr169Kjt37pRFixbJzz//bLXcvn37JDQ0VFq0aFGunXTu3Fk+++wzMRqNlbZJS0uThIQESU1NlZKSkhrF5eiuXr0qLVu2lD59+pR998WLF5erGwDSsWNH6du3r8yYMUOOHDkiq1evlgkTJsigQYOkZcuWZeXuu+8+ASAPPvigGAyGsuUGg0HWrFlj09hNbaq2ABwVaznH2oq6eGmJLdzsc5S1RKclwnkAfM2W+aN0aNZqgtSS8FEAR9u3b2+TCiXHc+3aNZk/f778/e9/l02bNsmaNWvkzTfflO+//75S2XXr1knDhg0FgAQFBcncuXNl/vz5ZcumTJkiRUVFdo/5t99+k5SUFPntt9/sfqw7FRMTIwCkXbt2smjRInnppZekc+fOYjAYxMPDQyZNmiSpqakiInLmzBlZtGiRFBQUlG0fFxcnrq6ucv/998uLL74oAORPf/qTXL58WYxGoyQnJ8vZs2fLyhcXF8vDDz8s3t7eEhsbWylhffvttzJ+/HiJioqSb775RtavXy9PP/20uLq6lv0x9vDwkLFjx4qPj4+4u7tL165dZdiwYfLkk08KAPHx8ZGJEyfKJ598Il999ZV8+umn4u/vLwAkJCREEhIS5ObNmxIfHy9PPfVUuUTx0EMPyeLFiyUpKanSjwZnUVhYKJ999pkMHz5cWrduLUopSU5OLltvNBolMTFRli5dKgsXLpQzZ85Uub+SkhL56quv5NtvvxWj0Shr164VPz8/mT59uly4cEF+/vlnGTBggBgMBomIiJCPP/5YLly4cNt9ZmVlSWpqqty6davcuuzsbJk9e7Y0a9ZMjhw5cucVoXG2BBlWRfkQlA6lliVIs/2wB0nVdubMGZkzZ44EBQWV/fEbNmyYTJ06VQBImzZt5IUXXpBp06bJ+PHjJSgoSB555BFZt26dxV5GdRUVFcnWrVtl5MiR0rhxYwEg/fv3l7y8PBt+uztnNBplyZIl8swzz0hERIQYDAZ5/PHHpXv37gJAGjduLMOHD5e33npLXnjhBWncuLEYDAYZMWKENGrUSABIjx49JCUlReLj48XDw0O6du1atv1f//rX2/74SEpKki5dupT9gImPj5fCwkLZtm2buLm5SYMGDcolLS8vL4mMjJT4+Hg5fvy4DBs2TJo0aSLDhw+XyMhIGTFihPj5+cl9990ns2bNstjzuHXrlixZsqRsJMLFxUUASNOmTeXtt9+WhIQEiYmJkcDAwLLjuru7y9SpUyUhIUHi4uIkLS3N4X7sfPfddxISEiKhoaFy7tw52bp1q7Rq1UoASIcOHWTs2LGydetWu8eRn58vw4cPL/sR6uHhIR9++KGkpaXJiRMn5J///KdERkbKpk2b5JVXXhEPD4+yen7ggQdkwoQJMmnSJOnZs6copcRgMMjo0aMlLS2t1rHpliC1xGXt5VnTIVatzFWt5xgOIEx7barYu7T0YoIkS3766Sc5cOBAWeLbtm2bhIWFSfPmzaVJkybi7e0tAwcOlEcffVQAyGOPPSYHDx6s0TFu3bola9euFV9fXwEgrVu3lpdfflnmzJkjLi4u0qNHD8nMzLTH16uW69evy4EDB2TEiBFlQ2VKKRk4cKAUFBRISUmJJCcnV0ouV65ckYiICHF1dZWRI0fKf/7zH7n33nvL/ri1bdtWfvrpJykpKZHDhw9X6g1YU1xcLDExMXL//feXS4b+/v5y6dIlSU1NlR07dsjx48erPdRbHQUFBbJ69WqZNm2abN26VW7cuFFuvdFolLS0NImNjZVx48bJPffcUy6+Bg0ayNNPPy3bt2+XxMREWbNmjURERMjKlStr9cOqOm7duiXbtm2ThIQESUpKkpEjRwoAadWqlTRp0qTsB5m/v7/s27ev2v8vbMloNEp6ero88cQT5erNYDCU/cBydXWV5557TqKiomTlypUyePBgad26tdx3333y2GOPydtvvy2nT5+2WUyO3IP0BBBl9nmfhTLhFXqNZyomUS05sgdJdlVSUiIxMTFl51lefPHFcsOJ1pw4cUJ69+4tAMTPz0+2bNkixcXFZet37dolHh4e4uHhIQsWLJBt27ZJdna2zeM/f/68TJ8+Xfz8/OSBBx6QXr16ydy5cyUsLKysx+Ti4iILFy4Uo9EohYWFd/RH/cyZMxIVFSUbN26Uixcv1irmwsJCiY6Oljlz5sjSpUvl2rVrtdqfrZ07d07i4uIkMTFRVq9eLREREdKsWbNKSROAPPvss1WeH60Jo9Eo33zzjYwdO1aGDBkiK1eulODg4HLH9fDwkHfeeUdu3LghP/zwg4SEhMj06dNt+oOiNvEnJSXJ6tWrZcWKFXLx4kW5efOmfPvtt3Zp+1Vx2ARZGhvCtKHTEPMkByBJS6CeZutnmvc4tXL+Wtko9iCpLty4cUNmzJghSinp2rWrLFmyRA4fPiznz5+XwsLCcmVXrVolLi4u0rx5c4mJibH6q/3s2bPSv3//sj9unp6e8vXXX9ss5qioKHF1dZV77rlHBg0aJGPHjpUePXqUDVFOmzZNduzYIZcuXbLZMe9Wv/76qyQmJsr+/fslJSVFiouLZc6cOWIwGEQpJQEBATJ48GCZPHmy5OTk3HZ/P/30k7z11lvy0EMPSadOnaRfv37i5eVVlgR9fHzKhnyXL18uu3fvlmXLllVr31R1gtT9Psi6xNs8yJb27duHF198sWx+WJMmTZqgc+fO6NSpE9atW4eQkBCsX78ezZs3r3J/IoLMzEycP38ekyZNwrlz5zB48GAopfCPf/yj2vcGHjt2DN999x0MBgMMBgNSU1PxySefYPDgwVi6dCk6dOhQVjY7OxteXl6cz7YOnDp1Cps3b8b+/fuRn5+PlJQUtGrVCu+//z68vb3h6+uLBx54oNztJj/88AP69u2LK1euIDg4GE2bNsXFixfxhz/8Af369cPo0aPh7u6Ow4cPo23btrj//vt1/IbOyWHvg6xrTJBkayKCrKwsHD9+HJcvX8Yvv/yCnJwcJCcn4/DhwxgzZgw+/fRTNGjQoEb7zcnJwYQJE5CZmYmff/4ZJSUl2Lt3L4KCgqxuc/nyZUydOhXr16+vtG7SpElYtmwZXFxcavwdyT6Sk5MxZswYZGRklC1r1qwZnnrqKYSFhaGoqAiRkZH49ddfkZiYiG7duukYbf3FBKlhgqS6ZDQabTLTyLlz5zBgwABcunQJQ4cORdeuXZGdnY1BgwZhzJgxAEoT9ZAhQ5CYmIjp06cjPDwcLi4uEBG4uLigdevWtY6DbK+wsBDp6em4du0avv/+exw8eBCbN29Gfn4+AKBp06ZISEjgzEJ2xASpYYIkZ3X+/Hm8/fbb2LVrF7Kzs+Hm5gaDwYD09HS0a9cOn332GZ5//nksWrQIEREReodLtVBQUIDDhw/Dw8MDHTt2RLNmzW6/Ed0xJkgNEyQ5OxFBcXExLly4gC5dumDkyJF45ZVXMGzYMHTt2hUHDhyoF/NjEtUVPs2DqJ5QSqFBgwbo0KEDXn/9daxfvx69e/eGu7s7YmJimByJbIhn7Imc1KxZs5CUlAQ/Pz/Mnj0bjRs31jskonqFCZLISbm5uWH79u16h0FUb3E8hoiIyAImSCIiIguYIImIiCy4q27zUErlAMi6bcHbawHgFxvsp644W7yA88XsbPECjLkuOFu8gPPFXNt4HxCRlpZW3FUJ0laUUket3TfjiJwtXsD5Yna2eAHGXBecLV7A+WK2Z7wcYiUiIrKACZKIiMgCJsg7E613ADXkbPECzhezs8ULMOa64GzxAs4Xs93i5TlIIiIiC9iDJCIisoBTzdWAUioMQB4ATwAZIpKsc0iVaDF6AQgAsElE4rTlUQCitGKBIuIwwyjWYnPU+lZKeQKYLSKzLKxziHrWYgwHkGceg7U6dYS6vk3MDtemq4jXYduzpZidpD1bawN2bc9MkNWkNaJQEZmsfd4HIFTfqMpTSvmjtDHEap+vAjA9TM4LwCYAcabv4EAqxebg9e0LIFz7RwiUxj9XRObDceq50mXv1urUgeraUsyO3Kat3VrgyO3ZUswO3Z6ttYG6aM9MkNX3NIAzZp/zlFL+jtKr0XgBGAPA1JBzzWKMAnBUK+NoLMXm0PUtImVPsVVKhZn+8cJB6llE4pRSvhUWW6xTlP7R1L2urcTssG3aSryAA7dnazE7eHu22AZgpd1aW34ndc1zkNXnidIuu0kuHCzZiEjFX3leZo3CF6Xxeiql5tV9dFWyFJvD1rf5PzTtV3ec2WpHrmdrderIde2MbZrt2YaqaAN2b8/sQdaOp94BWKM15mDT5wrnSFYopXxFJEOX4CqwFJuVog5V39pQjpeIlP1jdOR6tsJanTpUXQPO06bZnu2nYhuwwKbtmT3I6jOd8DXxAqD7P0ZLtF+BG8xPWCulZpoVydUnssqqiM0Z6ns2SoefADh2PWus1anD17WztGm2Z/up2AZQB+2ZCbL6NgLoaPbZ01HOh5nTxuCTRSRZKeWr/XrNQPlhEy9H+KWtsRabM9R3CMoP5ThyPQPW69Sh69rJ2jTbsx1YaQN2b8+cKKAGzC4dBlA6Nq5jOJVojSgev//S8zKdfDe7Qs0XQKyD/DEBYD02J6jvJBEJqLDMIepZKRWC0osaPAFEmV0BaLFOHaGuLcXsyG36NnVcKS5HrWOzdQ7ZnqvRBuzWnpkgiYiILOAQKxERkQVMkERERBYwQRIREVnABElERGQBEyQREZEFTJBETkSb7cTh90lUHzBBEjk4pZT5xMtP2+EQlfZZ4ZhEdyXeB0nk4JRSnqb5MS3dzG2D/Vu6QbzsmER3K/YgiRyYNotIkvY+BKVPVQeuODgAAAGLSURBVAjTlpvKhJvmzVRKmdafUUqFKKU2mYZQtc8hSql55ssq7tP8mNrnMG27MK28aV9XlVL+2nvTQ3WJ6g0mSCIHps0hmaG9j0Ppk+BjzSbt9gfQUZsyLBrAPO19Bkqn5ppl1hOcrO1jA0onpra4T/Njavv31R45FAtgjNa7jIM2sbVpH6bkSVRfMEESObdnAFwxJTL8/sR4LxFJNp87U0TGmNah+o//mQzAfKLnMyid1NrEtP8rNQ2cyNExQRI5IbNnDHpCe8qB9jKdS6z0aCJtaDXM0roK+7ydsofP8jwl1WdMkETOxZTcTMlsEwDz85H+lbYoXR4O4EqFpzeYylbcp7ly+wcQhNLHCRHVe0yQRA7MNHSqJTgAiDJ7X+78n/kFNAACzcuh9Ll+HSucJzT1BMvt0/yYFfYfhtJHJOWZlZmp9TxDoZ2ftG0NEOmHt3kQERFZwB4kERGRBUyQREREFjBBEhERWcAESUREZAETJBERkQVMkERERBYwQRIREVnABElERGTB/wN0C6J5tDas5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 467.995x155.998 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "latexify(fig_width=fig_width, fig_height=fig_width/3, columns=2)\n",
    "plt.plot(val_cost, c='k', label=\"Loss\")\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('cost')\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"supply_chain_training.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.64  0.3   0.02 -0.06]\n",
      " [ 0.3   1.44  0.32  0.3 ]\n",
      " [ 0.02  0.32  1.14  0.06]\n",
      " [-0.06  0.3   0.06  1.01]]\n",
      "[[-3.05]\n",
      " [-2.92]\n",
      " [-2.97]\n",
      " [-2.99]]\n",
      "\\begin{bmatrix}\n",
      "  0.64 & 0.3 & 0.02 & -0.06\\\\\n",
      "  0.3 & 1.44 & 0.32 & 0.3\\\\\n",
      "  0.02 & 0.32 & 1.14 & 0.06\\\\\n",
      "  -0.06 & 0.3 & 0.06 & 1.01\\\\\n",
      "\\end{bmatrix}\n",
      "\\begin{bmatrix}\n",
      "  -3.05 & -2.92 & -2.97 & -2.99\\\\\n",
      "\\end{bmatrix}\n"
     ]
    }
   ],
   "source": [
    "print(np.round(P_sqrt_train.T.dot(P_sqrt_train), decimals=2))\n",
    "print(np.round(q_train, decimals=2))\n",
    "def bmatrix(a):\n",
    "    \"\"\"Returns a LaTeX bmatrix\n",
    "\n",
    "    :a: numpy array\n",
    "    :returns: LaTeX bmatrix as a string\n",
    "    \"\"\"\n",
    "    if len(a.shape) > 2:\n",
    "        raise ValueError('bmatrix can at most display two dimensions')\n",
    "    lines = str(a).replace('[', '').replace(']', '').splitlines()\n",
    "    rv = [r'\\begin{bmatrix}']\n",
    "    rv += ['  ' + ' & '.join(l.split()) + r'\\\\' for l in lines]\n",
    "    rv +=  [r'\\end{bmatrix}']\n",
    "    return '\\n'.join(rv)\n",
    "print(bmatrix(np.round(P_sqrt_train.T.dot(P_sqrt_train), decimals=2)))\n",
    "print(bmatrix(np.round(q_train.T, decimals=2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0. 0.]\n",
      " [0. 1. 0. 0.]\n",
      " [0. 0. 1. 0.]\n",
      " [0. 0. 0. 1.]]\n",
      "[[-3.]\n",
      " [-3.]\n",
      " [-3.]\n",
      " [-3.]]\n"
     ]
    }
   ],
   "source": [
    "print(P_sqrt_baseline.numpy().T.dot(P_sqrt_baseline.numpy()))\n",
    "print(q_baseline.numpy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mean values over episodes\n",
    "n_episodes = len(u_behav)\n",
    "u_print = np.sum([np.vstack([u.detach().numpy().flatten() for u in u_b]) for u_b in u_behav], 0)/n_episodes\n",
    "x_print = np.sum([np.vstack([x.detach().numpy().flatten() for x in x_b]) for x_b in x_behav], 0)/n_episodes\n",
    "u_print_bl = np.sum([np.vstack([u.detach().numpy().flatten() for u in u_b]) for u_b in u_behav_bl], 0)/n_episodes\n",
    "x_print_bl = np.sum([np.vstack([u.detach().numpy().flatten() for u in x_b]) for x_b in x_behav_bl], 0)/n_episodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAACeCAYAAACl1tpAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfXgb1Z0v8O+ZkWRbL7Ysy44dSgiOcVhCeVFeoC2UUJwGmnSLIQllmybQ3uvw3KUte+mSve3tdsN9nt2mC2Fz297W3gU2KfdSIOBSQim1y3tbiF+AlgTqNxKSOIlfZNmS9WJp5tw/JBnZnpElWy8j6fd5nvPE0oxG58T6+aczc+YcxjkHIYQQUmiEbFeAEEIIyQZKgIQQQgoSJUBCCCEFiRIgIYSQgkQJkBBCSEGiBEgIIaQgUQIkhBBSkHTZrgAhhBAyH8aYFUATABfnvEVlny0AXACsAAY4593xjkkJkBBCSC5YE29jJEFu4JzvijxuA7Ah3mvoFCghhBDN45y3I9y7U7MNQH/MYxdjzBHvmHF7gHa7nS9fvjzhChKS77q6ukY455XJvIbiiOSzZGNitWDiE1xS3NaHwFEA/pinWtROdyqwYmaCdAKwxXtB3AS4fPlydHZ2JvjehOQ/xtiJZF9DcUTyWbIx4WYyfmS+UHHbTe4P/JzzuKc6k2SNt5GuARJCCMkcARBLROVt7kUdOTr4JcoGYGCeqhBCCCGZwQQGsURQLIv0JIAVMY+tNAqUEEKIdghYULJjjDUgPKrTyhhzcs4PRZ7vAnAD59zFGGuL7AcAe+c7JiVAQgghGcMYIBqST4CRUaDtCs+vjvn5UDLHpARICCEkcxiDoFe5BphhlAAJIYRkDGOAqNfG8BNKgIQQQjKIQRBZtisBgBIgIYSQDGIMdAo0Fw0NDeHgwYPo6emBx+OB2WxGfX09du7cicrKpCYH0ZR8bRfRhnz5fOVLO7KOMc2cAtVGLTSuo6MD27dvx8qVK/H+++/D4XBg06ZNcDgcOHbsGOrr67F9+3Z0dHRku6pJydd2EW3Il89XvrRDKxgDBJ2oWDKOc65aVq9ezQtdc3Mzr66u5g899BB3Op2K+zidTr5v3z5eXV3Nm5ubM1zDhcnXdqUbgE4eJ2aUSiHGUb58vvKlHemUbExcUmbh7/319YplIfG1mFIwgStLEg+5hrgcCiX8mubmZl5bW8t7e3sT2r+3t5fX1tZmNAhkWeZTbheXQsGEX5Mr7QpMurkUnMrYeyaikBKgLMt85LUjPDCi/IdfjRY/X7Is84nJAA+FpIRfo8V2cB5uy/jkFJckOa3vk6hkY2KV1cKP3dKgWDKdAAvmGmDwTB+CHx0D0xdBf149dFXLwUT1LndHRwe+//3v4/XXX0ddXR0AwOVyoaWlBVarFU1NTXNeU1dXhxdffBHXXnstHA4H1qxJ5ZyuyqZcw/Cc7gMTRRRXLEVxRQ0EUf3XqtSuqPb2dthstumJm6NtzEq7vB6MnewDmABTuR0mWxUEnT7t70s+NvTCq+j80i6IZiOW370DtffcAUNFedzXaDVuxtwBfHTODVFgqCwvQWVZCURR/QqQluNkaMyPUyM+6ESGalsJKsuKIAjaGFWZEMYg6LRx9U0btcgAHgxM/zt1/M/wvf1bBM/0g0vKy3Ls378fu3fvnvHhT2RG/7q6Otx3333Yv39/aio+D1kKAQC4JME3dBKuni54h05OPz+bUruA8B+p3bt3w+FwoKmpCbt3756xPVvtApcx6RzCUP8xuIdOQw4FM/L+BAiOjAEAJI8X/T/4GV6q+xw++N5DmBodU32NVuMmJMkAAEnmODvqxbHjTpx1TkKKPD+bluMkJPPwvxLHqWEv/vyhC+fG/JAjz2udlq4BFkwCnC1eIhwaGsLzzz+PnTt3znhNQ0MDrNa4q2sAAHbu3InDhw9jeHg45fWeT7xEqNYuALBarejq6gIADAwMKH57zWa7KBFm33yJMJfiJl4izLU4yblEGBkFqlQyLa2nQHloCtL4MKCB3wf3eZSfjyTCqZMfQLSfB8PyT+LgwYNobGxEeXn80z1qbDYbbr75ZjzS/FPc843/tphqzyvkn1R8PpoIfcOnYSi1wVRTm1C72tvb0dbWhubm5jnbou16uOVn+NY37k5ZG5RM+ZTbFU2Ek2PDKDKVwrLkE9DpDWmtS7ZwWcboK2/F7XGl01jHu4rPRxPhh//2KKq+sB6X/uj7KKqqSFnc/OzfH8bd37hnMVWfwxdQPiMSTYRDTi/KTEU4r8qcsjj52b8/gr9NcTsAwBdQPmsVTYSDo17YLAacX2mEIGivjxPtAWpBWhNg4MM/QRo5lc63SB0pCOnccQSZgJ6eHqxbt25Rh1u3bh3eerkNnpM9KargAnEZU+Mj4LKUULsaGsITqTc3N2Pv3rmTqa9btw5/fPUljA8eT0dtE8c5Ap5xSMEp2C+8OLt1SZPhtjfQsfm/ZrsaqmR/AGefeRGBM0P49Gu/SFncvPTamzhxdnELwyVL5sCYJwCJ85TFye9eewsfnlH5IpdGsgyMjE8BAC5YYs74+8+PFUYC1FctBziHFrqA8uQ4uEpvaZq+CLqqC+DxeGCxWBb1fhaLBd5gCIbSikUdZz5SwAcp4I27DxN1KLZVJ9yuhoYGbN26FbfddhscDseMbRaLBb6pIIot85/SWgwpOIWgf552CQJMtvy9Adm69jIsa/oygk5XVt7fe/w0xjv/HHcfXakZy7+xAwBSFjfBgBdWc9GijjObfyoE/5RyzylKFBgqy0pSFifBgBfl5tSfnfAFQvAHla9dRulEBntpccrfOyVYeE1ALUhrAhTL7BDL7Ol8i4QFjv8ZoTP9itsESwUM518ModQOxhjMZjPc7sV9A3W73SivqoFl2cpFHWc+vpFBeM8eV9wmFptgrDofekv5vO1qaWlBf3//9LdZm82muJ/b7Yatcgms512Ykvqr8U2MqfYyRX0RzPZqFJeG25WvDDYrPvmTPVl7/1MHW/Hu1/9BcZvp4lpc9N2/xdKtN02Ppk5V3CyptGF5TemijjPb0JgXgyPKX4BLinSothlRajKkNE6qK8tRuzT1PbDTI16cdfoVt5mKdVhaUQKLUafh2GBxR+BnUsHcBqFkduKLqq+vx5EjR3DXXXfN2D963t/lcsFms2HLli2qx+7o6MCqVavSVvd4Zie+KLV2AcC2bdvQ2dmJ9vZ2dHd3Y9euXXO+1QJZbleBJD4tU0p8UbkWN7MTX1QuxkluJL4wLV0DZOH7GJWtWbOGJzKEORdMnTiK4GAvAPXEFzU0NISVK1diYGBgQRf0nU4nVqxYgZ6enrTPEegbPQPvmQ8BqCe+qFxql9/tgut0pF0aSnyMsS7OeVI3eOVqHJ1+/Dm8s+PbAOInvigtf76GXV6cHg73ANUSX5SW2wEAgyNenIn0ALWQ+JKNicurK/iLX71JcVvNA/836fhajILpAeqqLgCXQtBVLFVNfFFVVVXYtGkTDhw4gHvuSX4U14EDB7B58+aMTJBbVFoBKeCDwVwGvcWWN+0yGC0wlldCX2zUROIrRFU3fhbL7/4qrFddETfxTe+v4c+X1VyEwJQEi9GgmviitNwOAKgoLUJI4rCaDTnR45st3ANc2OhUxtgWAC4AVgADnPNuhX0aItsBwMXDK8krizdNTK5O4ZQKR44c4dXV1QlPgxTV29vLq6ureUdHR5pqtjj52q5MQQFNhbYQ+fL5ypd2ZEKyMXF5dQUf+u4diiXesSJJrTnmcZvKPltiHt8Xry7au0lEI9auXYs9e/Zg48aN6OvrS+g1fX192LhxI/bs2ZOR6cIWIl/bRbQhXz5f+dIOTWJsoTPBbAMQO5LRxRibewEW+B+MsdrIz/GH4cfLjoX0zVVNdDb4ffv2qc4GPzo6yv/1Xx/gtool/Mc/+WmGa7gwibbrwQcfLNhZ7pWAeoAJSfTz9cADD2j680VxMr9kY+LypXY++r+aFAuA4wA6Y0pT9HUA7pv1uBlAw+zjA2gAMAbgqfnqUjDXABeqqakJDocD+/fvx/3334/GxkasXbsWFosFbrcbHR0daG1txbKLPodL1/8E9uXrs13lhCTars2bN+O5556jb7QkKfN9vo4cOYJftraiYf21eOrxx3DN+huyXWVFFCepx1jc2yBGeHKDYJRuSK4FsBrAXsZYG+d8g2pdeIGMAk2F4eHh6RWh3W43LBYL6uvrsWPHDgy5TPjuQ6dQZGD4P/+0HBXW3Fm5INqud//0Pnp6h3HJxVVYtepi7Nixg1a6nqWQRoGmilrcbL15MwxBL0S9AfYL/wpMg9N2xaI4UZZsTFz5iSr+yrduU9xmve/HqsdijDUBsHLOfxh5/BSAf+ExA2Eig2SmB8cwxvYifK1QcSAM9QCTUFlZiXvvvVdlG3D1FWa8+Y4Hj/1qFN/aUZ3h2i1ctF13fLMTuopJWD9RjnvvvSzb1SJ5Qi1uOOcY/fADhKb88LpGNT+rz3ScfCMcJ+XnU5wsSPweYDxPAoidd87K544CtQGIfa4NwIDaAbX9lSvH7Gy0Q6djePnNCfR9pDxTg1Z9dMqLvg/D90l1/8mFMddUlmtE8h1jDObKGgDA5OhZyCpLk2nJiZNe9B0Px0nXn1wYG6c4SVb0FKhSiYdz7gLQxhhriNzqMJ0MGWNdjDEr57wFwBbGWFOkxwjOOSXATKipNGDz+vAp6UcPDSPe6WWtOfDEiemfQxLH460ns1gbUiiKzGXQFxshSyF4x7KwzFaSDjwZEychjsdbc2Syf60RBOUyD875Ic55e7TEPL86kiDBOf8h57wlUtTvAQQlwJTbcqMNpWYRR/t8ePMd5SWYtOajU160vTo047lnnh+kXiBJO8YYLFVLAQCTznOaXufxxEkv2mfHyeHT1AtM1gJ7gOlACTDFzEYRX94UvvXkP1tHEJxn1nYtOPDECcxeQ9MfkKkXSDLCYLTAYCoFl2V4Rs9luzqqDjypFifUC0wWJcA8tvGaMpxfY8C5kSCefyU7S9kkSqn3F0W9QJIplsi1QK9rBFJQe585pd5fFPUCk8QYmE5ULJlGCTANRJHhzlvCI9qefMGJcbfyatRaoNT7i6JeIMmU6Jyv4ByekTPZrs4cSr2/KOoFJokBEEXlkmGUANPEscqEKy8xwuuX8YvnR7NdHUXxen9R1AskmWK2h3uBvnEnQgFflmvzsXi9vyjqBSaDrgEWhDtvqYQgAC++MY6TZwLZrs4cr/5xRPVbbZQ/IOMPnc7MVIgUNJ2hCEZreAFt97B2eoGvvZlYnPyxg+IkIYxppgdIN8Kn0bKlRfj8Z8rwm9fH8ejTw/jHuz+R7SrNcNPnlsDrCyEQCA/UeefoOHr6PbjsklJcXGcBABhLRFx7dfz5ZAlJFZO9Gr5xJwKecUz5JmEoMWW7SqpxcvmqUqxcERsn9mxWM2cwBloRvlDcvrkCr3W40X3Mi+6jk3Csyn5AR9krirBrR+304x8/3I+efg+uucqOv7nl/CzWjBQqUaeH0VaJydFzcA8NwrasLuvr3c2Okx/FxMntjRQnyWOAoI0ESKdA06zMosPWm2wAgEefGYYk5c7N8YRkg8lWBSaICPo8mJp0Z7s6JNUYA0SdcskwSoAZsHm9FUvsepw8M4Xf/n4829UhRNMEUQdzxRIAgHt4MKdmVCIJ0sg1QEqAGaDXC9jZGL4+8PjhUXi82p/zkJBsMpZXQtDpEQr44Hdr+15akiQNDYKhBJghn7rCjEvqSjDhkXDoNzRajJB4mCDAbA+vqOKhXmD+YUy5ZBglwAxhjOFrt4Zvjj/8igtnR+ieIULiKSmrgGgoghScgs81ku3qkFSha4CFqe6CYlx/VSlCIY4DrRTQhMTDGIPFHp4o2zNyFrJMlw7yhiAql0xXI+PvWOC2f6kCBj3DH9/24GivN9vVIUTTiixl0OXQckkkAXQNsHBVWPVo3FAOAHjk6WHI800xQUgBY4zBUhlZLml0CLKk3Xl1SWI4GLggKpZMowSYBY0bbLCViej/KIBXj9B9ToTEU2SywGC0gMsSJjW8XBJJHCXAAlZcJOCrXwoPiPn5s8PwB7S/ZiAh2TTdCxwb1uRySSQJjNE1wEJ33ToLViwrgnNcQmsb3RZBSDz6EiOKLdbIcklns10dsigMXBQVy7yvZGwLY6wh8q8jzn5N0f3iHY8SYJYIAsPXtoR7ga1tYxh1BbNcI0K07ePlkkYRCvizXBuyYAzggk6xxH0ZY1YAGzjn7ZzzQwD2quz3FIAnOeftAHbFOyYlwCxaVWfEp640YyrI8diz2lwzkBCt0BUVo6QsvDKJW4OL5pJELXgQzDYA/TGPXbN7gbGPGWO1nPMN8Q5ICTDLdjbaodMxvPzWBPpO0LdaQuIx26sBxhBwuxD0TWa7OmQh4l8DtDPGOmNKU8wrrQBi58VzArDNOvqayL82AFbGmGIvMYoSYJZV2w3YvN4KAHj40DBN+URIHKLeAFN5+NKBlhbNJYnjiDsKdIRzviamtMxzOKvSY875AOe8G4CDMVY792VhlAA1YOtNNpSaRbzf78Mf3/FkuzqEaJqpYgmYIGLK60ZgciLb1SELwMEUyzxcmJnwbAAGZu0zMOs5FwBKgFpmKhFx++bwtY0DrSMIBum2CELUCKIOpooqAOFeIJ01yTVsQYNgADwJYEXMY2uklxerHTMTXi2ATrUDUgLUiM9/pgzLagw4NxLE4Vdo+RdC4jGWV0IQdQj5vQjQckk5hTMGWRAVS9zXce4C0Ba5vaEBMaNAGWNdjDFrZJ8nIrdBNAH4l8hzijI//TZRJIoMd95aiT0/Po2nXnDi+qtLYbXQr4cQJYIgwmyvwcS5k3CPnEGRxQqWheV0yMIsdNaXyO0PSs+vnm8fJdQD1JArLzHBcYkRXr+MXxym2yIIiafEWgFRb4A0FYBvnOIldzDITFQsmUYJUGPuvLUSggD89o1xfDQYyHZ1CNEsxhjMlR8vl8RlunaeExgDZ6JiyTRKgBpzfk0RNl5TBpkDjz5Dy78QEk+xxQpdUQnkUBCTtFxSTuDAgq4BpgMlQA26fbMdxhIBbx/zouso3exLiJqZyyWdo+WScgKdAiVxlJpFbLsxPMHBo08PIyTRMG9C1BhMFhiM5vBySc6hbFeHzIMDkJmgWDKNEqBGbVpvRXWlHqfOTuG3b4xnuzqEaFbstcBJ5xCkEE0sr2mMeoBkHnq9gJ032wEAjx8egccrZblGhGiXocSEInMZLZeUIygBknldfYUZq+pK4J6U8dQLtGYgIfFYKiPLJblGEJqiEdRaxekaIEkEYwx3bqkEY8Dzr4zhzBCthE2IGl1RCUrKwtfOPbRckqZxxhRLplEC1Li6ZcVYf1UpQhJw4Jcj2a4OIZpmttcAjME/MYag35vt6hBFDDIXFUumUQLMAV/9azuKDAxvvuPBez0U1ISoEfUGGK3ha+fu4cEs14Yo4QBkCIol0ygB5gCbVYdbPh8+tfPI08OQZbotghA15opqMEHA1KQbgUl3tqtD5mCUAElyvnRDOSqsOgycDODlt2gNNELUCDodTLYlAADP8CAtl6QxHIDMBcWSaZQAc0RxkYDtXwqf2nnsVyPwB2jeQ0LUGG3h5ZKCfi8CHrqPVmuoB0iSdt1aC+qWFWFsXEJrG90WQYgaQRBhqqgGQIvmag0Hg8QFxZJplABziCAwfG1LeCXs1rYxjIzRjBeEqDGWR5dL8sM3Tl8YNYPTKVCyQJfUleDTV5oxFeT4+bN0WwQhahgTwrdFIHxfIC2XpBWMEiBZuB2Nduh0DK8ecaPnuC/b1SFEs4pLy6ErKoYcCsLroi+MWkCDYMiiVNsN+OL1VgDh1SLo+gYhymYsmjt6FrJEc+pqgcSZYpkPY2wLY6wh8q8jkX3j7UMJMEdtvdGGMrOI9/v9+MPbnmxXhxDNKjKVQl9iApdouSQt4As8BcoYswLYwDlv55wfArB3nn13zVcXSoA5ylgi4m++WAEAONA6jKkgXd8gREnsorleWi5JEziYYpnHNgD9MY9dcXqBawC0zXdASoA5rOHTZVhWY8DQaAiHX3ZluzqEaJbBaEaRqRScy5gcPZft6hQ0zuOeArUzxjpjSlPMS60AYv/QOQHYZh8/khQ7E6kLJcAcJooMd95aCQB46jdOuNyhLNeIEO2KXgv0jtFySdkmy0yxABjhnK+JKS3zHMqq9CTnPKEeASXAHHflJSasXmWCzy/j8edGs10dQjRLX1yC4tJyAJyWS8qiRdwI78LMhGcDMBC7Q6THWMsY2wJgLYANjLFatQNSAswDd9xihyAAbb8fx/HT9M2WEDXh+wKjyyXRLUTZInOmWObxJIAVMY+tnPPu2B045y2c80ORQTIDANo45zOSZCxKgHng/Joi3HhtGWROt0UQEo/OUARjeXhOXQ8tl5QVnAOSzBRL/NdxF4C2yG0QDYgZBcoY64qM/Iw+dgBoALA1Xg9Qt+jWEE348iY7XjnixrsfeNF1dBJrLjVnu0qEaJKpYgl8rlEEJicw5fXAYKRYybRE7vlTEunZKT2/etbjbgCrlfaNRT3APFFqFnHbF8K3Rfzn0yMISdQLJESJqNPDaAvPqeum5ZIyjkN5AIw8Tw8wHSgB5pEvXGdFTaUep85N4cXX6bYIQtSYbFVgooigbxKBSVpfM6M4IMnKJdMoAeYRvY5hZ2P4+sYvnh+Fx0vTPhGiRBBFmCPLJdGiuZnFsfCp0FKNEmCeuepyM1ZdVAL3pIwnX6AlYAhRY7TaIej0CAX88E+MZbs6BUWWlUumUQLMM4wxfH1LJRgDfv3KGAaHprJdJUI0iQmzlkviNJ1gJix0FGg6UALMQ7XnF+NzV5ciJIXnCSWEKCsps0FnKIYUnILXRRNJZArnyiXTKAHmqa980Y7iIoa33p3En//izXZ1CNGk8HJJ4V7g5MhZyDJdN083DuoBkjSzWXVo3BCeJ/aRp4chyXSRnxAlReYy6IuNkKUQvE46Y5J2NAqUZMLNDeWosOrw4akAXnlr5lBvSoiEhDHGYKkKT5Q96TwHOfTxpPKyJNHySSnGQYNgSAYUGQTsuDl8W8Rjz47A55cx4ZHQ8sQQbrunD4dfppFvhACAwWiBwWQBl+XpleM9I2cw3PceRgaOgWfjr3O+0lAPkKZCy3PXrrHg8Csu9B734wctg+g74cekL/xJO9bnw+bry7NcQ0K0wVK5FKOTf4F3bBi+cSd49HogB2QpBFEwZLeCeYIDkDRyqZV6gHmOMeCqy0wAgHc/8E4nPwDwBehbLSEAwDlHKOALBwzwcfKb3k6xkiqceoAkEwZO+vHIoWG816u87IvPT0FNSGDSDffQKYQCftV9aKaY1JI1MgaBEmCeOjEYwLf3fhT3wjIlQEKA8TMnIM830IUSYMpwTqdASZpZjCKsFjHuPnQKlBDAUnUewOL/KaQeYGpp5RQoJcA8ZbPqsP9/Lsd16yyq+1APkBCgpLQc9gtXQldsVN+JrgGmTLgHyBVLplECzGMWk4i/u6MG/9C0FGUKvUFKgISE6QzFqLigHqaKJYrbqQeYWlqZCo2uARaAq68w45K6ErQ8MYQ3utzTz4ckIBji0OsyPwURIVrDGIOlcimKTBa4BmdeF6T7AFMnfBvEwrIdY2wLABcAK4CByMrvSvvYEF4R/inOebva8SgBFohSs4hvf70Gn7rSjJ88dg7eSO/P65NQZqGPASFRBqMF9gsvxvjZkwi4wwtLh6YCWa5VHuF8QTNRMcasADZwzndFHrcB2DBrHwfCifFQ5PEYANWbnekvXxKGhoZw8OBB9PT0wOPxwGw2o76+Hjt37kRlZWW2q5eQzzgsWHVRCb73b6fg88swGYXpdj336y6c/MiJR1vsGBy4MqfaRbQrF+NGEHWwLl0Oz8gZ+CfGUGwpnxsnzXac7qc4SRYHIC+sB7gNQH/MYxdjzDGrF2gDsBXArshjp8I+0+gaYAI6Ojqwfft2rFy5Eu+//z4cDgc2bdoEh8OBY8eOob6+Htu3b0dHR0e2q5oQq0WHH31vOXbdPIo7du6Ybtft267DP+/ZgVtvvjYn20W0JdfjJnpK9LjTizu+9rU5cXILxcnCLHwQjBXh059RToQT3seH5rw92kOMsKklv+gLVMvq1at5oWtububV1dX8oYce4k6nU3Efp9PJ9+3bx6urq3lzc3OGa7gw+dqudAPQyePEjFIpxDjKl89XvrQjnZKNieoLHPw7D/sVC4DjADpjSlP0dQDum/W4GcAWtfcBsBeAI15dKHDjaG5u5rW1tby3tzeh/Xt7e3ltba3mgyBf25UJlADnly+fr3xpR7otJAHubvEplnjHAtAE4L6Yx0+pJTgAW+ZLfryQEuDU1BQfHh7mXq83of2PHDnCq6urE/7wR/X29vLq6mre0dGxkGomTQoGuPvMce6fcHJZlufdP1faJYeC3HfiGJ8aOZ1QuzKlkBLghDvI//d/9POX3hjikpTY70Crn6+x8SB/pHWEdx2dzPk4kfyT3Hf8PR4cH9ZEbCQbE0uWXcn//qdexTJPArQCaI553KaynwNAbeTn2ujPSqVgrgFOTExgYnwcZwYHMXj6NHw+5fkxo/bv34/du3ejrq5uxvOHDh1CS0sLdu3ahfb2uaNr6+rqcN9992H//v0prb+awMQYfM5zmDjZC9eHRxFwj0U/BIrU2hXr0KFDc9qW6XaFJkYQPPch/APvwPveawiODsZtF0m9N7ucePJXp/GPP/wAd3yrGy//fnjeORy1Gjd/eNeLF1534wcPD+E7+8+g+5g3Z+MkOHwKwXPH4fvgLfg+eBOhiZGciw1JlhVLPJxzF4A2xlgDY6wB4VOcAADGWBdjzBoZBfq7yH79ALo45wNqx2Tx/uPWrFnDOzs7k2yaNo2OjGB8fHzGc8XFxSi32VBSUjLj+aGhIaxcuRIDAwMoL/94BG13d/haqsPhAACUl5djbGzumnpOpxMrVqxAT09P2keHeUfPYvLcRzOe0xUbYaw8DwazFYx9fI+fWrtiuVwubN26Fbt370ZDQ8OMbZlsV9B5Bv7+t2c8JxSbYFh6EXS2mhntynmjEb8AAAgjSURBVCTGWBfnfE0yr8nVOHrhpXP45/09M567cJkRd355Ga77lB2CMPN3oOW4ee6Vcfz8uZnvueJ8A7Z+3oor/6okp+IkcPIDTJ3pn/GcaLHBcN5FEC0VGY+NZGNiyTIHv/3vX1Pctv+blqTjazEKpgeoxO/3K/YIDx48iMbGxjkffqfTiebm5unHNpttOrhj2Ww2NDY24uDBg+mrfBwhv1exR6jWrlidnZ3YsGGD4rZst0v2T1KPMMs+/Mir2iPMtbjpPzml2CPMxTiR3M6c6RFyziFJsmLJtLTeB+jz+eCemIAWfhWBgPqNrNFEKIoiqmtq0NPTg3Xr1s3Zr6GhYca3PafTOf2tdra1a9ei4/evY+JU3+IrH0e8JVyiiZCJOpR+YoVqu6K6u7uxZs0axT9OUWvXrkXn6y/BN6t3lmpyQP0UdTQR4sR7KL7gUugrlqa1Ltni9gTxH//vBCbcoay8/+BZ9c9WNBGaTSJ2330R1n+6MmVx8/gzb0K3ZHjxDYhx+pz6ag/RRGgxCfi7HertiEoqTvpSHyeSd0J9WyQRMn0RiuvXQmcqS/n7p4Ic0sbMOmlNgBMTE5j0eNL5FiklSRLcExPweDywWNQnkQaA3bt343e/+53qdovFgvHxMQQmnKmuZtK4FIJv9GxC7bJarXG3WywWuF1OhJxnUlnFhZFCmDo7kLcJ8Ohf3HjmeQ38P8fhmZTwxLOnsf7TlSmLm7ND4/j925Opruq83JMyfv1aYvGfUJyMORFyDqayignjwQCCQyegu/CyrLx/XBzQyHKA6U2AdrsdJpNJE2tpud3uuANfGGMwm82wVVTAbDbD7Xar7nvo0CHcdtttqt9io+9Xbq+C5bzaRdV7PlMeFwLjcZIsE1BkscJcvTxuu1paWmCz2TAwMICOjg6Mjo6itrYWtbUz6+92u1FWtRTFtZenshlzSB4XgkMn1HdgDGKpHcXLL01rPbJp7RXleOD7l2LCPc9adWny7tFxPPviWdXteh3DVattuPeu8ECRVMXNqots+OZX7AuvuILu9714o9urut2gZ1izqgT/5dYK/Pc3UxknV6S0HQAQdA5Ccg2p7yCI0NlqULRsVcrfOxU4OORsrH2kIK0JUBRFmM3mdL5FwgKBgGICFEQRVqsVpaWlEITwJdH6+nocOXIEd91115z9u7u74XA4UFtbi4GB8OCi2R9+IDwLxqpLL0NxWWoDeTY5FFJMgEwQUVJRjRLbEghi+Nccr11NTU3TP3d0dGDDhg3q7frkFdBXnJfCVihggnICZAL0VRfAUFMLQV+U3jpkmSgyXOVQvw6VbiGJKybA4iIBt2xaii/ffB7KywzTz6cqbq751Cpc40jt342xCUkxARqLGTZ9tgxf+KwFphJx3nYkFSeXXQG9PfVxIvvcyglQZ0BRTS30VReAiRqe5ZIjK9f7lGj4fym9lBJf1I4d4enBHnzwwTmj2W644QbYbOHZd5xOp+pottbWVuzdu3fOtnRTSnxRau2K1d3djfb2drhcrjnfbLPZrkJKfFqllviicilulBJfVM7FSa4kvgjOC+QaoJYIojj9r1rii6qqqsKmTZtw4MAB3HPPPdPPOxwOxcCd7cCBA9i8eXNGJsiNtite4otSa1csh8OBrq4uxW2ZbBcT9ZEfKPFlS2lklZD5El+UluPGYgzHerzEF6X1OGG6SGzkWOL7GNdMDzAvZ7BQIssyn5yc5JIkJbS/lmeCiCXLMg+4XVwKBRPaP5faFRwf4dKUPyPvlygU0EwwsizzrnfHuNMVSPg1Wv18SZLM3/2Ll3u8oYT212o7OOdcliQedA1xOcGYT7dkY8JWfRn/yndOK5aFxNdiSsHcB8gYg9FoVO31zbZ27Vrs2bMHGzduRF9fYrcy9PX1YePGjdizZw/WrMnMvZyMMRjMZaq9vtlyqV260grq9WURYwyOy6xxe32zafXzJQgMl9WXqPb6ZtNqOwCACQJ0ZZU51uuLFR4Eo1QyX5U8/OaaStHZ4Pft26c6G/zo6Ch/8MEHc2o2+HxtV7qhgHqAi5Evn698aUc6JRsT5Us+ybfde1yxLCS+FlNy9StExjQ1NcHhcGD//v24//770djYiLVr14bv83G70dHRgdbWVmzevBnPPfdcxnpIi5Wv7SLakC+fr3xph6Zw7dwGUTBzgabC8PDw9MrWbrcbFosF9fX12LFjR06vCJ2v7UqHQpoLNFXy5fOVL+1ItWRjonzJJ/n1t7Uqbmv90UUZnQuUEiAhSaAESMhMycaEtepSvn7LM4rbnv3pyowmQDoFSgghJHN4eNpJLaAESAghJKO4RiYDpQRICCEkYzjnkKkHSAghpPBwOgVKCCGk8ITnAqUESAghpNAs4hQoY2wLABcAK4ABzvmcVYkT2SeKEiAhhJCM4RyQFtADZIxZAWzgnO+KPG4DsCHZfWIVzFyghBBCtCDcA1Qq89gGoD/msYsxNnt15UT2mRa3B9jV1TXCGIuzLDchBeeCZF9AcUTyXFIx4XH95cXXf3mt2krhxYyx2FkjWjjnLZGfrQif2oxyArDNen0i+0yLmwA554U7vw8hKUJxRMjHOOc3pvBw1sXsQ6dACSGE5ILowJYoG4CBBewzjRIgIYSQXPAkgBUxj60KIzwT2Wda3MmwCSGEEK2IucUBAMA5b4883wXgBs65S20fxeNRAiSEEFKI6BQoIYSQgkQJkBBCSEGiBEgIIaQgUQIkhBBSkCgBEkIIKUiUAAkhhBQkSoCEEEIKEiVAQgghBen/A8LElSH1RdYlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 467.995x170.18 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Add wharehouse and store nodes (to plot edges)\n",
    "A_in_g = np.vstack([[0, 0, 0, 0, 0, 0, 0, 0],\n",
    "                    [0, 0, 0, 0, 0, 0, 0, 0], \n",
    "                    A_in,\n",
    "                    [0, 0, 0, 0, 0, 0, 1, 0],\n",
    "                    [0, 0, 0, 0, 0, 0, 0, 1]])\n",
    "A_out_g = np.vstack([[1, 0, 0, 0, 0, 0, 0, 0],\n",
    "                     [0, 1, 0, 0, 0, 0, 0, 0], \n",
    "                     A_out,\n",
    "                     [0, 0, 0, 0, 0, 0, 0, 0],\n",
    "                     [0, 0, 0, 0, 0, 0, 0, 0]])\n",
    "\n",
    "\n",
    "def draw_policy(ax, u_draw, plot_title, max_norm=None, fig_width=None, fig_height=None, print_colorbar=True):\n",
    "    # Create graph\n",
    "    DG = nx.DiGraph()\n",
    "    DG.add_nodes_from(range(n + k + c))\n",
    "\n",
    "    weights = np.sum(u_draw, axis=0)\n",
    "    if max_norm is None:\n",
    "        weights /= np.max(weights)\n",
    "    else:\n",
    "        weights /= max_norm\n",
    "        \n",
    "    for z in range(m):\n",
    "        i = np.where(A_out_g[:, z] == 1)[0]\n",
    "        j = np.where(A_in_g[:, z] == 1)[0]\n",
    "        # Check if they are not empty\n",
    "        if len(i) > 0 and len(j) > 0:\n",
    "            DG.add_edge(i[0], j[0], weight=weights[z])\n",
    "\n",
    "    # Layout\n",
    "    layout = {0: np.array([0.0, 0.5]),\n",
    "              1: np.array([0.0, 0.0]),\n",
    "              2: np.array([1.0, 0.5]),\n",
    "              3: np.array([1.0, 0.0]),\n",
    "              4: np.array([2.0, 0.5]),\n",
    "              5: np.array([2.0, 0.0]),\n",
    "              6: np.array([3.0, 0.5]),\n",
    "              7: np.array([3.0, 0.0])\n",
    "             }\n",
    "    # Labels\n",
    "    labels = {0: \"\",\n",
    "          1: \"\",\n",
    "          2: \"1\",\n",
    "          3: \"2\",\n",
    "          4: \"3\",\n",
    "          5: \"4\",\n",
    "          6: \"\",\n",
    "          7: \"\"}\n",
    "    node_color = ['red'] * k + ['blue'] * n + ['green'] * c\n",
    "    edges, edge_color = zip(*nx.get_edge_attributes(DG,'weight').items())\n",
    "    nodes = DG.nodes()\n",
    "    plt.sca(ax)\n",
    "    cmap = plt.cm.coolwarm\n",
    "#     cmap = plt.cm.Greys\n",
    "#     nx.draw_networkx_nodes(DG, layout, \n",
    "#                            node_color='k',\n",
    "#                            nodelist=[0,1],\n",
    "#                            node_shape='s')\n",
    "    nx.draw_networkx_nodes(DG, layout, \n",
    "                           node_color='w',\n",
    "                           edgecolors='k',\n",
    "                           nodelist=range(2, 6),\n",
    "                           node_shape='o')\n",
    "#     nx.draw_networkx_nodes(DG, layout, \n",
    "#                            node_color='k',\n",
    "#                            nodelist=[6, 7],\n",
    "#                            node_shape='^')\n",
    "    nx.draw_networkx_labels(DG, layout, labels=labels)\n",
    "    ec = nx.draw_networkx_edges(DG, layout, \n",
    "                                edge_cmap=cmap, \n",
    "                                edge_color=edge_color, \n",
    "                                edgelist=edges,\n",
    "                                edge_vmin=0,\n",
    "                                edge_vmax=1,\n",
    "                                arrowsize=20, width=2)\n",
    "\n",
    "    # Plot colorbar\n",
    "    if print_colorbar:\n",
    "        sm = plt.cm.ScalarMappable(cmap=cmap, \n",
    "                                   norm=plt.Normalize(vmin=0., vmax=1.0)\n",
    "                                  )\n",
    "        cbar = plt.colorbar(sm)\n",
    "#     plt.title(plot_title)\n",
    "    \n",
    "\n",
    "latexify(fig_width=fig_width, fig_height=fig_width/2.75)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, gridspec_kw={'width_ratios': [0.85, 1]})\n",
    "\n",
    "\n",
    "\n",
    "# Maximum order u between baseline and tuned\n",
    "max_norm = np.maximum(np.max(np.sum(u_print, axis=0)), np.max(np.sum(u_print_bl, axis=0)))\n",
    "\n",
    "draw_policy(ax1, u_print_bl, \"Untuned policy\", \n",
    "            max_norm=max_norm, \n",
    "            print_colorbar=False)\n",
    "draw_policy(ax2, u_print, \"Tuned policy\", \n",
    "            max_norm=max_norm, \n",
    "            print_colorbar=True)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"supply_chain_policies.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcgAAACvCAYAAAB5J068AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVxN+f8H8NeptGovilGoZE9qYtJQ9jEzZBsUyRpjmLHNMGYQxtiXYSj7kOxqapC9aIYvGWthWiZGQiVEi+q+f39U51e6LTfVvZf38/G4D/fc8znnvu+5R+/7OeezCEQExhhjjJWkIu8AGGOMMUXECZIxxhiTghMkY4wxJgUnSMYYY0wKTpCMMcaYFJwgGWOMMSnU5B1AcSYmJtS4cWN5h8HeQ1euXEklItOqbs/nLpOntz1/mXQKlSAbN26MqKgoeYfB3kOCINx7m+353GXy9LbnL5OOL7EyxhhjUnCCZIwxxqTgBMkYY4xJwQmSMcYYk4ITJGOMMSYFJ0jGGGNMCk6QjDHGmBScIBljjDEpOEEyxhhjUnCCZIwxxqTgBMkYY4xJwQmSMcYYk4ITJGOMMSYFJ0jGGGNMCk6QjDHGmBRKkSBv3ryJCxcuyDsMxmRGRIiIiMDdu3flHQpjTEZKkSCnTZuGb775Rt5hMCYzQRDQp08fbNq0Sd6hMMZkpBQJ0srKCvHx8fIOg7EqqVevHp48eSLvMBhjMlKKBGltbY20tDSkp6fLOxTGZGZqasoJkjElpDQJEgDXIplSqlevHlJSUuQdBmNMRkqRIK2srAAAcXFxco6EMdlxDZIx5aQm7wAqo1mzZjh16hTs7e3lHQpjMps2bRrGjBkj7zAYYzJSigSpoaGBbt26yTsMxqqkTZs28g6BMVYFSpEgAeD8+fP477//4OHhIe9QGJPJw4cPER4ejt69e8PIyEje4TCGK1euqKuoqExUVVUdRUT6AAR5xyQHEgC38vLyxjo4OEi9B6I0CXL79u0ICwvjBMmUzo0bN+Dp6Yk///wTzs7O8g6HMaipqW3W09Pr1KBBg1fq6uppgvD+5UeJRCKkpKS0fPTo0RYAfaWVUYpGOkBBS9bk5GS8evVK3qEwJhNTU1MA4IY6TJG4WFpaPtfQ0Mh9H5MjAKioqJCpqelzAK3LLFOL8byVoq4eCQkJco6EMdnUq1cPALirB1MkqioqKiTvIOSt8BiUmQdrNEEKgjBIEITxgiD4C4LQ/W32xV09mLLiGiRjstu+fbth0fPU1FRVZ2dnm+DgYN2yyqempqqWt74qaixBCoLQHkACEW0iIh8AB95mf5wgmbLS1NSErq4u1yAZq6QVK1aYjBo1Shw6zcTEJB8A3N3dM8rapqhMTEyMenXFUZONdIwADAbgU7j8VBCE9kT0d1V2ZmBggH/++QeWlpbVFiBjtSU8PBzm5ubyDoMxhZeamqpa1W3d3d0zJk6c2HDjxo1J1RFLjSVIIjoF4FSxl4ykJUdBEMYDGA8AFhYW5e7TxsamOkNk7K3Icu62b9++NkJiTOmtWbPG5Jtvvkkt/lpMTIy6vr5+fnBwsG5cXJxGx44dM11cXDJrOpZaaaQjCMJSAFJ7+hdegnUkIseiezVlOXHiBObOnVsTITImM1nO3TNnzmDHjh21E5gCIiIQvfdtQlglJCQkaBZdLi1y9OhRPQcHh1fu7u4Z1tbWOWFhYWXea7SyssqJjIzUro5YarwfpCAIgwDsq+ql1eL++usvLFq0CHPmzIGGhkY1RMdY7QgMDMTRo0fh7e0t71BqHRHB1dUVTZo0ea9/JCiy0aNHN7p161a1JJWytG7dOnPbtm3/VVQuPT291CXW06dP6+3ateseAMTFxWlYWVm9Lmt7Y2Pj/NjYWA0XF5fM1NRU1dDQUL34+Hj1RYsWPZY15ppuxdoewN9E9LcgCE0FQWj6NvuzsrICESExMbF6AmSslhTN6PE+1qIuX76Mc+fOYefOnbh9+7a8w2FK6Pnz56pFtcrDhw8bfv755y+K1nl4eJRqmGJoaJgHFDTc6dChQ5U7z9dYDbIwOZ5GQeMcoOAepGH5W5WvqC9kXFwcbG1t3zpGxmqLqakp8vLy8OzZMxgavtV/A6Xj5+cHHR0dEBF++uknBAQEyDsk9obK1OzkJTU1VVVfXz//zdcjIyO1jYyM8hITE9VjYmLUW7Zs+RoA0tLSVG1sbHKq471rrAZJRH8TkSERWRU+3vqvAs8LyZRVWYMFfPnll5g8ebI8QqoV6enp2Lt3Lzw9PTFx4kTs2bPnne6q9eOPP2Lv3r3yDkOpOTg4vCrekjUyMlK7W7duYo3Rzs4uMzIyUtvFxSWzXr16+V27dn1RlBwBID4+XqO6GvAozUg6AGBiYgI9PT08evRI3qEwJpOiBPnmYAHHjx/HiRMn5BGSzIgIsbGxMm2za9cuZGVlYcKECZgxYwbU1dXx888/11CE8hUfH49Fixbhq6++QkZGmd31WAU8PDzSd+zYIVao3N3dM2bMmCG2at24cWNSUX/I0NBQPUdHx8y36RpSHqVKkIIg4MmTJ1i8eLG8Q2FMJp06dcK9e/fQoUOHEq/Hx8fj7t271fpeRISrV68iOjoamZnV1xL+t99+Q7NmzXDs2LFKx+Hn5wcnJyfY29vDzMwM48aNw86dO9/JdgTbt2+HIAhIS0vDunXr5B2O0ipeG6yIoaFhXlxcnMaTJ09UASA4OFh38uTJJS7THD16VO/KlSs6VRlAQFCkRgOOjo4UFRUl7zDYe0gQhCtE5FjV7aty7ubl5WHUqFFQU1PD5s2boab2dk0CMjMzERgYiHXr1uHGjRvi6+bm5mjatCmsrKzQtGlTNG3aFH369IGxsXGl9y2RSNCyZUvcvXsXLVq0wI0bNyqM99y5c+jSpQu2b98utt598OABrKysMHr0aGzcuLFKn1MR5efnw9LSEnZ2dhAEAX/99RcSExOhp6dXK+8v6/l7/fr1RDs7u9SKS8pPcHCwbnkj57wpNTVVNTIyUluWbQDg+vXrJnZ2do2lrVOqGiQAHDt2DAMHDkReXp68Q2Gs0ogIixcvxvHjx8XX9u/fj4CAAOzYsQNeXl5V3vf9+/cxa9YsNGrUCOPGjQMA+Pv7Y8+ePVi0aBF69+4NNTU1nDlzBr6+vvDy8kK3bt0gkUgq/R4hISG4e/cuPDw8cPv2bWzdurXCbfz8/GBgYIAvvvhCfO2DDz7AqFGjsG3bNiQlVctgJwrhxIkTSEpKwujRo+Hr64v09HSsXbtW3mEpNVkTnYmJSb6s21RE6RLkw4cPcfjwYfz3n8I2umKsFEEQsGzZMvzxxx/ia1u3boW2tjZMTEywf/9+yFoDPXfuHAYNGoQmTZpg+fLlcHNzQ3h4OK5du4bx48dj6NChmDNnDrZt24bw8HD8999/yMzMhL+/P65fv17pxiREhKVLl6Jx48b47bff0LlzZ8ydOxcvXrwoc5snT57g4MGDGDlyJLS1S3avmzVrFvLz87Fs2TKZPq8i27p1K0xNTfH555/DwcEBffv2xapVq/Ds2bNq2X9aWlq17IfJRukSZPGuHowpk6K+kABw6dIlnDlzBj/++CPu3r2LBg0aYOjQoeUmnSJEhNmzZ6NLly44e/YsZs6ciYSEBBw8eBBdunRBefP7aWpqYuzYsWjbti3mzp2L3NzcCt/vzz//xMWLFzF9+nSoqalh5cqVePLkSbkJbseOHcjNzYWPj0+pdY0bN8aIESOwadMmPH4sc99thZOSkoKQkBCMGDEC6uoFt7nmz5+PZ8+eVUstMiYmBk2aNEFgYOBb74vJRmkTJHf1YMrG1NRUbMXq6+sLQ0NDeHh4wMjICJs3b0ZiYiLGjx9f7mACRISpU6diyZIlGDduHB48eIAlS5bINIi/iooKFi9ejPj4+EpdKl26dCmMjY0xatQoAICjoyM8PT2xcuVKqVdyJBIJ/P390aVLF7Ro0ULqPr///nu8fv0aK1eurHTciiogIAC5ubkYPXq0+Jq9vT369++PVatWIT09vZyty5eVlYUhQ4ZAU1MTbm5u1REuk4HSJUhzc3NoampyDZIpnaIa5KVLl3D06FFYW1vD0dER69atQ+/evTF79mzs27cPW7Zskbq9RCLBxIkTsXbtWkyZMgX+/v7Q0tKqUix9+vRBp06dsGDBgnJbukZHR+OPP/7A5MmToaOjI77+008/gYgwZ86cUtucOnUKCQkJmDBhQpn7tbGxwdChQ7Fhwwakpip0W5FyERG2bt2KDh06oFWrViXWzZ8/Hy9evMDq1aurvP9p06bh1q1b2LlzJ88GIwdKlyBVVFTw4YcflnsZiTFFZGpqipSUFCxYsABGRkZYuHAhVq5cCRcXFyxduhRTp05Fjx49MGXKFNy6davEtvn5+RgzZgz8/f3x3XffYc2aNW/1f0AQBCxZsgTJycnldklYsWIFtLS0MGnSpBKvW1pa4ptvvsGuXbvw998lh1n28/ODqakp+vfvX24Mc+bMQWZmJtasWVNmmdzcXMTFxVXqUrA8XL58GdHR0SVqj0Xatm2LQYMGYc2aNXj69KnM+z506BD8/Pwwc+ZM9O7duzrCZTLibh6MoXa6ebx69Qo3btyAs7MzfvrpJ3z//felyjx+/Bh2dnYwNjbG5cuXoa2tjdzcXHh5eWHv3r2YP38+5s6dW20/ED/99FNcuHABCQkJMDAwKLHuwYMHaNq0KSZMmIBffvml1LbPnz+HtbU1WrdujTNnzkAQBCQlJcHS0hIzZszAkiVLKnz/wYMH48SJE7h37x4MDAxARIiPj8eJEydw/PhxnD17FhkZGdDQ0EDr1q1hb28Pe3t7tG/fHm3bthUbAGVnZ+P+/ftITEzEvXv3xH9VVVXRqFEjfPDBByX+NTQ0rJZj6OPjg4CAACQnJ0vt0nHr1i20bdsWs2fPxk8//VTp/SYmJqJdu3awtbXF+fPnxXubZXkXu3nUlvK6eYjT0CjCw8HBgRiTBwBRVAvn7meffUZGRkb0/Plzio6OphcvXhAR0atXr+i///4jIqKTJ0+SIAg0ZswYysnJof79+xMAWrJkSfV94EJXr14lAPT999+XWjd9+nRSVVWlf//9t8zt169fTwAoJCSEiIh8fX0JAMXHx1fq/a9du0YAaPjw4eTj40NNmjQhAASAGjduTD4+PuTv708zZsygbt26kZGRkbheRUWFbGxsyMzMTHyt6KGmpkZNmjShRo0akaqqaqn12tra9Nlnn9HLly+rdNyIiF6+fEm6uro0cuTIcssNGTKE6tatSykpKZXa7+vXr+mjjz4iPT29Sh9HWc/fa9euJRJRlCI/tm3bFl/0PCUl5epHH330PCgo6G5Z5VNSUq6Wt76sR+GxkHqc5J4Uiz8q+0fmxIkT5OjoSElJSZUqz1hFaiNB7tmzhwDQjBkzKDs7mwDQwoULiYioY8eO5ObmJpadM2cOAaA2bdoQAFqzZk31fdg3DBs2jLS1tSk5OVl87enTp1S3bl3y8PAod9vXr1+Tra0t2draUlZWFn3wwQfUq1cvmd6/X79+BIB0dXWpX79+9Ouvv1JsbCxJJJJSZSUSCSUmJlJQUBDNnTuXBgwYQKNHj6YFCxbQzp076dy5c3T//n3Ky8sTt8nLy6MHDx7QhQsX6MCBA7Rq1SqaNGkSqaioUK9evSgnJ0emeIv89ttvBIAiIiLKLRcdHU2CINCsWbMqtd/Zs2cTANq7d2+lY3nXEuTy5ctLxffRRx89r2i7oKCgu9HR0Tdkea/yEmSNzwdZU6KiohAXF4cGDRrIOxTGKqWoyX/Xrl0BFAwU0LJlSwDA7NmzS4xMM3/+fERERCAyMhJ+fn5Su0tUlwULFuDAgQNYtGgR1q9fD6DgPuLLly8xc+bMcretU6cOli1bhn79+mHQoEF48OCBzMOs7dixA3fv3kX79u1Rp06dcssKggBLS0tYWlrC3d29UvtXVVVFw4YN0bBhwxKv29vbY+zYsRg1ahR27doFFRXZmmRs3boV1tbW+Pjjj8st17JlSwwbNgzr1q3DtGnTUN7k2idPnhRbKA8ZMkSmeN4VbzOuqru7e8bEiRMbbty4sVpGoVC6RjpAwbyQAPeFZMrjypUruHjxIoCCe5EaGhoYPHiw2PKxb9++6NOnj1heTU0NR44cwd9//12jyREo6Do1ZswYbNq0Cf/++y+ys7Oxdu1a9OzZE+3atatw+88//xyurq44cuQIGjRogM8++0ym9zcwMECHDh0qTI7VbcyYMfj5558RGBiIqVOnFlxSq6TY2FicO3cOo0ePrtS9zLlz5yIrKwvTp0/H1atXkZNTejamx48fY8SIEWjRokW5DZfedWvWrDHx9vYu0TcmJiZGXV9fPz84OFh3xYoVJpGRkTU6uXMRpUyQFhYWUFNT4wTJlMaCBQugr68PoKBjeXJyMm7cuCEOmZifn4/o6OgSHef19PRgb29fK/HNnTsXqqqqmDdvHnbu3InHjx/ju+++q9S2giBgxYoVEAQBPj4+bz2mbG367rvvMHXqVPzyyy8yTYKwfft2qKioYOTIkZUqb2tri/Hjx2PXrl1o37496tatizZt2mD48OFYvnw5Tpw4AS8vLzx//hz79u0rNfrQ+yQhIUGzaHLkIkePHtVzcHB45e7unmFtbZ0TFhamW9b2VlZWOdWVQJXnTC5GTU0NTZo04cECmFL4+++/ERISgnnz5sHX1xdPnjzB7t27MXPmTDx79gz6+vpISUlB69atxT6Ota1BgwaYMmUKli9fjlOnTsHBwUGmjukODg6IiYkRr+4oi6LknpKSgh9++AEmJiYV1tjz8vKwY8cO9OnTR6ZbPBs2bMDXX3+N69evi4+IiAjs3r1bLOPv74/WrVtX+fO8DScnp1Kz0A8YMODprFmzUjIyMlS6detm8+b64cOHp06ZMiUtOTlZrV+/fqW+/HHjxj0ZN25celxcXB1ra+tK9dVJT08vdYn19OnTert27boHAHFxcRpWVlZlzvhhbGycHxsbq+Hi4pJZNIB5XFychrW1dY6sY7UqfIIkIhw5cgQaGhro0aOH+HqvXr3e619ZTHksWLAABgYGmDp1Kvz8/JCVlQUvLy80btxY7BpQv359BAYGwtnZWW5xfvfdd/D390dycjLWrl0rczeI5s2b11BkNUtFRQXbtm3D06dPMXHiRBgbG2PQoEFllj9+/DiSk5Ol9n0sjyAIaN68OZo3b17i/mJaWhpu3LiBFy9eoG/fvlX+HO+y58+fqxbVKg8fPmwYEhKSULTOw8PDMjAw8F7x8oaGhnkAsGPHDkNvb+90d3f3DGdnZ5t3MkF+9913ICLcvHkTqqoFPy54vjWmDB49eoQTJ05g1qxZ0NfXR3Jysph4mjRpIpYTBAHDhg2TV5gAACMjIyxZsgSHDh3CgAED5BpLbatTpw4OHDiAHj16wNPTE0ZGRmJjqjdt3boV9erVk/lea1mMjY0VYhi5S5culTkxqa6urqS89ebm5nnlra9s7VGa1NRUVX19/fw3X4+MjNQ2MjLKS0xMVI+JiVEvmkcyLS1N1cbGJgcAiiZajoyM1Lazs5N5clSFvwepoqKCefPm4fbt29i/f7+8w2FMJmZmZkhISMA333wDAGJy/Ouvv0rdInj8+DFCQ0NlaixS3SZMmICTJ0+KP0TfJ9ra2ggNDYWNjQ26d+8OLS0tqY+goCCMGDGi1hsVvS8cHBxeFW/JGhkZqd2tWzdxFH87O7vMyMhIbRcXl8x69erld+3a9UXxSZbj4+M1XFxcSiTD3bt3G1alZavC1yABYNCgQWjdujUWLFiAL774Aqqqqvjzzz8xcOBABAcHo2PHjvIOkbEymZmZic83btyI6OhohIaGokuXLti5c6e4LigoCBMnTsT9+/fRqFEjeYT63jMyMsKpU6ewceNGZGdnSy1Tp04dTJ48uZYje394eHik79ixw7Co9ld4WVS8NFo80YWGhuo5Ojpmpqamqr7ZsKfI9u3bDRcuXPioeC2zspQiQRbVIgcPHoy9e/fC09MTxsbGePz4MeLi4jhBMqVx48YNHDx4EH/88Uepe+j9+vVDu3btUK9ePTlFx4CCHzS+vr7yDuO91bJly9dHjx6tVFlDQ8O8uLg4jWbNmuWYmJjkBwcH606ePDmlaH1wcLDu5s2bTfbv32/YuHHjHFlrkUqRIAFgwIABaNu2LRYsWIAhQ4agSZMmEASBu3owpWJqaoq0tDQ4OjqWuoxpbm7OMzYwhoJ7h8HBwboVNaopXrssuixbvJbo7u6eIWvDnOJkvgcpCELjwn/1ip7XhqJa5D///IM9e/ZAQ0MDFhYW3NWDKZV69eqBiLBv3z6ps81HRkYiLCxMDpExplhkTWwmJib5b5MMpZEpQQqC4AfAXxCEGQAcAXSvzmAq4u7uDjs7OyxYsAB5eXmwsrLiGiRTKkXDjHl6eko9dxctWiR1jkXGWO2T6RIrEU0QBKEJgKcoSI7/1khUZVBRUcH8+fPRv39/7N69GwMGDBBnaGdMGZibm6Nhw4bw9fUVx2Et7tdff5U6bRJjrPZV5R6kPhH9C+BQdQdTGUUNGRYuXIg7d+4o1bBWjHXu3BkPHjwoc72yjUTD2LuswkusgiDYC4IwQxCEAYIg6AFIEARBeg/aWiAIAubPn4/4+HgEBAQgOztb6sC/jCmq8PBwnD17Vuq6ly9fYuPGjbh27VotR8UYe1Nl7kE6ArgKwBjAQQAHAPQod4sa1rdvX7Rv3x4//vgjtLW18fvvv8szHMYqjYgwePBgjBkzRup6QRDw5ZdforLN3BljNUfq9UlBELoS0ZnCxfjC56cBbK61yMpRVIssGreQW7IyZSEIAiQSCTp06CB1vY6ODpKSkri7B2MKoKwa5CxBEI4LgnAcgE/hJVaFajnw2WefiX3J/vnnH3mHw1ilmZmZITe37KEpGzRoIPNA4Yy9a7Zv325Y9Dw1NVXV2dnZJjg4uMxprlJTU1XLW18VZSXIwUTUi4h6AZgF4DmALYVJc58gCGOrM4iqKKpF5ufn46+//pJ3ODUuPz8fWVlZ8g6DvaXMzEzk5ubi3r17ZZa5cuUKZsyYgdevZRoVi7F3xooVK0xGjRolTppcNIxcef0ci8rExMSoV1ccUhMkET0v9vxfItpMRF8UJszxqOXuHWXp06cPTExMEBcXh0ePHsk7nBpDROjXrx9sbW25W4uSu3fvHmJjY5GUVPaIV//88w9+/fXXcpMoY++q4gOVy8rd3T1j3bp1ptUVi8wj6RDRcyI6XV0BvA1BEDBz5kxIJBJMnjxZrrMg1KQtW7bgyJEjePDgATw9PZGfL3VM3ndCfn7+O/35mjVrhmnTpqFVq1Zllhk4cCBevnwJG5tS89My9s5bs2aNibe3d3rx12JiYtT19fXzg4ODdVesWGESGRlZK5MBK/x0VxUJDQ0FABw8eBC7du2SczQVy8yUbUqye/fuYdq0aejWrRv8/Pxw6tQp/PTTTzUUnXwREfr374+2bdsiNTVV3uHUCFVVVaxcuRInT54ss4y6uvp7Od0UYwCQkJCg+ebMHEePHtVzcHB45e7unmFtbZ0TFhZW5r1GKyurnOpKoErdy56I8O+//6J169ZIS0vDmDFj4OzsDGtr67fab1ZWFi5evIjs7GyYmZnBzMwMpqam4qAERITc3Fyoq8t2qTsyMhLdu3fHlClTsHTp0gobYkgkEnHW8q1bt8LCwgKRkZGYP38+OnXqhG7dulXtAyqoAwcOiD94+vXrh1OnTkFLS0tcv3fvXqxYsQL79u1T2g71ERERuHPnDsaPH1/u979+/XpkZGRg9uzZtRgde585OTnZVlSmd+/ezxYsWPC4qPzw4cNTp0yZkpacnKzWr1+/cv9TljehcnHp6emlfh2ePn1ab9euXfcAIC4uTsPKyqrMG/TGxsb5sbGxGi4uLpmpqamqkZGR2nFxcRp9+vR5Iet0V0pdg4yKikJSUhIaNGiA5ORk5OXloX///jI3bnj9+jXOnz8PX19fuLq6wsDAAF27dkWfPn3Qvn17NGjQAOrq6qhfvz7atGkDU1NTmJub4/nz5xXvvFB6ejo8PT0BAMuXL8fPP/9cYn1OTg7OnTuHlBRxphb4+fnhzJkzWLVqFSwtLSEIAjZu3IgWLVrAw8MDDx8+lOlzKrKMjAxMnToV9vb22Lt3Ly5cuICRI0dCIpEAAJ4+fYrx48fjypUr6NGjh9J+9v379+O7775Du3btcPv27TLLXbhwAefOnavFyBhTXM+fPxfnezx8+LDh559/Lk6g7OHhYflmeUNDwzygYLJlExOT/I4dO2ZW5d6kUtcgExISAACDBw9G9+7dMXv2bNy6dQvz5s0rlYDelJmZCT8/P4SFhSEyMhJZWVkQBAH29vaYMmUK3NzcYGhoiEePHomPBw8eICQkBGlpaQAKJr+dNWtWhXESEcaPH4+HDx/i/PnzWL9+PebMmQNDQ0OMGDECXl5eOHHiBF69eoWGDRsiPDwcKioqmDlzJnr27ImxY/+/0bCOjg4OHDiADz/8EMOGDcPp06ffieH2fH198fDhQxw6dAgdO3ZEUlISpk+fDktLSyxfvhxGRkbYs2cPMjMzMXr0aPTs2RMREREwNjau1P63bt0KV1dXudc8165di+7du2PAgAF4+PAhWrRoIbVcQEAAd/VgtaqyNTxp5c3NzfNk3b6yUlNTVfX19Us1TIiMjNQ2MjLKS0xMVC8+GXJaWpqqjY1NDlDQaCcmJkZ93bp1psXniaysGvvLKgiCAQpavD4jok018R5Fta1PPvkEr1+/xty5c5Gfn48lS5agR48e6NpV+oh4YWFh+PLLL8XLs+PGjYObmxu6dOkCQ0NDqdtkZmaib9++SEtLw86dO3Hs2DHMnj0bFhYW8PDwKDfOrVu34uDBg1iyZAnq1q0La2tr2NraYtKkSdDT00NaWhpGjBiBDh06YObMmXB1dYWFhQXU1NSwZcuWUn8oW7ZsCX9/f4wYMQI//vgjFi9eDEEQ8PDhQ4wYMQIAxG0EQcDkyZPFQRVkkZeXBwBQU1PD7du3cfLkSSQnJ4sPFRUVrFu37nRd4nAAACAASURBVK0vad+6dQtr1qzB2LFjxcmvp06din///RcrVqxA/fr1MWPGDHz66acAABMTE3h5eSEpKanCBCmRSDB79mwsW7YMkyZNwvr1698q1relpqaG5s2bA0C5LZJrIjlmZmYiOzsbRkZG1b5vxqqLg4PDq9TUVLHGGBkZqd2tWzexxmhnZ5cZGRmp7e7unpGamqratWvXEpdO4+PjNWbMmCE2YmjZsuXrXr16vQgMDDRctGjRY5mCIaIaeaBgto9vAYyv7DYODg4kixEjRpCBgQFJJBIiIjp37hzVqVOH6tSpQ+bm5pSamlqifHJyMg0dOpQAkK2tLZ09e7bS77Vw4UISBIF27txJt2/fpvXr11Pz5s0JAAUGBpa53e3bt0lbW5u6detGGzZsIAAkCAL17t2bunTpQqqqqhQaGkrh4eEUFxdHN27coP79+xMA2rZtW7kxjR07lgBQr169iIgoKSmJXFxcqFOnTtSpUydydnYmZ2dnOnToED169IiuXr1Ks2bNov/++6/Mfebn51NkZCR9+eWXZGJiQv/73/+IiGjbtm0EgNTU1MjMzIzatGlDjo6O9Pz580ofQ2kkEgl17tyZjIyMKCUlpcS6vLw86t69OwGg6dOnl1iXlZVVopw02dnZ4vc9YcIEys3NLTMOAFH0Fud7Zc/dxYsXU3BwMAGgtWvXllkuPT2dhg4dSr///nul9itNfHw8xcTEEBFRTk4OaWlp0axZs4iI6PHjx7Rnz54q75spFlnP32vXriUSUZQiPqKjo28sX768UvFt27YtPigo6G5KSsrVotcmTJiQXPR8zpw5/xWt++ijj55L20fhsZB6nGrsHiQRnQJQekbYapKbm4uAgADo6emJv7bbtGkDLS0tmJub48mTJxg7diyICBKJBP7+/mjevDkOHz6M+fPn4/r163B1da30+3377bcICwvDyZMn8euvv2Lu3Llwc3ODiooKvLy8cPdu6asL2dnZGDp0KLS0tODq6oovv/wSn376KZKTk3Hs2DGEhITA3t4egwYNwqBBgzB16lRoamrijz/+gLOzM9q0aYPY2Fip8aSlpeH+/fsAgLNnzyI+Ph4NGjTAqVOnsG7dOjg7OyMnJwcaGhr46quvYGlpCTc3NyxduhSWlpb4+OOPERAQgOfPn+PZs2e4dOkShg8fDjMzM7i4uGDTpk0wNDTEr7/+ijlz5uD27dvYsGEDgoKC8Pr1a1haWiI8PBx16tSBn58funfvjgMHDsjc1Wb37t04d+4clixZAhMTkxLrBEFAbm4uVFRU8Ouvv+Ly5cviOk1NTQDAvHnzMHDgQLG2W+Tp06fo2bMn9u7di6VLl2LDhg1yvxSdl5eHefPm4erVq1BRUSm3Bqmrq4tr1669VWveHj164IcffgBQ0DJ2+fLl4pWEn3/+GWPGjEFycnKV989YTZClIY2hoWFeXFycxpMnT1QBIDg4WLf4pdTevXtn3LlzR+OHH36o/+2338reWV6WXx2yPlBwibXcGmRhmSgAURYWFmX8PiotIyOD1NXVadCgQSVenzRpEgEQax7ff/89OTs7EwBydXWlO3fuVPo9MjMzadKkSfT48WMKDw8nJycnGjVqFPn6+tKIESMIALVp04YAkJGREb148aLE9l9//TUBoICAADI0NCR3d3fKyckpUSYlJYVatGhB2tradOLECerYsSOpqKhQ//79qVmzZtSwYUP6559/Smxz+fJlsrS0JHV1dVq4cCHVrVuXmjZtSi1atCBVVVUCINZU7e3tyc7OjgwNDcnQ0JA0NDTE9RU96tatSwYGBmWu19LSEmtozZo1IwDUp08f+vfffyt1fNPT06l+/frUoUMHys/PL7W+qMa9fPlyatKkCdWrV48SEhJKlFm3bh0BIC8vL3Ef//77LzVv3pzU1dXLrd0XhyrUIKty7ubm5lJmZib169eP/Pz8KrVNZR09epQGDBggHoeTJ0/SrVu3pJbNycmhqKioan3/qsjNzRWvALGqk/X8VeQaZNEjKCjorizlU1JSrsq6DVVQg5R7giz+kOUSa0JCAgGgTZs2lXj91atXZGpqSoaGhmRlZUUAyNjYmHbs2CHTf8TMzEzq0aOHeAkzODiY2rVrJybYnJwcatSoEQGg9u3bi/8WXe47cuQIAaDJkycTUcGl1tevX4v7l0gkNG/ePNq0aRM9ePCAGjduTHXq1CEAtH79enr16hXdvHmTjI2NydzcXEySRZ/PwsKCLl26RBKJhFauXEl6enqkqqpKQ4YMIU9PTwJADx48ICIif39/MjU1LZXg6tWrR82aNaO2bdtSgwYN6PDhwxQVFUWmpqakqqoqHi9jY+NS26qqqpKOjo74uePi4mj16tWkqalJmpqatHTpUqlJr7jJkyeTIAh05coVqd+vjo4O9ejRgyQSCd2+fZsMDQ3JwsKCZs6cSSEhIfT06VMiKrj8DYAmTZpEBw8epPr165OBgQGtWbOGDh06RIGBgXTkyJFyY6lKgqQqnruy+OOPP2jVqlVlHkuJRCJ+T4GBgdSqVSvxe6+soKAgmjRpUpmXqqtTfn6++D5Hjx4lPT09unv3LhER7du3jzp06EDJycmV2teRI0do586d4vLdu3cpKSmpzPJ5eXklLs3LepwU2buYIGvr8U4myIsXL9KaNWsoPj6+1LrQ0FACQHXq1KE5c+aUuLclkUgoMTGRDh8+TPPnz6cpU6bQmDFj6NWrV0REtHHjRurQoQNZWFgQAOrRowcZGhrSli1bSiXY7Oxs+vTTT0kQBDIyMiIA9PXXX9PDhw/JxMSE6tevT76+vlLjz83NpU8++YTGjh1LRAX/uevXr08DBw4U3ycnJ4fat29PmpqaZGJiIibJiIgIOnv2LLm4uJCZmZmYtPr160cxMTEUFxdHwcHB4mcq+twPHz6ko0eP0uLFi+mTTz4ha2trsaapqakp1oC3bNlCfn5+4h/lu3fvUmRkJP3++++0bds2WrFiBc2ePZtatmwp3pesW7curV27lrS1tcnS0pLs7e3Fz/Hy5ctSn//vv/8mFRUV+vLLL6Uen99++42MjIzo/v374vH6/vvvyczMrEQtuU2bNjRx4kRyc3MjAKSiokIWFhZ05swZsrGxEcu1bNlS6vsUqY0E+b///Y/mzJkjJvbK8PHxocaNG4vLy5cvp3nz5hER0dOnT8nV1ZW2bt1KRCWTjyx+/PFH6tChg9TvSRa5ubnk6+tL58+fJ6KC/x8LFy6kCxcuEBFRdHQ0mZqaUkhICBEV3CMdN24cxcbGEhFRcHAw9ezZU/wMP/30E3Xt2lVcnjhxIn322Wfi+3Xr1o169OghLn/44Yfi/XgioiFDhtDs2bPFZQsLC/L29haXGzRoQPPnz3+rz1wZOTk5dPbsWTEhZ2RkUGBgYLUmaE6QnCBLaNu2bYn/DMVJJBKys7MjAPTJJ5/Qnj17KD09nYiIli1bVqImpKmpSVpaWuTs7EzTp08nV1dXMjAwEGtHzZo1o3v37pUZR2ZmJnXs2JEAkIODAwEgS0tLUlNTIwDk7e1dIrFmZWVRRkYGPX36lH7//XeaO3cu9ezZk/T19alp06ZiUvr2229p3rx5NGvWLFq2bBnp6+vTqFGjaOTIkdSqVSsxfgsLC9qyZQulpaVV+tgV9/Lly1KXfSsrKyuLZs6cSSoqKtS2bVsCQM7OznT58mV68uQJxcXFkZ+fH+nr65OHh4fY4Cc/P586duxIpqam5SaLZ8+eUX5+Ph04cIC0tLTE76vos6urq5Obmxtpa2uLr9na2tLt27cJADk6OlKnTp2ofv36pKOjU+5nqY0E6e/vT6qqqpSWlkbTpk0jJyenCrchohINoby9valfv35EVHCe9+3bt0QtqqqKalY5OTklflhJU3Q+5+bmkqurKy1dulTcBwD6+eefxbgB0MqVK8X1I0eOpIsXL1Yqpo0bN9LIkSPF5VWrVomNjIgKGqVlZGSIy2fOnKGIiAhxedy4cbRo0SJx+ddff6XQ0FAiKjgHFy5cSMePHyeigqQVERFRLZd7c3Nz6cSJE3Tjxg0iIrp3716JRllXrlwhAHT48GEiIrp27RpZWFhQeHg4ERGlpaVRREREhd9DcZwglSxBoqAV6wEAJwEMqsw2siTINWvW0K5du8pcv3fv3hKJcOrUqRQREUF2dnbUrVs3AiA1WZqYmIj36ry9vSt1kiYnJ5ORkRGpqamJtSoANH78eMrPzyeJREJxcXG0adMmatq0qZh8i2o8bdu2JR8fH9q+fbv4H9TT05MmTpwovkfxy5wuLi60fv16sXYlb/Hx8ZSfn09r166lOnXqkIGBAQUEBNCUKVOoTp06NGTIEKpbty4BIHt7e/H+7Y4dOygtLY3i4+PFx19//UVhYWGUn59PixcvFhNvvXr1aMqUKZSXl0e3b9+mdevWicfn5cuX1KJFC2ratCllZGTQ69evydzcnDQ0NKhdu3b0xRdfiJe6y1Jbl1jz8vJIIpHQ5MmTSV9fX/aDTVSj9+yGDx9O3bt3F5d9fHxowoQJ4rKHhwd5eXmJy56enuTv7y/GlZ+fL/7Ik0gklJubW+GldkWwdu1aAkDXrl2r9DbZ2dni80WLFok/VHJzc6lu3bolro4cP36cnj17RkQFP0Kio6PF5ZiYGBoxYoR4++bw4cMEgC5fvkxERJcuXRKPcVk4QdZMgqyxZn1U0Ir1VE3t39/fH82aNcPw4cOlrh8yZAhCQ0Oxe/duDB48GKtXr0bHjh2hr6+PCRMmwMvLC5988glGjBgBTU1N1K1bt8qtHM3MzBAfH4+uXbuKo6N4enrigw8+wLJly7B9+3a4u7tj2bJl0NbWRvPmzTFo0CB07NgRjo6O0NUtPaxgQECA+Dw/Px9z5sxBamoq+vbtW+Zku/LStGlTAIC7uzumT58OLS2tEt/Lvn37xOdXr17F1atX4eTkBC8vL8yfPx8LFiwosT8DAwPMmjUL33//PRo1aoSAgAAMHTpUHJ+0efPmYl9CoGDwhJ9//hnR0dGoW7cuAODGjRswMjKCiopiDRZV9BlMTU3x/PlzvH79WuYhC2tyAIHBgwfjwoUL4nLxVuJAwbluaWkpLhc/TwVBKFFWEAS5txyurHHjxqFhw4aws7MDUDDa1e3bt7Ft2zYAwIQJEyCRSLBpU0GX7o8++ghmZmYICgoCABw+fBj29vYYMWIE1NTUcPbsWdja/v/IbT179hSfq6uro2XLluJyixYtsHPnTnHZ1dUVx44dEwe0j46ORmJiYs18cFYuoeDHh2JwdHSkqKioCsvdvHkTMTExsLW1Rbt27cosl5aWho0bN8LFxQUNGzaElZVVjf7BTE1NxapVq/Dzzz9j586d8PLywtKlSxEfH49p06ZBEATY2Ni80yOknD59GnZ2doiMjMT27dtx5MgRAP9/pQIo+MOpoqKCvLw8qKio4JdffoGWlhbGjh0LFRUVxMTEoEGDBli7di2mT58udumoSYIgXCEix6puX5lzd/Hixahfvz7GjBkDf39/TJgwAQ8ePEDDhg2r+rasBhARvL29YWJigpUrVwIAZs2aBSLC0qVLARQMA6mjoyMOzCGRSOT6Y0zW8/f69euJdnZ27+aMADK6fv26iZ2dXWNp65Tj590blixZgoiICDx48KDccsbGxmI/sNpgYmKCKVOmoH79+ujSpQuePn1a5sg876qiAdTd3d2hra2N1q1bl1k2PDwclpaW8Pb2hoaGBk6ePAk1NTU0a9YMADBnzpxaibm2HD16FFZWVhgzZgxMTQuGhUxJSeEEqWAEQcBvv/1W4rUlS5aUWJ4wYUKJZUW7UsGqh1ImyCFDhsDCwkLsCK9IzMzM8PXXX8s7DIXQs2fPEpeWKrJnz54ajEb+IiMjxVq0jY0Nhg0bVmK2EsbY/9u+fbvhqFGj0oGC8Vj79u3b9Ntvv33k7u6eIa180cwdZa2vCqX82RMREYHVq1fznHlM6RQf9SkwMLDEfaqa9uLFC4SHh2PlypU4duxYrb0vY7JasWKFSVFyBICicVnLS35FZWJiYmS7qV8OpUuQd+/excmTJ9G8eXOlaQDAWFxcHEaOHIno6OgSrxdN51XdimqqEokE3t7eaNGiBQwMDODm5oYZM2agT58+GD58ONLT0yvYE2O1KzU1tco1H3d394yqTGtVFqVLkIcOHcLNmzdLtAJjTNE9fvwYZ8+exatXrwAUJDBjY2OMGTMGmzdvxu+//45Lly5Vy3t9//33GDJkCICCe2OJiYmwsbGBr68vjh07huTkZMyfPx/79u1D69atERYWVi3vW1V5eXk4f/48MjMz5RoHUwxr1qwx8fb2LvHLLSYmRl1fXz8/ODhYd8WKFSaRkZHatRGL0iVIb29vAICjY5UbHDJW6zp16oT79+/DyckJQMGlVnV1dRw8eBDjx4+Hu7s7Bg0aJJYfOHAg+vTpIy7fvXsXubm5UvdNRAgPDxdro6qqqiUa/oSHhyMkJAQ//vgjevfuDTMzM8ybNw8XL16EgYEBhg8fjhcvXkjdd03JzMwUu5Ooqqriq6++EicUB4Dz58+LPyaKIyI8fvwYFy5cEAdaf/r0Kfbv34/HjwtmMpJIJHizdT4R4eXLl3jw4AFu3bqFyMjIChv5MflISEjQLLpcWuTo0aN6Dg4Or9zd3TOsra1zwsLCSveNK2RlZZVTXQlU6RLknTt3ABTcw2FMmeno6ODly5dYs2YNrly5UqKRkqurq9giWCKRwMXFBUZGRvj888+xdu1aREdH49WrV/D390erVq3g5uaGkydPAgAWLlyI1atXV/j+Dg4OuHLlCk6ePAk9PT1IJBJUppuVNLGxsXBychL7vGZkZODatWvIyckRyxARLl68CB8fH5ibm6NHjx7IyCi4pbRx40ZMmzYNQEFtu3PnzjAwMICzszO+/PJLDBgwAHZ2dtDT04OZmRmcnZ0REhICoGA+0SFDhuDWrVsAgKCgINStWxdWVlawtraGiYkJ6tSpA11dXTRq1Aht2rTBxx9/jHv37gEA4uPjcePGjVJJtbqkpqZi48aN+Pjjj2Fubo42bdqIf8cuXLgAX19fvHz5EkDNXXKvDCcnJ9tffvnFGABycnIEJycn2w0bNhgBQEZGhoqTk5Pt5s2bDYGCSYmdnJxsf/vtNwMASE5OVnNycrINDAzUB4D79++rOTk52R48eFAPAOLi4upUNo709PRSl1hPnz6t980336QW7kvDysqqzBk/jI2N82NjY0u03vzhhx/qV+XSrVLdxHv8+DGWL18OAGjbtq2co2Gs8lauXImkpCSsWrVKfC03NxeCIMDb2xv6+volyk+ePFl8XtRB/eTJkzh16hT++OMPAAU1r/z8fLRv3x47d+6Uafq2IpqamrC3twdQMLG3j48PLl68CCcnJyQlJUFbW1tqV6V//vkH+/fvh6mpKXx8fGBhYSEmIAD4888/8cknn0BVVRW2trZo3bo1bt68idu3b0NbWxuDBg3CqFGjoKOjA0EQ4OzsLO5bX18fYWFhCA8PR0REBAICAsR+zK6urrCysoKVlRUcHBwAAB9++CFu3rwpDmDQtGlTTJgwAcnJyRAEAQYGBqUe+vr64uTcq1evxrZt25CSkgIdHR08f/681Pchq1evXiEkJAS7d+/G8ePHkZeXh5YtW6JPnz54+vQp9PT0AAAXL17E/PnzMXPmTADA7NmzcfToUbi5ucHNzQ2dO3eucFLw98Hz58/FCZQPHz5sGBISklC0zsPDwzIwMPBe8fKGhobi/HcxMTHqCQkJVepMrVQJ8saNGwgLC4OBgQHq168v73AYq7QHDx4gPj5eXJZIJEhJSYGWllaFf4zV1NTQv39/9O/fHwBw7949nDp1CtHR0RgwYAA6depULYNPeHh4ID8/Hx9++CEAYObMmdizZw/MzMzQsmVLtGzZEkZGRggJCcG1a9cAAF5eXvDx8YGGhgZOnz4t7qt9+/bYu3cvbt68iZs3byIqKgrm5ubYsmULBg8eLCYIaTQ1NdGrVy/06tWrUnFraWmV6G9rb28vJv3KmDdvHj7//HPo6OgAKOjLm5aWJo7KJM3HH3+MDRs2AAAWLFggdttJTU3F1KlTERQUhFevXuGDDz7A1KlT4enpibZt25b6nqZOnYqvvvoKdeoUVLBatWqFq1evYuvWrVi3bh0EQUDbtm0xatSoGu8+dunSJXFSWw0NDSq+rKurKym+bGxsnF982dzcPK/4soWFRYlla2tr6fcHKiE1NVVVX18//83XIyMjtY2MjPISExPVY2Ji1IvmkUxLS1O1sbERL138888/Gk2bNs2uynsrVYLs0aMH2rdv/9a/7hirbW9e8vzrr7+QlZUlJj1ZWFpaYsyYMdUVmkhHR6dEB/hJkyahffv2iImJQXR0NHbs2IGXL1+iY8eOWL16NQYNGoQPPvhA6r7q1auHIUOGiI2FFJmpqamYjIkIo0aNQkREBPLzS/1NBlDw46b45N779u1Dly5dMGzYMOjr6+N///sfPDw84OnpiY8//rjCQQSKkiNQ8IPDy8sLr1+/xuXLl3H27FmEh4fj6dOn1fBJlYODg8Or1NRUscYYGRmp3a1bN/EmuZ2dXWZRf8fU1FTVrl27vig+yXJ8fLzGjBkzUoGCCZTd3d0zoqKiqnRPUqkSpEQiQWxsbI38cWCsNu3fvx+ampqlRmxRJJ06dUKnTp3E5aKGLtLGDn5XCIKASZMmYdKkSZXepnjXnTp16uDu3btvXaNXV1cXj39tjgamCDw8PNJ37NhhWJTkCvs+iv0fN27cmFT0PDQ0VM/R0TGzeEItzsTEJD84OFj3ypUrOqGhoa+L962sDKVppENEGD9+PPbu3Yu5c+fKOxzGKu3Vq1fo1q0bQkNDxdcaNWqEcePGASjo5qAMBEF4p5NjdXmXx1quDcVrgxUxNDTMi4uL03jy5IkqUFBjnDx5ckrRehcXl0wXF5dMoODSq6yxKE0N8vHjxzh27Bg6duz43o1vypRbRkYGsrKySiTCmTNnIiQkBHp6eoiKihIbnDDGgBkzZqQWXR4tr1zx2mVRK9U3E6yJiUn+sWPHEqRsXiGlSZBmZmZYtmwZrl27BiLiX2lMaZiZmeGvv/4Sl+/cuQMrKyvxPlZKSkpZmzL23pJ1TFUTE5P86hyHFVCiS6wAcO3aNZw4cYKTI1NaeXl56Ny5M8aPH4969eoBAJ48eSLnqBhj0ihNgvTx8UHDhg3F5uWMKYstW7aga9eu4kg4mzZtwoQJE0pMecUYUzxKc4k1KSkJpqamXHtkSkdFRQVqampic353d3cABQ3P1NXVuQbJ5CFfIpEIKioqNTN8kJKQSCQCgDKHL1KaGuTChQsRHh7ONUimdEaPHo0TJ04gNzcXixcvFoc4EwQBvr6+6N69u5wjZO+hyHv37hnk5OTUqakh9hSdRCIRUlJS9AHcKquM0tQg//77b/z555/ljm7BmCI7c+YM5syZg1atWonDos2aNUvOUbH3UV5e3rhnz55NzMjI8CYiIyhRZakaSQDcysvLG1tWAaVIkN9++y2CgoKgra2Npk2byjscxmTSuXNnfPbZZ7hz5w709PTQu3dvcZ1EIsGff/6Jdu3acR9DVmscHBxeA1hb+GBlUIpfDQYGBsjKykLr1q0rHLaJMUVCRGjUqBH09PQQFBQEd3d3aGj8/0QDly9fRufOnXH48GE5RskYk0Ypss3s2bORnZ3NU1wxpSMIAnbv3o1GjRrh2bNn+OKLL0qsd3JyQtOmTREQECCnCBljZVGKBPno0SOkpaXxFFdMae3btw+Ghobo0aNHidcFQcDw4cNx+vRpPHz4UE7RMcakUYoEefPmTQA8STJTTtnZ2QgODkb//v2hrq5ear2npyeICIGBgXKIjjFWFk6QjNWwsLAwZGRklDn1U7NmzeDk5IQjR47UcmSMsfIoRSvWGzduwNzcvMQcbIwpi/3798PY2Bhubm7lljE3N6/FqBhjFVGKBNm4cWMMGDBA3mEwViVGRkYYNWpUiYlx31TUL5IxpjiUIkH6+vrKOwTGqmz9+vWVKhcYGAg/Pz+Eh4dzdybGFAD/L2RMQaioqOD8+fOIiIiQdyiMMXCCZExh9O3bF7q6utwnkjEFwQmSMQWhra2NgQMH4uDBg8jKypJ3OIy99zhBMqZAhg8fjhcvXiA0NFTeoTD23uMEyZgCcXV1xdixY9GoUSN5h8LYe08pWrEy9r5QVVXF5s2b5R0GYwxcg2RMIcXGxuLSpUvyDoOx9xonSMYU0ODBg/H111/LOwxWDYgIFy5cgJeXFz777DMcPnwYeXl58g6LVUKNJkhBEAYJgtC98N/2NflejL1Lhg8fjosXLyI2NlbeoZSLiMTnX3zxBUaNGoW9e/ciNTVVjlEplmfPnsHNzQ2///47rl+/joEDB6Jx48Y4f/68vENjFaixBCkIggGAHkR0iogOAlhaU+/F2Ltm2LBh4lySiiY+Ph4bNmxAv3790KZNGzFJ2tra4vfff8ewYcNQr149fPjhh/Dz85NztLUvNjYWU6dORZ8+fQAAhoaGOHr0KJKSkpCYmIjg4GC0adMGNjY2AIALFy4gMjKyxI8NphhqspHOFwDiiy0/EwShPRH9XYPvydg7oWHDhujatStWrlwJJycn9OnTB1FRURg7dmypsmvWrIGrqyvOnTuHKVOmlFq/efNmfPjhhzh27Bhmz55dav3u3bvRqlUrHDp0CAsXLiy1Pjg4GI0bN8auXbvw448/4t69ewAKxkju1asXsrOzoaWlhYULF2L+/PmIiorC8ePHcfz4cdy/fx8A8Pr1a3Tq1Am5ubkl9v3NN9/A29sb9+/fR9++fUu995w5czB48GDcuXMHQ4cOLbX+p59+wqeffoorV65gzJgxpdavXr0abm5uOH/+PCZPnlxq7X0SdwAABLhJREFU/aZNm+Dk5ISwsDDMmjWr1PqAgAC0bt0ahw8fxoIFC0qtDwoKQpMmTRAQEIAVK1YgNzcXMTExUFNTw8CBA5GdnQ1NTU107dpV3KZfv37o169fic9w5MgRNGnSBHp6egCAli1bitOfDRs2DC1atMDcuXNLvT+rWTWZIA0APCu2/BSA0ZuFBEEYD2A8AFhYWNRgOIxVr5o+d319fbFy5Uro6OgAADQ1NdG4ceNS5bS1tQEAWlpaUtdramoCAHR0dKSu19DQAADo6upKXV80yLqenh46dOiAGTNmoFevXrC2toYgCCXKqqqqokOHDujQoQPmzp0r1oru37+Phg0bltq3vr6++B7S3ltXVxcAoK6uLnW9Ih2bovVDhw7FuHHjYGZmVqq8NPv27cOePXtw7Ngx5OfnA0CJmV0aNGgAY2PjSu2LVS+hpqr1giB8C+AZEW0qXPYHcLLwcqtUjo6OFBUVVSPxMFYeQRCuEJFjVbfnc5fJ09uev0y6mmyk8wwFtcgiRgASavD9GGOMsWpTkwlyPwCrYssGfP+RMcaYsqixe5BE9EwQhJOCIHQvfIlbsTLGGFMaNXYPsioEQUgBcK+M1SYAuHNV2fj4VKy8Y2RJRKZV3XEF525F7834+FRGjZ2/TDqFSpDlEQQhim9Cl42PT8XkeYz4+ykfH5+K8TGqfTzUHGOMMSYFJ0jGGGNMCmVKkJvkHYCC4+NTMXkeI/5+ysfHp2J8jGqZ0tyDZIwxxmqTMtUgGWOMsVrDCZIxxhiToiYHK68WgiAMwv8PW5fAo/GIU4mNR7Gxbgtf52NVqPBYGAFwAHCAiE4Ve73WjhF/JyXxuVs5inL+vu8UOkEWm1PSp3D5JIAe8o1KIZTqC8XH6v8VTs6dUDQwviAI6QAMa/sY8XciFZ+7FVCU85cp/iVWqXNKyisYRVH4a/LZGy/zsfp/RgB8ii0/LTwWtX2M+Dt5A5+7laIo5+97T9ETZKXmlGQA+FiJiOhU0a/sQkaFl6Jq+xjxd1I5fJyKUaDz972n6AlSGoOKi7BC7/2xEgRhKYBu5RSp7WP03n8nlcTHCQp5/r5XFD1B8pySlcfH6g2FDRr2FWvIUNvHiL+TyuHjJIUCnL/vPUVPkDynZOXxsSqm8N7M30T0tyAITQVBaIraP0b8nVQOH6c3KMj5+95T+JF0ijVrBiDe5H+vFc6x6YOCX5P+xVq78bGC+MflNAru0QAF93AMC9fV6jHi76QkPncrpkjn7/tO4RMkY4wxJg+KfomVMcYYkwtOkIwxxpgUnCAZY4wxKThBMsYYY1JwgnxHFI7TyJhS4vOXKSJOkO+O0/IOgLG3wOcvUzicIN8BhZ2Io+QdB2NVwecvU1ScIJVcYafipYXPu8s5HMZkwucvU2QKPR8kq1jhUFRPASwlIh6XkSkVPn+ZIuMa5LvBkf+4MCXG5y9TSJwg3w1PKy7CmMLi85cpJE6QSq5o1P/C53wPhykVPn+ZIuMEqfwSAKQJgtCdR/ZnSojPX6aweDYPxhhjTAquQTLGGGNScIJkjDHGpOAEyRhjjEnBCZIxxhiTghMkY4wxJgUnSMYYY0wKTpCMMcaYFJwgGWOMMSn+D8hQZSk6+ImfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 467.995x187.198 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "latexify(fig_width=fig_width, fig_height=fig_width/2.5)\n",
    "# fig, (ax1, ax2, ax3) = plt.subplots(1, 3)\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True)\n",
    "\n",
    "\n",
    "# # Loss function\n",
    "# ax1.plot(val_cost, c='k', label=\"Loss\")\n",
    "# ax1.set_xlabel('iteration')\n",
    "# ax1.set_ylabel('cost')\n",
    "\n",
    "\n",
    "# Storage\n",
    "ax1.plot(x_print_bl[:, 0], c='k', label=\"$(h_t)_1$\")\n",
    "ax1.plot(x_print_bl[:, 1], linestyle=\"--\", c='k', label=\"$(h_t)_2$\")\n",
    "ax1.plot(x_print_bl[:, 2], linestyle=\"-.\", c='k', label=\"$(h_t)_3$\")\n",
    "ax1.plot(x_print_bl[:, 3], linestyle=\":\",c='k', label=\"$(h_t)_4$\")\n",
    "# ax3.legend()\n",
    "ax1.set_ylim([-0.1, 2.7])\n",
    "ax1.set_xlabel(\"$t$\")\n",
    "ax1.set_ylabel('$h_t$')\n",
    "# ax1.set_title(\"Untuned policy\")\n",
    "\n",
    "\n",
    "# Storage\n",
    "ax2.plot(x_print[:, 0], c='k', label=\"$(h_t)_1$\")\n",
    "ax2.plot(x_print[:, 1], linestyle=\"--\", c='k', label=\"$(h_t)_2$\")\n",
    "ax2.plot(x_print[:, 2], linestyle=\"-.\", c='k', label=\"$(h_t)_3$\")\n",
    "ax2.plot(x_print[:, 3], linestyle=\":\",c='k', label=\"$(h_t)_4$\")\n",
    "# ax2.legend()\n",
    "ax2.legend(bbox_to_anchor=(1.04,0.5), loc=\"center left\")\n",
    "ax2.set_ylim([-0.1, 2.7])\n",
    "ax2.set_xlabel(\"$t$\")\n",
    "# ax2.set_ylabel('$h_t$')\n",
    "# ax2.set_title(\"Tuned policy\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"supply_chain_storage.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}