RL_SAC_train.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f1033e6c-3a8d-4022-bcd6-87f6e190253a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# import argparse\n",
    "import datetime\n",
    "import numpy as np\n",
    "import itertools\n",
    "import torch\n",
    "from SAC import Environment, SAC\n",
    "from torch.utils.tensorboard import SummaryWriter\n",
    "\n",
    "from replay_memory import ReplayMemory\n",
    "\n",
    "from ROMWithMSE2C import ROMWithE2C\n",
    "\n",
    "# from pcc_model import PCC\n",
    "from data_preprocessing import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f69d1257-0c6d-4b4b-a924-10882120c192",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f2480ac5-b50b-4250-af1c-f74ade2c8338",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "output_dir = './saved_models/'\n",
    "\n",
    "data_dir = 'data/'\n",
    "state_file = 'states_norm_slt'\n",
    "ctrl_file = 'controls_norm_slt'\n",
    "yobs_file = 'rate_norm_slt'\n",
    "\n",
    "state_data = state_file + '.mat'\n",
    "ctrl_data = ctrl_file + '.mat'\n",
    "yobs_data = yobs_file + '.mat'\n",
    "cond = 'SC'   # standard condition for rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5c39560c-5a90-4444-98e3-096b09d6fa4e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(600, 21, 64, 64)\n",
      "[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19]\n"
     ]
    }
   ],
   "source": [
    "nsteps = 2\n",
    "_, _, _, _, _, num_t_slt, Nx, Ny, num_well, num_prod, num_inj = prepare_data(nsteps, data_dir, state_data, ctrl_data, yobs_data, cond)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "00f715ac-be04-4601-9d58-b366f6517ad7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_channels = 2\n",
    "method = 'E2C'\n",
    "ode_steps = 0\n",
    "\n",
    "#################### model specification ##################\n",
    "learning_rate = 2e-4\n",
    "latent_dim = 20\n",
    "\n",
    "u_dim = num_well  # control dimension, gaussian 9 wells\n",
    "\n",
    "input_shape = (n_channels, Nx, Ny)\n",
    "perm_shape = (1, Nx, Ny)\n",
    "prod_loc_shape = (num_well, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "84be7f17-a99d-41dc-a257-92edf9ed7f0c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "epoch = 200\n",
    "num_train = 9000\n",
    "num_eval = 3000\n",
    "encoder_file = output_dir + 'e2co_encoder_dt_'+'_nt%d_l%d_lr%.0e_ep%d_steps%d_cond%s.h5' % (num_train, latent_dim, learning_rate, epoch, nsteps, cond)\n",
    "decoder_file = output_dir + 'e2co_decoder_dt_'+'_nt%d_l%d_lr%.0e_ep%d_steps%d_cond%s.h5' % (num_train, latent_dim, learning_rate, epoch, nsteps, cond)\n",
    "transition_file = output_dir + 'e2co_transition_dt_' +'_nt%d_l%d_lr%.0e_ep%d_steps%d_cond%s.h5' % (num_train, latent_dim, learning_rate, epoch, nsteps, cond)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "efd7c6fa-9991-428e-b928-0f16121f86bf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_rom = ROMWithE2C(latent_dim, \n",
    "                    u_dim, num_prod, num_inj, \n",
    "                    input_shape, \n",
    "                    perm_shape, \n",
    "                    prod_loc_shape, \n",
    "                    learning_rate,\n",
    "                    method, \n",
    "                    ode_steps,\n",
    "                    nsteps, \n",
    "                    sigma=0.0).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "779d1ff1-1d59-4ea8-a553-92b4abda4163",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "my_rom.model.load_weights_from_file(encoder_file, decoder_file, transition_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8e95dd14-4835-49ac-b669-9ee5ee3b0b8a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hf_r = h5py.File(data_dir + state_data, 'r')\n",
    "sat = torch.tensor(np.array(hf_r.get('Mole_frac_norm_slt')).transpose((3,2,1,0)), dtype=torch.float32)\n",
    "pres = torch.tensor(np.array(hf_r.get('Psim_norm_slt')).transpose((3,2,1,0)) , dtype=torch.float32)\n",
    "hf_r.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fb030b75-7d01-4a4e-8b5b-b300e6b36b7b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sat_t_seq = sat[0:1, 0:1, ...].to(device)\n",
    "pres_t_seq = pres[0:1, 0:1, ...].to(device)\n",
    "\n",
    "state_t_seq = torch.cat((sat_t_seq, pres_t_seq),dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7eba2411-fcec-481d-8f13-8f9b8ed384fa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "z0 = my_rom.model.encoder(state_t_seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6ff8acd8-3573-4c4d-b2a5-839fea2d24d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Agent\n",
    "agent = SAC(latent_dim, u_dim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2059013e-6734-4c24-ab52-f4c3526266f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Tesnorboard\n",
    "writer = SummaryWriter('runs/{}_SAC_{}_{}_{}'.format(datetime.datetime.now().strftime(\"%Y-%m-%d_%H-%M-%S\"), 'ROM_RL', 'Deterministic', \"autotune\" if False else \"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ad2e4d47-5b08-48c9-b436-8e2991f9be29",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Memory\n",
    "memory = ReplayMemory(10000, 123456)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "31a50308-a55a-439d-86c4-aa81594c98d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "state = z0\n",
    "# ROM based Enviroment\n",
    "env = Environment(state, num_t_slt, num_prod, num_inj, my_rom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "28dc6aa6-1c4b-4c69-82b3-d0fed432a206",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Episode: 0, total numsteps: 20, episode steps: 20, reward: [-151.72453]\n",
      "Episode: 2, total numsteps: 60, episode steps: 20, reward: [-151.72453]\n",
      "Episode: 4, total numsteps: 100, episode steps: 20, reward: [-151.72453]\n",
      "Episode: 6, total numsteps: 140, episode steps: 20, reward: [-151.72453]\n",
      "Episode: 8, total numsteps: 180, episode steps: 20, reward: [-151.72453]\n",
      "Episode: 10, total numsteps: 220, episode steps: 20, reward: [-151.72453]\n",
      "Episode: 12, total numsteps: 260, episode steps: 20, reward: [-151.64601]\n",
      "Episode: 14, total numsteps: 300, episode steps: 20, reward: [-130.74089]\n",
      "Episode: 16, total numsteps: 340, episode steps: 20, reward: [-91.47116]\n",
      "Episode: 18, total numsteps: 380, episode steps: 20, reward: [-41.13715]\n",
      "Episode: 20, total numsteps: 420, episode steps: 20, reward: [-9.975745]\n",
      "Episode: 22, total numsteps: 460, episode steps: 20, reward: [-14.153573]\n",
      "Episode: 24, total numsteps: 500, episode steps: 20, reward: [-31.169514]\n",
      "Episode: 26, total numsteps: 540, episode steps: 20, reward: [-38.4878]\n",
      "Episode: 28, total numsteps: 580, episode steps: 20, reward: [-40.30822]\n",
      "Episode: 30, total numsteps: 620, episode steps: 20, reward: [-41.6811]\n",
      "Episode: 32, total numsteps: 660, episode steps: 20, reward: [-43.177273]\n",
      "Episode: 34, total numsteps: 700, episode steps: 20, reward: [-44.499264]\n",
      "Episode: 36, total numsteps: 740, episode steps: 20, reward: [-45.48634]\n",
      "Episode: 38, total numsteps: 780, episode steps: 20, reward: [-46.172016]\n",
      "Episode: 40, total numsteps: 820, episode steps: 20, reward: [-46.551075]\n",
      "Episode: 42, total numsteps: 860, episode steps: 20, reward: [-46.732616]\n",
      "Episode: 44, total numsteps: 900, episode steps: 20, reward: [-46.750217]\n",
      "Episode: 46, total numsteps: 940, episode steps: 20, reward: [-46.311493]\n",
      "Episode: 48, total numsteps: 980, episode steps: 20, reward: [-43.350746]\n",
      "Episode: 50, total numsteps: 1020, episode steps: 20, reward: [-23.836678]\n",
      "Episode: 52, total numsteps: 1060, episode steps: 20, reward: [-19.718]\n",
      "Episode: 54, total numsteps: 1100, episode steps: 20, reward: [-19.24349]\n",
      "Episode: 56, total numsteps: 1140, episode steps: 20, reward: [-18.034922]\n",
      "Episode: 58, total numsteps: 1180, episode steps: 20, reward: [11.150725]\n",
      "Episode: 60, total numsteps: 1220, episode steps: 20, reward: [22.896896]\n",
      "Episode: 62, total numsteps: 1260, episode steps: 20, reward: [23.27326]\n",
      "Episode: 64, total numsteps: 1300, episode steps: 20, reward: [23.463556]\n",
      "Episode: 66, total numsteps: 1340, episode steps: 20, reward: [24.03093]\n",
      "Episode: 68, total numsteps: 1380, episode steps: 20, reward: [61.95368]\n",
      "Episode: 70, total numsteps: 1420, episode steps: 20, reward: [101.93341]\n",
      "Episode: 72, total numsteps: 1460, episode steps: 20, reward: [102.06423]\n",
      "Episode: 74, total numsteps: 1500, episode steps: 20, reward: [102.094475]\n",
      "Episode: 76, total numsteps: 1540, episode steps: 20, reward: [102.10586]\n",
      "Episode: 78, total numsteps: 1580, episode steps: 20, reward: [102.11045]\n",
      "Episode: 80, total numsteps: 1620, episode steps: 20, reward: [102.11246]\n",
      "Episode: 82, total numsteps: 1660, episode steps: 20, reward: [102.11353]\n",
      "Episode: 84, total numsteps: 1700, episode steps: 20, reward: [102.11413]\n",
      "Episode: 86, total numsteps: 1740, episode steps: 20, reward: [102.11451]\n",
      "Episode: 88, total numsteps: 1780, episode steps: 20, reward: [102.11473]\n",
      "Episode: 90, total numsteps: 1820, episode steps: 20, reward: [102.11487]\n",
      "Episode: 92, total numsteps: 1860, episode steps: 20, reward: [102.114914]\n",
      "Episode: 94, total numsteps: 1900, episode steps: 20, reward: [102.11469]\n",
      "Episode: 96, total numsteps: 1940, episode steps: 20, reward: [102.12722]\n",
      "Episode: 98, total numsteps: 1980, episode steps: 20, reward: [90.00903]\n",
      "Episode: 100, total numsteps: 2020, episode steps: 20, reward: [146.19852]\n",
      "Episode: 102, total numsteps: 2060, episode steps: 20, reward: [146.29536]\n",
      "Episode: 104, total numsteps: 2100, episode steps: 20, reward: [146.30205]\n",
      "Episode: 106, total numsteps: 2140, episode steps: 20, reward: [146.3053]\n",
      "Episode: 108, total numsteps: 2180, episode steps: 20, reward: [146.30696]\n",
      "Episode: 110, total numsteps: 2220, episode steps: 20, reward: [146.30785]\n",
      "Episode: 112, total numsteps: 2260, episode steps: 20, reward: [146.30833]\n",
      "Episode: 114, total numsteps: 2300, episode steps: 20, reward: [146.30864]\n",
      "Episode: 116, total numsteps: 2340, episode steps: 20, reward: [146.30885]\n",
      "Episode: 118, total numsteps: 2380, episode steps: 20, reward: [146.30898]\n",
      "Episode: 120, total numsteps: 2420, episode steps: 20, reward: [146.30907]\n",
      "Episode: 122, total numsteps: 2460, episode steps: 20, reward: [146.30916]\n",
      "Episode: 124, total numsteps: 2500, episode steps: 20, reward: [146.30919]\n",
      "Episode: 126, total numsteps: 2540, episode steps: 20, reward: [146.30923]\n",
      "Episode: 128, total numsteps: 2580, episode steps: 20, reward: [146.30927]\n",
      "Episode: 130, total numsteps: 2620, episode steps: 20, reward: [146.3093]\n",
      "Episode: 132, total numsteps: 2660, episode steps: 20, reward: [146.30931]\n",
      "Episode: 134, total numsteps: 2700, episode steps: 20, reward: [146.30931]\n",
      "Episode: 136, total numsteps: 2740, episode steps: 20, reward: [146.30934]\n",
      "Episode: 138, total numsteps: 2780, episode steps: 20, reward: [146.30936]\n",
      "Episode: 140, total numsteps: 2820, episode steps: 20, reward: [146.30937]\n",
      "Episode: 142, total numsteps: 2860, episode steps: 20, reward: [146.30937]\n",
      "Episode: 144, total numsteps: 2900, episode steps: 20, reward: [146.30936]\n",
      "Episode: 146, total numsteps: 2940, episode steps: 20, reward: [146.30937]\n",
      "Episode: 148, total numsteps: 2980, episode steps: 20, reward: [146.30937]\n",
      "Episode: 150, total numsteps: 3020, episode steps: 20, reward: [146.30937]\n",
      "Episode: 152, total numsteps: 3060, episode steps: 20, reward: [146.30937]\n",
      "Episode: 154, total numsteps: 3100, episode steps: 20, reward: [146.30939]\n",
      "Episode: 156, total numsteps: 3140, episode steps: 20, reward: [146.30939]\n",
      "Episode: 158, total numsteps: 3180, episode steps: 20, reward: [146.3094]\n",
      "Episode: 160, total numsteps: 3220, episode steps: 20, reward: [146.3094]\n",
      "Episode: 162, total numsteps: 3260, episode steps: 20, reward: [146.3094]\n",
      "Episode: 164, total numsteps: 3300, episode steps: 20, reward: [146.3094]\n",
      "Episode: 166, total numsteps: 3340, episode steps: 20, reward: [146.30942]\n",
      "Episode: 168, total numsteps: 3380, episode steps: 20, reward: [146.3094]\n",
      "Episode: 170, total numsteps: 3420, episode steps: 20, reward: [146.30942]\n",
      "Episode: 172, total numsteps: 3460, episode steps: 20, reward: [146.3094]\n",
      "Episode: 174, total numsteps: 3500, episode steps: 20, reward: [146.3094]\n",
      "Episode: 176, total numsteps: 3540, episode steps: 20, reward: [146.3094]\n",
      "Episode: 178, total numsteps: 3580, episode steps: 20, reward: [146.30942]\n",
      "Episode: 180, total numsteps: 3620, episode steps: 20, reward: [146.3094]\n",
      "Episode: 182, total numsteps: 3660, episode steps: 20, reward: [146.3094]\n",
      "Episode: 184, total numsteps: 3700, episode steps: 20, reward: [146.30942]\n",
      "Episode: 186, total numsteps: 3740, episode steps: 20, reward: [146.30942]\n",
      "Episode: 188, total numsteps: 3780, episode steps: 20, reward: [146.30942]\n",
      "Episode: 190, total numsteps: 3820, episode steps: 20, reward: [146.3094]\n",
      "Episode: 192, total numsteps: 3860, episode steps: 20, reward: [146.30942]\n",
      "Episode: 194, total numsteps: 3900, episode steps: 20, reward: [146.30942]\n",
      "Episode: 196, total numsteps: 3940, episode steps: 20, reward: [146.30942]\n",
      "Episode: 198, total numsteps: 3980, episode steps: 20, reward: [146.3094]\n",
      "Episode: 200, total numsteps: 4020, episode steps: 20, reward: [146.30942]\n",
      "Episode: 202, total numsteps: 4060, episode steps: 20, reward: [146.30942]\n",
      "Episode: 204, total numsteps: 4100, episode steps: 20, reward: [146.30942]\n",
      "Episode: 206, total numsteps: 4140, episode steps: 20, reward: [146.30942]\n",
      "Episode: 208, total numsteps: 4180, episode steps: 20, reward: [146.30942]\n",
      "Episode: 210, total numsteps: 4220, episode steps: 20, reward: [146.30943]\n",
      "Episode: 212, total numsteps: 4260, episode steps: 20, reward: [146.30942]\n",
      "Episode: 214, total numsteps: 4300, episode steps: 20, reward: [146.30942]\n",
      "Episode: 216, total numsteps: 4340, episode steps: 20, reward: [146.30943]\n",
      "Episode: 218, total numsteps: 4380, episode steps: 20, reward: [146.30942]\n",
      "Episode: 220, total numsteps: 4420, episode steps: 20, reward: [146.30942]\n",
      "Episode: 222, total numsteps: 4460, episode steps: 20, reward: [146.30942]\n",
      "Episode: 224, total numsteps: 4500, episode steps: 20, reward: [146.30943]\n",
      "Episode: 226, total numsteps: 4540, episode steps: 20, reward: [146.30943]\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[16], line 43\u001b[0m\n\u001b[1;32m     39\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(memory) \u001b[38;5;241m>\u001b[39m batch_size:\n\u001b[1;32m     40\u001b[0m     \u001b[38;5;66;03m# Number of updates per step in environment\u001b[39;00m\n\u001b[1;32m     41\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(updates_per_step):\n\u001b[1;32m     42\u001b[0m         \u001b[38;5;66;03m# Update parameters of all the networks\u001b[39;00m\n\u001b[0;32m---> 43\u001b[0m         critic_1_loss, critic_2_loss, policy_loss, ent_loss, alpha \u001b[38;5;241m=\u001b[39m \u001b[43magent\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mupdate_parameters\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmemory\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbatch_size\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mupdates\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     44\u001b[0m         \u001b[38;5;66;03m# critic_1_loss, critic_2_loss, ent_loss, alpha = agent.update_critic_parameters(memory, batch_size, updates)\u001b[39;00m\n\u001b[1;32m     45\u001b[0m         \u001b[38;5;66;03m# policy_loss = agent.update_policy_parameters(memory, batch_size, updates)\u001b[39;00m\n\u001b[1;32m     46\u001b[0m         writer\u001b[38;5;241m.\u001b[39madd_scalar(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mloss/critic_1\u001b[39m\u001b[38;5;124m'\u001b[39m, critic_1_loss, updates)\n",
      "File \u001b[0;32m/scratch/user/jungangc/E2C-NODE/CCS_MSE2CO_9wells/SAC.py:124\u001b[0m, in \u001b[0;36mSAC.update_parameters\u001b[0;34m(self, memory, batch_size, updates)\u001b[0m\n\u001b[1;32m    122\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcritic_optim\u001b[38;5;241m.\u001b[39mzero_grad()\n\u001b[1;32m    123\u001b[0m \u001b[38;5;66;03m# qf_loss.backward()\u001b[39;00m\n\u001b[0;32m--> 124\u001b[0m \u001b[43mqf_loss\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbackward\u001b[49m\u001b[43m(\u001b[49m\u001b[43mretain_graph\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m    125\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcritic_optim\u001b[38;5;241m.\u001b[39mstep()\n\u001b[1;32m    127\u001b[0m pi, log_pi, _ \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mpolicy\u001b[38;5;241m.\u001b[39msample(state_batch)\n",
      "File \u001b[0;32m/scratch/user/jungangc/mylab/lib/python3.8/site-packages/torch/_tensor.py:487\u001b[0m, in \u001b[0;36mTensor.backward\u001b[0;34m(self, gradient, retain_graph, create_graph, inputs)\u001b[0m\n\u001b[1;32m    477\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m has_torch_function_unary(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m    478\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m handle_torch_function(\n\u001b[1;32m    479\u001b[0m         Tensor\u001b[38;5;241m.\u001b[39mbackward,\n\u001b[1;32m    480\u001b[0m         (\u001b[38;5;28mself\u001b[39m,),\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    485\u001b[0m         inputs\u001b[38;5;241m=\u001b[39minputs,\n\u001b[1;32m    486\u001b[0m     )\n\u001b[0;32m--> 487\u001b[0m \u001b[43mtorch\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mautograd\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbackward\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    488\u001b[0m \u001b[43m    \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgradient\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mretain_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcreate_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minputs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43minputs\u001b[49m\n\u001b[1;32m    489\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/scratch/user/jungangc/mylab/lib/python3.8/site-packages/torch/autograd/__init__.py:200\u001b[0m, in \u001b[0;36mbackward\u001b[0;34m(tensors, grad_tensors, retain_graph, create_graph, grad_variables, inputs)\u001b[0m\n\u001b[1;32m    195\u001b[0m     retain_graph \u001b[38;5;241m=\u001b[39m create_graph\n\u001b[1;32m    197\u001b[0m \u001b[38;5;66;03m# The reason we repeat same the comment below is that\u001b[39;00m\n\u001b[1;32m    198\u001b[0m \u001b[38;5;66;03m# some Python versions print out the first line of a multi-line function\u001b[39;00m\n\u001b[1;32m    199\u001b[0m \u001b[38;5;66;03m# calls in the traceback and some print out the last line\u001b[39;00m\n\u001b[0;32m--> 200\u001b[0m \u001b[43mVariable\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_execution_engine\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mrun_backward\u001b[49m\u001b[43m(\u001b[49m\u001b[43m  \u001b[49m\u001b[38;5;66;43;03m# Calls into the C++ engine to run the backward pass\u001b[39;49;00m\n\u001b[1;32m    201\u001b[0m \u001b[43m    \u001b[49m\u001b[43mtensors\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgrad_tensors_\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mretain_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcreate_graph\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minputs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    202\u001b[0m \u001b[43m    \u001b[49m\u001b[43mallow_unreachable\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maccumulate_grad\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# Training Loop\n",
    "# batch_size = 256\n",
    "batch_size = 256\n",
    "total_numsteps = 0\n",
    "updates = 0\n",
    "updates_per_step = 1\n",
    "total_episodes = 500\n",
    "max_steps = 20\n",
    "\n",
    "max_reward = 0\n",
    "# eval = True\n",
    "# for i_episode in itertools.count(1):\n",
    "for i_episode in range(total_episodes):\n",
    "    episode_reward = 0\n",
    "    episode_steps = 0\n",
    "    done = False\n",
    "    state = env.reset(z0)\n",
    "    \n",
    "    # while not done:\n",
    "    for j in range(max_steps):\n",
    "        # if total_numsteps<2000:\n",
    "        if total_numsteps<0:\n",
    "            action = env.sample_action()  # Sample random action    --- exploration\n",
    "        else:\n",
    "            action = agent.select_action(state)  # Sample action from policy   --- exploitation\n",
    "\n",
    "        # print(episode_steps)\n",
    "        next_state, reward, done  = env.step(action) # Step\n",
    "        \n",
    "        # print(next_state)\n",
    "\n",
    "        # Ignore the \"done\" signal if it comes from hitting the time horizon.\n",
    "        # (https://github.com/openai/spinningup/blob/master/spinup/algos/sac/sac.py)\n",
    "        # mask = 1 if episode_steps == num_t_slt else float(not done))\n",
    "\n",
    "        # memory.push(state, action, reward, next_state, mask) # Append transition to memory\n",
    "        memory.push(state, action, reward, next_state)\n",
    "        \n",
    "        if len(memory) > batch_size:\n",
    "            # Number of updates per step in environment\n",
    "            for i in range(updates_per_step):\n",
    "                # Update parameters of all the networks\n",
    "                critic_1_loss, critic_2_loss, policy_loss, ent_loss, alpha = agent.update_parameters(memory, batch_size, updates)\n",
    "                # critic_1_loss, critic_2_loss, ent_loss, alpha = agent.update_critic_parameters(memory, batch_size, updates)\n",
    "                # policy_loss = agent.update_policy_parameters(memory, batch_size, updates)\n",
    "                writer.add_scalar('loss/critic_1', critic_1_loss, updates)\n",
    "                writer.add_scalar('loss/critic_2', critic_2_loss, updates)\n",
    "                writer.add_scalar('loss/policy', policy_loss, updates)\n",
    "                writer.add_scalar('loss/entropy_loss', ent_loss, updates)\n",
    "                writer.add_scalar('entropy_temprature/alpha', alpha, updates)\n",
    "                updates += 1\n",
    "                \n",
    "        state = next_state\n",
    "        episode_steps += 1\n",
    "        total_numsteps += 1\n",
    "        episode_reward += reward\n",
    "\n",
    "    # if total_numsteps > num_t_slt:\n",
    "    #     break\n",
    "\n",
    "    writer.add_scalar('reward/train', episode_reward.detach().cpu().numpy(), i_episode)\n",
    "    if i_episode % 2 == 0:\n",
    "        print(\"Episode: {}, total numsteps: {}, episode steps: {}, reward: {}\".format(i_episode, total_numsteps, episode_steps, episode_reward.detach().cpu().numpy()))\n",
    "\n",
    "    if episode_reward >= max_reward:\n",
    "        agent.save_checkpoint(env_name='ROM_RL')\n",
    "        max_reward = episode_reward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "dcdf7bbc-eedb-496d-a0a4-97f6ee746505",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading models from checkpoints/sac_checkpoint_ROM_RL_\n"
     ]
    }
   ],
   "source": [
    "ckpt_path = \"checkpoints/sac_checkpoint_{}_{}\".format('ROM_RL', \"\")\n",
    "agent.load_checkpoint(ckpt_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "25e1a792-5c08-4705-b727-85c7711732c3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "### choose the best sac for planning\n",
    "dt = torch.tensor(np.ones((1,1)), dtype=torch.float32).to(device)  # dt=20days, normalized to 1\n",
    "state = z0\n",
    "env = Environment(state, num_t_slt, num_prod, num_inj, my_rom)\n",
    "done = False\n",
    "ACTIONS = []\n",
    "Yobs_pred = []\n",
    "while not done:\n",
    "    action = agent.select_action(state, evaluate=True)\n",
    "\n",
    "    _, yobs_pred = my_rom.predict_latent(state, dt, action)\n",
    "    next_state, reward, done = env.step(action) \n",
    "\n",
    "    \n",
    "    state = next_state\n",
    "    ACTIONS.append(action)\n",
    "    Yobs_pred.append(yobs_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "629b21c1-1c03-4ed9-a126-0e46f3062ff7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Q_min = 0\n",
    "# Q_max = 1000*1.20\n",
    "# Q_diff = Q_max - Q_min\n",
    "p_max = 4069.2\n",
    "p_min = 2200\n",
    "p_diff = p_max - p_min\n",
    "\n",
    "Q_min = 0\n",
    "Q_max_w = 3151.0\n",
    "Q_max_g = 1.2e6\n",
    "Q_diff_w = Q_max_w - Q_min\n",
    "Q_diff_g = Q_max_g - Q_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0bdb1a59-d7c8-45f5-8ead-a69241956917",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([1, 4])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ACTIONS[0][:,num_prod:].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "d50a11ef-759a-470a-a327-da85e82014f2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "npv = 0\n",
    "for i_time in range(len(Yobs_pred)):\n",
    "    alpha = np.power((1-0.05),i_time*100/365)\n",
    "    Qpro_w = Yobs_pred[i_time][:,:num_prod]*Q_diff_w+Q_min\n",
    "    Qpro_g = Yobs_pred[i_time][:,num_prod:num_prod*2]*Q_diff_g+Q_min\n",
    "    Qinj_g = ACTIONS[i_time][:, num_prod:]*Q_diff_g+Q_min\n",
    "    lf3toton =0.1167*4.536e-4 # convert lf^3 to ton \n",
    "    pv = np.sum((50-10)*lf3toton*Qinj_g.detach().cpu().numpy().squeeze())- np.sum(5*Qpro_w.detach().cpu().numpy().squeeze()) - np.sum(50*lf3toton*Qpro_g.detach().cpu().numpy().squeeze())\n",
    "    npv += pv*100*alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "11a8187d-014e-48dd-8844-2e1071a26762",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14183097.589750618"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "cd97d782-3639-4776-95cf-20f17cf671df",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Controls_opt = np.zeros((len(ACTIONS), num_well))\n",
    "for i in range(len(ACTIONS)):\n",
    "    Controls_opt[i, :num_prod] = ACTIONS[i][:,:num_prod].detach().cpu().numpy().squeeze()*p_diff+p_min\n",
    "    Controls_opt[i, num_prod:num_well] = ACTIONS[i][:,num_prod:num_well].detach().cpu().numpy().squeeze()*Q_diff_g+Q_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "24fbf1e5-658c-4351-b4aa-bd6b76310f8f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAELCAYAAADJDcLzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA91ElEQVR4nO3de1xUdf4/8NcMw3C/icolATVYQMDLKrJe8kKseFnvbZqmYFuUilK6mvStJF2V1rzsL1utTF3dNHNXXSOVxEutVqaWqKmkqckWYkmCYlzn/ftDmBhgjsBc5PJ6Ph7zOJ75nPm8P+/PHIf3nDlzRiUiAiIiIiKqlfp+D4CIiIioMWOxRERERKSAxRIRERGRAhZLRERERApYLBEREREpYLFEREREpIDFEhEREZECzf0ewP2k0+nwww8/wMXFBSqV6n4Ph4iIiOpARHDr1i34+vpCrbb8cZ8WXSz98MMP8PPzu9/DICIiogbIzs5Gu3btLB6nRRdLLi4uAIDnANjd36EQERFRHRUDWIFf/45bWosulio/erMDiyUiIqKmxlqn0PAEbyIiIiIFLJaIiIiIFLBYIiIiIlLAYomIiIhIAYslIiIiIgWNtliKj4/HqFGj7vcwiIiIqIWrV7EUHx8PlUoFlUoFrVaLwMBALFiwAGVlZZYan1FFRUWIj49HREQENBoNCysiIiKyiHofWRo8eDBycnJw4cIFzJ49GykpKVi6dGmt25aUlJg8QGPKy8vh4OCAmTNnIiYmxmJxiIiIqGWrd7FkZ2cHb29vBAQEYOrUqYiJicGuXbsA/PrR2aJFi+Dr64vg4GAAwOnTpxEdHQ0HBwd4enoiISEBt2/f1vdZXl6OWbNmwd3dHZ6enpg7dy5ERHEcTk5OWL16NZ566il4e3vXNw0iIiKiOjH5nCUHBweDI0j79+9HVlYW9u3bh7S0NBQWFiI2NhYeHh44duwYtm3bhoyMDCQmJuofs2zZMmzYsAHr1q3D4cOHkZeXhx07dpg6tBqKi4tRUFBgcCMiIiJS0uBiSUSQkZGB9PR0REdH6+93cnLC2rVrERYWhrCwMGzevBlFRUXYuHEjwsPDER0djVWrVmHTpk3Izc0FAKxcuRLJyckYM2YMQkNDsWbNGri5uZmeXTVLliyBm5ub/sYf0SUiIqJ7qXexlJaWBmdnZ9jb22PIkCEYN24cUlJS9O0RERHQarX69XPnzqFLly5wcnLS39enTx/odDpkZWUhPz8fOTk5iIqK0rdrNBr06NGjgSkZl5ycjPz8fP0tOzvb7DGIiIioean3D+kOHDgQq1evhlarha+vLzQawy6qFkWNjZ2dHezs+JO5REREVHf1PrLk5OSEwMBA+Pv71yiUahMaGorMzEwUFhbq7zty5AjUajWCg4Ph5uYGHx8fHD16VN9eVlaGEydO1HdoRERERGZn8YtSTpw4Efb29oiLi8OZM2dw8OBBzJgxA5MmTYKXlxcAICkpCampqdi5cyfOnz+PadOm4ebNm/fs++zZszh58iTy8vKQn5+PkydP4uTJk5ZNiIiIiFqUen8MV1+Ojo5IT09HUlISIiMj4ejoiLFjx2L58uX6bWbPno2cnBzExcVBrVbjiSeewOjRo5Gfn6/Y99ChQ/Hdd9/p17t16wYA97zsABEREVFdqaQFVxYFBQVwc3PDPAA8k4mIiKhpKAaQCiA/Px+urq4Wj9dofxuOiIiIqDFgsURERESkgMUSERERkQIWS0REREQKWCwRERERKWCxRERERKSAxRIRERGRAhZLRERERApYLBEREREpYLFEREREpIDFEhEREZECFktEREREClgsERERESlotMVSfHw8Ro0adb+HQURERC1cvYql+Ph4qFQqqFQqaLVaBAYGYsGCBSgrK7PU+Iw6dOgQRo4cCR8fHzg5OaFr16549913rT4OIiIiat7qfWRp8ODByMnJwYULFzB79mykpKRg6dKltW5bUlJi8gCN+fTTT9G5c2f8+9//xqlTpzBlyhRMnjwZaWlpFotJRERELU+9iyU7Ozt4e3sjICAAU6dORUxMDHbt2gXg14/OFi1aBF9fXwQHBwMATp8+jejoaDg4OMDT0xMJCQm4ffu2vs/y8nLMmjUL7u7u8PT0xNy5cyEiiuN44YUXsHDhQvTu3RsPPvggkpKSMHjwYGzfvr2+KREREREZZfI5Sw4ODgZHkPbv34+srCzs27cPaWlpKCwsRGxsLDw8PHDs2DFs27YNGRkZSExM1D9m2bJl2LBhA9atW4fDhw8jLy8PO3bsqPdY8vPz0apVK6PtxcXFKCgoMLgRERERKWlwsSQiyMjIQHp6OqKjo/X3Ozk5Ye3atQgLC0NYWBg2b96MoqIibNy4EeHh4YiOjsaqVauwadMm5ObmAgBWrlyJ5ORkjBkzBqGhoVizZg3c3NzqNZ73338fx44dw5QpU4xus2TJEri5uelvfn5+DUueiIiIWox6F0tpaWlwdnaGvb09hgwZgnHjxiElJUXfHhERAa1Wq18/d+4cunTpAicnJ/19ffr0gU6nQ1ZWFvLz85GTk4OoqCh9u0ajQY8ePeo8poMHD2LKlCl4++23ERYWZnS75ORk5Ofn62/Z2dl1jkFEREQtk6a+Dxg4cCBWr14NrVYLX19faDSGXVQtiqzh448/xvDhw7FixQpMnjxZcVs7OzvY2dlZaWRERETUHNT7yJKTkxMCAwPh7+9fo1CqTWhoKDIzM1FYWKi/78iRI1Cr1QgODoabmxt8fHxw9OhRfXtZWRlOnDhxz74PHTqEYcOG4dVXX0VCQkJ9UyEiIiK6J4tflHLixImwt7dHXFwczpw5g4MHD2LGjBmYNGkSvLy8AABJSUlITU3Fzp07cf78eUybNg03b95U7PfgwYMYNmwYZs6cibFjx+LatWu4du0a8vLyLJ0SERERtSAWL5YcHR2Rnp6OvLw8REZG4pFHHsHDDz+MVatW6beZPXs2Jk2ahLi4OPTq1QsuLi4YPXq0Yr//+Mc/cOfOHSxZsgQ+Pj7625gxYyydEhEREbUgKrnXBY2asYKCAri5uWEeAJ7JRERE1DQUA0jF3UsGubq6Wjxeo/1tOCIiIqLGgMUSERERkQIWS0REREQKWCwRERERKWCxRERERKSAxRIRERGRAhZLRERERApYLBEREREpYLFEREREpIDFEhEREZECFktEREREClgsERERESlgsURERESkoNEWS/Hx8Rg1atT9HgYRERG1cPUqluLj46FSqaBSqaDVahEYGIgFCxagrKzMUuMzKisrCwMHDoSXlxfs7e3RsWNHvPjiiygtLbX6WIiIiKj50tT3AYMHD8b69etRXFyM3bt3Y/r06bC1tUVycnKNbUtKSqDVas0y0OpsbW0xefJk/Pa3v4W7uzsyMzPx1FNPQafTYfHixRaJSURERC1PvT+Gs7Ozg7e3NwICAjB16lTExMRg165dAH796GzRokXw9fVFcHAwAOD06dOIjo6Gg4MDPD09kZCQgNu3b+v7LC8vx6xZs+Du7g5PT0/MnTsXIqI4jo4dO2LKlCno0qULAgICMGLECEycOBH//e9/jT6muLgYBQUFBjciIiIiJSafs+Tg4ICSkhL9+v79+5GVlYV9+/YhLS0NhYWFiI2NhYeHB44dO4Zt27YhIyMDiYmJ+scsW7YMGzZswLp163D48GHk5eVhx44d9RrHxYsXsXfvXvTv39/oNkuWLIGbm5v+5ufnV/+EiYiIqEVpcLEkIsjIyEB6ejqio6P19zs5OWHt2rUICwtDWFgYNm/ejKKiImzcuBHh4eGIjo7GqlWrsGnTJuTm5gIAVq5cieTkZIwZMwahoaFYs2YN3Nzc6jSO3r17w97eHkFBQXjooYewYMECo9smJycjPz9ff8vOzm5o+kRERNRC1PucpbS0NDg7O6O0tBQ6nQ4TJkxASkqKvj0iIsLgPKVz586hS5cucHJy0t/Xp08f6HQ6ZGVlwd7eHjk5OYiKivp1UBoNevTocc+P4gBg69atuHXrFjIzMzFnzhy89tprmDt3bq3b2tnZwc7Orr4pExERUQtW72Jp4MCBWL16NbRaLXx9faHRGHZRtSiyhsqP0jp16oTy8nIkJCRg9uzZsLGxseo4iIiIqHmq98dwTk5OCAwMhL+/f41CqTahoaHIzMxEYWGh/r4jR45ArVYjODgYbm5u8PHxwdGjR/XtZWVlOHHiRH2HBp1Opz/iRURERGQO9T6yVF8TJ07E/PnzERcXh5SUFPz444+YMWMGJk2aBC8vLwBAUlISUlNTERQUhJCQECxfvhw3b95U7Pfdd9+Fra0tIiIiYGdnh+PHjyM5ORnjxo2Dra2tpdMiIiKiFsLixZKjoyPS09ORlJSEyMhIODo6YuzYsVi+fLl+m9mzZyMnJwdxcXFQq9V44oknMHr0aOTn5xsfuEaDV199Fd988w1EBAEBAUhMTMRzzz1n6ZSIiIioBVFJXc6ibqYKCgrg5uaGeQB42jcREVHTUAwgFUB+fj5cXV0tHq/R/jYcERERUWPAYomIiIhIAYslIiIiIgUsloiIiIgUsFgiIiIiUsBiiYiIiEgBiyUiIiIiBSyWiIiIiBSwWCIiIiJSwGKJiIiISAGLJSIiIiIFLJaIiIiIFDTaYik+Ph6jRo2638MgIiKiFq5exVJ8fDxUKhVUKhW0Wi0CAwOxYMEClJWVWWp8dXLx4kW4uLjA3d39vo6DiIiImp96H1kaPHgwcnJycOHCBcyePRspKSlYunRprduWlJSYPMB7KS0txWOPPYaHHnrI4rGIiIio5al3sWRnZwdvb28EBARg6tSpiImJwa5duwD8+tHZokWL4Ovri+DgYADA6dOnER0dDQcHB3h6eiIhIQG3b9/W91leXo5Zs2bB3d0dnp6emDt3LkSkTuN58cUXERISgkcffbS+qRARERHdk8nnLDk4OBgcQdq/fz+ysrKwb98+pKWlobCwELGxsfDw8MCxY8ewbds2ZGRkIDExUf+YZcuWYcOGDVi3bh0OHz6MvLw87Nix456xDxw4gG3btuGNN96o01iLi4tRUFBgcCMiIiJS0uBiSUSQkZGB9PR0REdH6+93cnLC2rVrERYWhrCwMGzevBlFRUXYuHEjwsPDER0djVWrVmHTpk3Izc0FAKxcuRLJyckYM2YMQkNDsWbNGri5uSnGv3HjBuLj47Fhwwa4urrWacxLliyBm5ub/ubn59fQ9ImIiKiFqHexlJaWBmdnZ9jb22PIkCEYN24cUlJS9O0RERHQarX69XPnzqFLly5wcnLS39enTx/odDpkZWUhPz8fOTk5iIqK0rdrNBr06NFDcRxPPfUUJkyYgH79+tV57MnJycjPz9ffsrOz6/xYIiIiapk09X3AwIEDsXr1ami1Wvj6+kKjMeyialFkSQcOHMCuXbvw2muvAbh7pEun00Gj0eCtt97CE088UeMxdnZ2sLOzs8r4iIiIqHmo95ElJycnBAYGwt/fv0ahVJvQ0FBkZmaisLBQf9+RI0egVqsRHBwMNzc3+Pj44OjRo/r2srIynDhxQrHfzz77DCdPntTfFixYABcXF5w8eRKjR4+ub1pEREREtbL4RSknTpwIe3t7xMXF4cyZMzh48CBmzJiBSZMmwcvLCwCQlJSE1NRU7Ny5E+fPn8e0adNw8+ZNxX5DQ0MRHh6uvz3wwANQq9UIDw+Hh4eHpdMiIiKiFsLixZKjoyPS09ORl5eHyMhIPPLII3j44YexatUq/TazZ8/GpEmTEBcXh169esHFxYVHh4iIiKhRUEldL2jUDBUUFMDNzQ3zAPBMJiIioqahGEAqgPz8/Dp/I94Ujfa34YiIiIgaAxZLRERERApYLBEREREpYLFEREREpIDFEhEREZECFktEREREClgsERERESlgsURERESkgMUSERERkQIWS0REREQKWCwRERERKWCxRERERKSAxRIRERGRgkZbLMXHx2PUqFH3exhERETUwtWrWIqPj4dKpYJKpYJWq0VgYCAWLFiAsrIyS43PqCtXrujHUvX2+eefW30sRERE1Hxp6vuAwYMHY/369SguLsbu3bsxffp02NraIjk5uca2JSUl0Gq1ZhmoMRkZGQgLC9Ove3p6WjQeERERtSz1/hjOzs4O3t7eCAgIwNSpUxETE4Ndu3YB+PWjs0WLFsHX1xfBwcEAgNOnTyM6OhoODg7w9PREQkICbt++re+zvLwcs2bNgru7Ozw9PTF37lyISJ3G4+npCW9vb/3N1ta2vikRERERGWXyOUsODg4oKSnRr+/fvx9ZWVnYt28f0tLSUFhYiNjYWHh4eODYsWPYtm0bMjIykJiYqH/MsmXLsGHDBqxbtw6HDx9GXl4eduzYUaf4I0aMQNu2bdG3b1990WZMcXExCgoKDG5EREREShpcLIkIMjIykJ6ejujoaP39Tk5OWLt2LcLCwhAWFobNmzejqKgIGzduRHh4OKKjo7Fq1Sps2rQJubm5AICVK1ciOTkZY8aMQWhoKNasWQM3NzfF+M7Ozli2bBm2bduGDz/8EH379sWoUaMUC6YlS5bAzc1Nf/Pz82to+kRERNRC1PucpbS0NDg7O6O0tBQ6nQ4TJkxASkqKvj0iIsLgPKVz586hS5cucHJy0t/Xp08f6HQ6ZGVlwd7eHjk5OYiKivp1UBoNevToofhRXOvWrTFr1iz9emRkJH744QcsXboUI0aMqPUxycnJBo8pKChgwURERESK6l0sDRw4EKtXr4ZWq4Wvry80GsMuqhZF1hYVFYV9+/YZbbezs4OdnZ0VR0RERERNXb0/hnNyckJgYCD8/f1rFEq1CQ0NRWZmJgoLC/X3HTlyBGq1GsHBwXBzc4OPjw+OHj2qby8rK8OJEyfqOzScPHkSPj4+9X4cERERkTH1PrJUXxMnTsT8+fMRFxeHlJQU/Pjjj5gxYwYmTZoELy8vAEBSUhJSU1MRFBSEkJAQLF++HDdv3lTs9x//+Ae0Wi26desGANi+fTvWrVuHtWvXWjolIiIiakEsXiw5OjoiPT0dSUlJiIyMhKOjI8aOHYvly5frt5k9ezZycnIQFxcHtVqNJ554AqNHj0Z+fr5i3wsXLsR3330HjUaDkJAQbN26FY888oilUyIiIqIWRCV1vaBRM1RQUAA3NzfMA8AzmYiIiJqGYgCpAPLz8+Hq6mrxeI32t+GIiIiIGgMWS0REREQKWCwRERERKWCxRERERKTA4t+Ga8wqz20vvs/jICIiorqr/Lttre+otehi6datWwCAFfd5HERERFR/N27cuOdvyZpDi750gE6nww8//AAXFxeoVCqLxKj8/bns7GyLfb3R0jGaQw6M0Xj6Z4zGFaM55MAYjad/a8XIz8+Hv78/fv75Z7i7u1skRlUt+siSWq1Gu3btrBLL1dXV4teCsHSM5pADYzSe/hmjccVoDjkwRuPp31ox1GrrnHrNE7yJiIiIFLBYIiIiIlLAYsnC7OzsMH/+fNjZWe4HVSwdoznkwBiNp3/GaFwxmkMOjNF4+m9OMapq0Sd4ExEREd0LjywRERERKWCxRERERKSAxRIRERGRAhZLRERERApYLBEREREpYLFEREREpKBF/9xJc1VWVoavv/4a165dAwB4e3ujU6dOsLW1vc8jq7vmkAMAXLt2DUePHjXIIyoqCt7e3vd5ZHVnjRyawzwR0f1XXFwMAGa//hKLJSuz5B8FnU6Hl19+GW+88Qby8/MN2tzc3JCYmIhXXnnF5N/SsWQhY60cAMvmUVhYiKeffhrvvfceVCoVWrVqBQDIy8uDiOCxxx7Dm2++CUdHR5NjWSoPa+RgzXmqlJ+fbzBX5vzFcmsV+ZbMgTHqrrk839bIw5Ix9u3bhxUrVuCzzz5DQUEBgLu/S9erVy/MmjULMTExJseAkFXcvn1bJk6cKDY2NqLRaKRt27bStm1b0Wg0YmNjI48//rgUFhaaFGPOnDnSpk0bWbNmjVy+fFnu3Lkjd+7ckcuXL8ubb74pbdu2lblz5za4//Lycvm///s/cXd3F5VKZXBzd3eXF198UcrLyxt1DtbK409/+pMEBQXJ3r17paysTH9/WVmZpKeny29+8xt58sknG3Ue1sjBGjEqvf322xIaGipqtdrgFhoaKmvXrjWpb2vsU5bOgTHqrrk839bIw9IxNmzYIBqNRsaPHy/r16+X3bt3y+7du2X9+vXy2GOPia2trWzcuNGkHEREWCxZiTX+KHh5ecnevXuNtu/du1fatm3b4P6tUchYOgcR6+Th7u4uR44cMdp++PBhcXd3NymGpfOwRg7WiCEi8te//lUcHR1l3rx5cvDgQTl79qycPXtWDh48KMnJyeLk5CRLly5tcP/W2KcsnQNj1F1zeb6tkYelYwQFBcmqVauMtr/xxhsSGBjY4P4rsViyEmv8UXB0dJRTp04Zbc/MzBQnJ6cG92+NQsbSOYhYJw9XV1c5duyY0fYvvvhCXF1dTYph6TyskYM1YoiI+Pv7y9atW422v/fee+Ln59fg/q2xT1k6B8aou+byfFsjD0vHsLOzk/PnzxttP3/+vNjb2ze4/0r8NpyV6HQ6aLVao+1arRY6nc6kGAMGDMCf//xn/PTTTzXafvrpJzz//PMYMGBAg/u/desWfH19jbb7+PigsLCwwf0Dls8BsE4ef/jDH5CQkICvvvqqRttXX32FqVOnYvjw4SbFsHQe1sjBGjEA4Pr164iIiDDaHhERUes+V1fW2KcsnQNj1F1zeb6tkYelY4SFheGdd94x2r5u3Tp06tSpwf3rmVxuUZ1MmDBBunXrJl9++WWNti+//FK6d+8uEydONCnG1atXJTw8XDQajXTr1k0GDx4sgwcPlm7duolGo5HOnTvL1atXG9z/0KFDZdCgQfLjjz/WaPvxxx9l8ODBMmzYMFNSsHgOItbJIy8vTwYPHiwqlUpatWolISEhEhISIq1atRK1Wi1DhgyRn3/+2aQYls7DGjlYI4aIyEMPPSSTJ0+W0tLSGm1lZWUyefJk6devX4P7t8Y+ZekcGKPumsvzbY08LB3j4MGD4uTkJBEREfLcc89JamqqpKamynPPPSedO3cWZ2dn+fjjj01JQUREVCIippdcdC8///wzJkyYgPT0dHh4eKBt27YA7r57uHnzJmJjY7F582a4u7ubFEen0yE9PR2ff/65wbcOevXqhUGDBpn0LbLs7GwMHToU58+fR0REBLy8vAAAubm5OH36NDp16oS0tDT4+fk12hysmQcAnDt3rtY8QkJCTO7bWnlYMgdrxTh16hRiY2NRWlqKfv36GczVJ598Aq1Wi48++gjh4eEN6t8az4Wlc2CMumsuz7c18rBGjCtXrmD16tW1voY888wzaN++fYP7rsRiycqs8YfHkixdyFgL82h5bt26hX/+85+1ztWECRPg6upqUv/WeC4snQNj1F1zeb6tkUdzeJ1isdQMffHFF/jss88MdsrevXsjMjLyPo+s7ppDDiUlJdi5c2eteYwcOVLxHLbGwho5NId5IqL7q/p1nHx8fBAaGmq2a0WxWLIiS/9RuH79OsaOHYsjR47A39/f4HDn1atX0adPH/z73//WfwTYUJYsZKyVA2DZPC5evIjY2Fj88MMPiIqKMsjj6NGjaNeuHfbs2YPAwECTY1kqD2vkYM15AmpeFNbHxwc9e/Y025XCrVHkWzoHxqi75vJ8WyMPS8Ww2oWMTT7rierkwoUL0rFjR7G3t5f+/fvLo48+Ko8++qj0799f7O3tJTAwUC5cuGBSjLFjx0qvXr1q/Rrl+fPnpXfv3vLII480uP/c3Fzp27evqFQqCQgIkJ49e0rPnj0lICBAVCqV9O3bV3Jzc01JweI5iFgnj5iYGBk5cqTk5+fXaMvPz5eRI0fKoEGDTIph6TyskYM1YohY/qKw1tinrHFhW8aom+byfFsjD0vHsMa1okR4nSWrscYfBWdn51q/bVfp+PHj4uzs3OD+rVHIWDoHEevk4eDgIKdPnzbafurUKXFwcDAphqXzsEYO1oghYvmLwlpjn2ouV1RvDjGay/NtjTwsHcMa14oSYbFkNdb4o+Dp6SmHDh0y2n7w4EHx9PRscP/WKGQsnYOIdfLw8fGRDz74wGj7rl27xMfHx6QYls7DGjlYI4aI5S8Ka419qrlcUb05xGguz7c18rB0DGtcyFiEF6W0Gnd3d1y5csVo+5UrV0y+bMC4ceMQFxeHHTt26H9MEAAKCgqwY8cOTJkyBY899liD+7ezszPot7pbt26Z/EvPls4BsE4eTz75JCZPnowVK1bg1KlTyM3NRW5uLk6dOoUVK1YgPj4eCQkJJsWwdB7WyMEaMQDLXxTWGvuUNS5syxh101yeb2vkYekY1riQMQCes2QtL730knh4eMjy5cslMzNTrl27JteuXZPMzExZvny5tGrVSubPn29SjKKiInnmmWdEq9WKWq0We3t7sbe3F7VaLVqtVqZOnSpFRUUN7n/atGkSEBAg27dvN/g4MT8/X7Zv3y7t27eXxMREi+SgUqnMkoO18hARSU1NFR8fH1GpVPofwFSpVOLj4yOvvvqqyf1bIw9L52CtGJa+KKw1ngtrXNiWMeqmuTzf1sjD0jGscSFjEX4MZ1XW+KMgcncnPHDggGzevFk2b94sBw4cqPVcqfqydDFWVX5+vuzfv1+fw/79+82Sg4h18xARuXTpknz66afy6aefyqVLl8zWrzXzsFQO1oph6SuFW+O5aC5XVG8OMZrL822NPKwRo7y8XHbv3i0vv/yyJCQkSEJCgrz88suyZ88eKS8vN6nvSrx0wH1w+fJlg69PdujQ4T6PqH4KCgpw4sQJgxy6d+9ulgukGaPVapGZmYnQ0FCz9Xk/8rCEgoICHD9+HLm5uQCaXh45OTlYvXo1Dh8+jJycHKjVanTs2BGjRo1CfHw8bGxszBbL0heFtcY+1RyuqA4A58+fr/FV8qYWwxrPtzXmyRp5NPXXKRZLjUR2djbmz5+PdevWmdTPL7/8ghMnTqBVq1Y1fjywqKgI77//PiZPntzg/itfRCv/s54/fx5/+9vfUFxcjMcffxzR0dEmjX/WrFm13v+3v/0Njz/+ODw9PQEAy5cvNylOVYWFhXj//fdx8eJF+Pr6Yvz48fo4DfXll1/Cw8NDXwhv2rQJa9aswdWrVxEQEIDExESMHz/epBgzZszAo48+ioceesikfpSsWrUKX3zxBYYOHYrx48dj06ZNWLJkCXQ6HcaMGYMFCxZAo9E0uP/jx48jJiYGgYGBcHBwwGeffYYJEyagpKQE6enp6NSpE/bu3QsXFxczZkVEzY3FrxVlluNTZLKTJ0+KWq02qY+srCz9tSvUarX069dPvv/+e337tWvXTIqxZ88e0Wq10qpVK7G3t5c9e/ZImzZtJCYmRqKjo8XGxkb2799vUg4qlUq6du0qAwYMMLipVCqJjIyUAQMGyMCBA02KERoaKjdu3BCRu593t2/fXtzc3CQyMlJatWolbdu2NfljoM6dO8u+fftEROTtt98WBwcHmTlzpqxevVqeffZZcXZ2lnfeecekGJXPc1BQkKSmpkpOTo5J/VW3cOFCcXFxkbFjx4q3t7ekpqaKp6en/OUvf5HFixdLmzZt5OWXXzYpRp8+fSQlJUW/vmnTJomKihKRux9DdO3aVWbOnGlSjErFxcWydetWefbZZ2X8+PEyfvx4efbZZ+X999+X4uJis8Qw5tq1a/LKK6+Ypa/s7Gy5detWjftLSkrM8oOhP/30kxw4cED/f+THH3+U1NRUeeWVV+Ts2bMm929Mhw4d5JtvvrFI3zqdTg4cOCBvvfWWfPDBB1JSUmJSf9nZ2QY/DPvJJ5/IhAkTpG/fvjJx4kT59NNPTR2yvPbaa3LlyhWT+7mXDz74QF566SU5fPiwiIjs379fhgwZIrGxsfLmm2+aJcadO3fknXfekSlTpsjgwYNl6NChkpiYKBkZGSb3bY1rRYnwnCWr+c9//qN4W7FihcnF0qhRo2TYsGHy448/yoULF2TYsGHSoUMH+e6770TE9GKpV69e8n//938iIrJlyxbx8PCQF154Qd8+b948+f3vf29SDkuWLJEOHTrUKLo0Go18/fXXJvVdSaVS6f/zTJw4UXr37i03b94UEZFbt25JTEyMPPbYYybFcHBw0L/QdevWTd566y2D9nfffVc6depkUgyVSiUZGRmSlJQkrVu3FltbWxkxYoR88MEHZvmc/sEHH5R///vfInK3mLexsZF//vOf+vbt27dLYGCgSTEcHBzk22+/1a+Xl5eLra2tXLt2TUREPvroI/H19TUphoh1LgqrxBxvhn744QeJjIwUtVotNjY2MmnSJIOiydT/3yIiR48eFTc3N1GpVOLh4SHHjx+XDh06SFBQkDz44IPi4OAgJ06cMCnG3/72t1pvNjY2kpycrF83xZAhQ/T/p2/cuCFRUVGiUqmkTZs2olarJSQkRK5fv97g/nv27Km/5MXOnTtFrVbLiBEj5Pnnn5fRo0eLra2t4iUx6kKlUomNjY3ExMTIe++9Z5GCfs2aNaLRaKR79+7i6uoqmzZtEhcXF3nyySfl6aefFgcHB1m5cqVJMS5cuCABAQHStm1b8fPzE5VKJcOGDZOoqCixsbGRP/7xj1JaWtrg/q1xrSgRFktWU3kUQKVSGb2Z+kLXtm1bg+tN6HQ6eeaZZ8Tf31++/fZbk19MXV1d9X9QysvLRaPRGHxT4/Tp0+Ll5dXwBCp88cUX8pvf/EZmz56tfwdoqWKpY8eO8tFHHxm0HzlyRPz8/EyK4enpKcePHxeRu8/LyZMnDdovXrxo8nW1quZRUlIiW7duldjYWLGxsRFfX1954YUXTCoAHBwc9IW2iIitra2cOXNGv37lyhVxdHRseAIiEhAQoH9HK3K3IFCpVHLnzh0REbl8+bLY29ubFEPE8heFzczMVLxt3brV5P/fkydPlqioKDl27Jjs27dPunfvLj169JC8vDwRuVssqVQqk2LExMTIk08+KQUFBbJ06VJp166dwYUPp0yZIqNGjTIphkqlknbt2kn79u0NbiqVSh544AFp3769dOjQweQYlf83pk6dKp06ddIfLc7Ozpbu3bvLM8880+D+nZyc9P1FRUVJamqqQfvrr78u3bp1a3D/IndzWL9+vYwcOVJsbW3F09NTkpKSFK/XV1+dOnXSv5E7cOCA2NvbyxtvvKFvX79+vYSGhpoUY8iQIfL000+LTqcTkbtfdBoyZIiIiHzzzTfSvn17k74Jbo1rRYmwWLIaX19f2blzp9H2r776yuQXUxcXl1oPk0+fPl3atWsnn3zyicnF0sWLF/Xrzs7OBkcFrly5YpY/bCJ3j/BMnjxZOnfuLKdPnxZbW1uzFkuV7yp9fX1rvPiYI4/HH39c/vSnP4mIyB//+Ed58cUXDdoXL14sERERJsWo+gehqu+++07mz58vAQEBJj3fHTp0kD179ojI3Rc1tVot77//vr79ww8/lPbt2ze4fxGRpKQkCQ8Plz179siBAwdk4MCBMmDAAH373r175cEHHzQphojlLwqr9Gao8n5T/3/7+vrK0aNH9etFRUUyfPhw6dq1q9y4ccMsR5Y8PDz0ryElJSWiVqsNYp44cUIeeOABk2I8/fTT0rVr1xqvVZZ6QxQcHCz/+c9/DNozMjJMKsjc3NwkMzNTRO6+Gar8d6WLFy+a/Eaiag65ubny6quvSkhIiKjVaomMjJS33npLCgoKTIpR2xuiqv9PLl++bHIejo6OBh+vFhcXi62trfz0008icvfInCmvI9a4kLEIiyWrGT58uLz00ktG20+ePGnyu8LIyEjZuHFjrW3Tp08Xd3d3k15MO3furP/jKXL3SFLVw6effPKJye8Iq9uyZYt4eXmJWq026wtpRESEdOvWTZydneVf//qXQfvHH39s8h+E77//Xtq3by/9+vWTWbNmiYODg/Tt21eeeuop6devn2i1Wvnwww9NimGsWKqk0+lqHDWrjxdffFHatGkjTz75pHTo0EHmzZsn/v7+snr1almzZo34+fnJc8891+D+Re4WxY8++qhoNBpRqVTSu3dvg/PF0tPTDQq0hrL0lcI9PT3lnXfekStXrtR6+/DDD00uZJycnGqc01NaWiqjRo2Szp07y6lTp8wS4/Lly/r16m+IvvvuO7O8Idq+fbv4+fnJ66+/rr/P3MVS5Ruitm3bGhwRFbn7hsjOzq7B/Y8YMULmzZsnIiKxsbE1PjZ8++23JSgoqMH9ixj///3JJ59IXFycODk5mXxl6so30SJ3X7NUKpXB69KhQ4ekXbt2JsXw9fU1+Oj2559/FpVKpS/0Ll26ZNJzYa3r5rFYspJPPvnEoNCo7vbt24rVcV0sXrxYf3izNlOnTjWpIFu9erWkpaUZbU9OTtYfTTGn7Oxs2blzp9y+fdss/aWkpBjcqv+u0J///GcZP368yXF+/vlnef7556VTp05ib28vWq1WAgICZMKECXLs2DGT+2/fvr3+3ZkllJeXy6JFi+QPf/iDLF68WHQ6nWzZskX8/PzE09NT4uPjzfac/PLLL7WetGwulr4o7KBBg2ThwoVG283xZigiIqJGYS/ya8Hk7+9vcrEUEhJicL5gWlqa/iNREZHPP//c5D+elf73v/9JdHS0DB48WHJycsxeLA0dOlRGjx4tHh4eNQrlzz//3KRTBs6ePSuenp4yefJkWbhwoTg7O8vjjz8uixYtksmTJ4udnZ2sX7/epBzUarXim6H8/Pwa50LW1/Tp0yUoKEj+8pe/SM+ePSUuLk5CQkJkz549snfvXomIiJAnnnjCpBhxcXHSv39/OXfunFy6dEnGjRtn8BHloUOHTDrtwVrXm2OxREQtgiUvCrt9+3bZtGmT0fa8vDzZsGGDSTHmzp1r9Lyq0tJSGTFihMkFWUpKimzZssVo+wsvvCBjxowxKUZVOp1OFi9eLN7e3mJjY2O2Yik+Pt7gtnXrVoP2OXPmSGxsrEkxLl68KOPHjxcXFxf9R662trbSu3dv2bFjh0l9i9z7yLE53L59W5566ikJDw+XhIQEKS4ulqVLl4pWqxWVSiUDBgwweQy5ubnyu9/9Tv//LiAgwOAco23btsn/+3//z9RULHYx5kq8zhIRtShN9aKwZWVluHPnjtGL+JWVleH7779HQECAxcZw584d2NjYmPx7YdWdOHEChw8fxuTJk+Hh4WHWvmtTWFgIGxsb2Nvbm9yXiOD69evQ6XRo3bo1bG1tzTDC+6uoqAilpaVmvb7ZhQsXUFxcjJCQEJOuzXa/8Id0iahF6dChA3r16oVevXrpC6Xs7Gw88cQTFotpjv41Go3i1Y5zcnLwyiuvmBTjXm7cuIGpU6eavd/u3bsjKSkJHh4eFn8uACAvLw/Tpk0zS18qlQpeXl7w8fHRF0rWyMGSMezt7eHi4mLWGEFBQQgPD69RKJkjxi+//ILDhw/j7NmzNdqKioqwceNGk/oHeAVvIiJkZmbit7/9LcrLy5tk/4zRuGI0hxyaSoxvvvkGgwYNwtWrV6FSqdC3b19s2bIFvr6+AIDc3Fz4+vqanEPTOxZGRFRPu3btUmy/dOlSo+6fMRpXjOaQQ3OJ8fzzzyM8PBzHjx/HzZs38eyzz6Jv3744dOgQ/P39Teq7Kh5ZIqJmT61WQ6VSQenlTqVSNfjdp6X7Z4zGFaM55NBcYnh5eSEjIwMREREA7p5DNm3aNOzevRsHDx6Ek5OTWY4s8ZwlImr2fHx8sH37duh0ulpvX375ZaPunzEaV4zmkENzifHLL78YnAelUqmwevVqDB8+HP3798c333xjagoAWCwRUQvQvXt3nDhxwmj7vd753u/+GaNxxWgOOTSXGCEhITh+/HiN+1etWoWRI0dixIgRDe67Kp6zRETN3pw5c1BYWGi0PTAwEAcPHmy0/TNG44rRHHJoLjFGjx6NLVu2YNKkSTXaVq1aBZ1OhzVr1jS4/0o8Z4mIiIhIAT+GIyIiIlLAYomIiIhIAYslIiIiIgUsloiIiIgUsFgiIiIiUsBiiYiIiEgBiyUiIiIiBS3+opRFRUUoKSm538MgIiKietBqtbC3t7dKrBZdLBUVFcHBoQ2A2/d7KERERFQP3t7euHz5slUKphZdLN09onQbwHMAnCruta1Yaowsq7Yba6tPHw1tv9djKtfvweYeIdW1tFe22Rh5jLnbbausG2urfv+92uvaT9Ux1TWGufLXbyeApuIXsyuWKpuyu6sV6zb6pa5iWWZ4v7piidqXmoqlukZ7GTTQ1foYdZVtqvZhU217dY32uo/B2Db1H0Pdtqtsry32r9uWVbu/5pzdjdHQMZQrzFldx3CPua34BXabsoo45VKxfvcGAKq7oVARquayvNqy+v11Xa9c6hRiNKQvpTHXtl7fvuo7D0r9mbOvWh5XWrGseLoN1svKDR9SaqQrY/fXt73qdqbEKAaw4to1lJSUsFiyHjsAlZNtrOiobd1YW336aOh6XR9jhKqOS3Uty+r3VS8eqi/rWgzVpd2Ux5qyXW21cX0ea1Js48WSyrZiXVN9ebddXbGuVhv+sfx1WWawXtsf/nv/Qb5XMdTwdvPFqFza3GM7tb7deB+qim3UBo/5danS91F1/dclKtor18sN7rcBYAOp2AY12gzHWVHkGF3qqsWqWJZXLMsqxlpbsWTsD/S9/oCb8jhz9lWf+22qtFWf7Mr7a3strK3P6q+d1V9TK7evjaraevUfIxMj9xsbW8V2pRVFVcXTjVLVr+vVC5Taihpzrle9v/o2GiOPqfo0VV23Jp7gTURERKSAxRIRERGRAhZLRERERApYLBEREREpYLFEREREpIDFEhEREZECFktEREREClgsERERESlgsURERESkgMUSERERkQIWS0REREQKWCwRERERKWCxRERERKSAxRIRERGRAhZLRERERApYLBEREREpYLFEREREpIDFEhEREZECFktEREREClgsERERESlgsURERESkQHO/B9A4FOPXqSivWFaulymsl1Zrs6221BhZmqP9Xo+pXDdC6rjU1bKs/LeqYllebb1yCSPr1cdgbEy6au1VY9sY2ba8Wnv1+6s/nTbV2qvfX3Waq7dVLo09XbX1UZ92/XYCaCoGWLEUm7uJSOW6fqmrWN5t11Xcr1JXLFH7UqefgOrLMqhqTGJFLP2yrNr9OoN2XUW7VItVfVnZXl6xtEE5bKr8u+q2letq/XpZtcfqFLfTGNmusr222L9uW1bt/urLyhgNHUO5fluN0THUdSkV/UjFGO4ubcorlmUV7RVPrU3Z3RsAqCqf9jIjy5q7SsPWK5c6hRgN6UtpzLWt17ev+s6DUn/G+rrXsrbX51qWpRWvkxVPt8F69eGVVlu/1/31ba+6nSkximFdLbpYEhE4Ozvj9u0V93so90/1/3zUJFTWiNVfgIiIWgpnZ2eIVH/XbRktulhSqVS4ffs2srOz4erqer+H06QVFBTAz8+Pc2kizqP5cC7Nh3NpHpxH86mcS5XK2EcX5tWii6VKrq6u3HHNhHNpHpxH8+Fcmg/n0jw4j00PT/AmIiIiUsBiiYiIiEhBiy6W7OzsMH/+fNjZ2d3voTR5nEvz4DyaD+fSfDiX5sF5NB9rz6VKrHUqOREREVET1KKPLBERERHdC4slIiIiIgUsloiIiIgUsFgiIiIiUtCki6UlS5YgMjISLi4uaNu2LUaNGoWsrCx9e15eHmbMmIHg4GA4ODjA398fM2fORH5+vkE/V69exbBhw+Do6Ii2bdtizpw5KCsrM9jm0KFD+O1vfws7OzsEBgZiw4YN1kjRasw1lyqVqsbtvffeM9impc8lADz99NN48MEH4eDggDZt2mDkyJE4f/68wTYtfb801zxyn6zbXFYSEQwZMgQqlQo7d+40aGvp+yRgvrls6ftlXeZxwIABNebomWeeMdjGavukNGGxsbGyfv16OXPmjJw8eVKGDh0q/v7+cvv2bREROX36tIwZM0Z27dolFy9elP3790tQUJCMHTtW30dZWZmEh4dLTEyMfPXVV7J7925p3bq1JCcn67e5dOmSODo6yqxZs+Ts2bPy+uuvi42Njezdu9fqOVuKOeZSRASArF+/XnJycvS3X375Rd/OubzrzTfflI8//lguX74sJ06ckOHDh4ufn5+UlZWJCPdLEfPMowj3SZG6zWWl5cuXy5AhQwSA7NixQ38/98m7zDGXItwv6zKP/fv3l6eeespgjvLz8/Xt1twnm3SxVN3169cFgHz88cdGt3n//fdFq9VKaWmpiIjs3r1b1Gq1XLt2Tb/N6tWrxdXVVYqLi0VEZO7cuRIWFmbQz7hx4yQ2NtYCWTQODZlLEan1RaEqzmXtMjMzBYBcvHhRRLhf1qYh8yjCfbI2xubyq6++kgceeEBycnJqzBv3ydo1ZC5FuF9WV9s89u/fX5KSkow+xpr7ZJP+GK66yo+EWrVqpbiNq6srNJq7P4v32WefISIiAl5eXvptYmNjUVBQgK+//lq/TUxMjEE/sbGx+Oyzz8ydQqPRkLmsNH36dLRu3Ro9e/bEunXrDH4VmnNZU2FhIdavX48OHTrAz88PAPfL2jRkHitxnzRU21zeuXMHEyZMwBtvvAFvb+8aj+E+WbuGzGUl7pe/Mvb/+91330Xr1q0RHh6O5ORk3LlzR99mzX2y2fyQrk6nw7PPPos+ffogPDy81m1++uknLFy4EAkJCfr7rl27ZjDRAPTr165dU9ymoKAAv/zyCxwcHMyZyn3X0LkEgAULFiA6OhqOjo746KOPMG3aNNy+fRszZ84EwLms6u9//zvmzp2LwsJCBAcHY9++fdBqtQC4X1bX0HkEuE9WZ2wun3vuOfTu3RsjR46s9XHcJ2tq6FwC3C+rMjaPEyZMQEBAAHx9fXHq1Ck8//zzyMrKwvbt2wFYd59sNsXS9OnTcebMGRw+fLjW9oKCAgwbNgydOnVCSkqKdQfXxJgyly+99JL+3926dUNhYSGWLl2qfwFoaZTmcuLEifj973+PnJwcvPbaa3j00Udx5MgR2Nvb34eRNm6mzCP3SUO1zeWuXbtw4MABfPXVV/dxZE2PKXPJ/fJXxv5/V30zHhERAR8fHzz88MP49ttv8eCDD1p1jM3iY7jExESkpaXh4MGDaNeuXY32W7duYfDgwXBxccGOHTtga2urb/P29kZubq7B9pXrlYdPjW3j6urarKp7wLS5rE1UVBT+97//obi4GADnsio3NzcEBQWhX79++Ne//oXz589jx44dALhfVmXKPNaG+2TNuTxw4AC+/fZbuLu7Q6PR6D9aHzt2LAYMGACA+2R1psxlbVrqfnmv/99VRUVFAQAuXrwIwLr7ZJMulkQEiYmJ2LFjBw4cOIAOHTrU2KagoACDBg2CVqvFrl27arxr79WrF06fPo3r16/r79u3bx9cXV3RqVMn/Tb79+83eNy+ffvQq1cvC2R1f5hjLmtz8uRJeHh46H/skHNp/DEion+h5H5pnnmsDffJmnM5b948nDp1CidPntTfAGDFihVYv349AO6Tlcwxl7VpaftlQ/5/V86lj48PACvvk/U6HbyRmTp1qri5ucmhQ4cMvlp4584dERHJz8+XqKgoiYiIkIsXLxpsU/0r2oMGDZKTJ0/K3r17pU2bNrV+9XDOnDly7tw5eeONN5rVVzhFzDOXu3btkrfffltOnz4tFy5ckL///e/i6OgoL7/8sj4O51Lk22+/lcWLF8vx48flu+++kyNHjsjw4cOlVatWkpubKyLcL0XMM4/cJ++611zWBkYuHdCS90kR88wl98t7z+PFixdlwYIFcvz4cbl8+bL85z//kY4dO0q/fv30fVhzn2zSxRKAWm/r168XEZGDBw8a3eby5cv6fq5cuSJDhgwRBwcHad26tcyePdvg6/CVfXXt2lW0Wq107NhRH6O5MMdc7tmzR7p27SrOzs7i5OQkXbp0kTVr1kh5eblBrJY+l99//70MGTJE2rZtK7a2ttKuXTuZMGGCnD9/3qCflr5fmmMeuU/eda+5NPaY6l9tb+n7pIh55pL75b3n8erVq9KvXz9p1aqV2NnZSWBgoMyZM8fgOksi1tsnVRWDJiIiIqJaNOlzloiIiIgsjcUSERERkQIWS0REREQKWCwRERERKWCxRERERKSAxRIRERGRAhZLRERERApYLBEREREpYLFERM3egAED8Oyzz97vYRBRE8UreBNRs5eXlwdbW1u4uLjc76EQURPEYomIiIhIAT+GI6ImYcCAAUhMTERiYiLc3NzQunVrvPTSS6h8v/f3v/8dQUFBsLe3h5eXFx555BGDx1b9GG7Tpk3o0aMHXFxc4O3tjQkTJuD69evWTomImgjN/R4AEVFd/eMf/8Cf/vQnfPHFFzh+/DgSEhLg7++Pbt26YebMmdi0aRN69+6NvLw8/Pe//zXaT2lpKRYuXIjg4GBcv34ds2bNQnx8PHbv3m3FbIioqeDHcETUJAwYMADXr1/H119/DZVKBQCYN28edu3ahb/85S+YMmUK/ve//9V6XtKAAQPQtWtXrFy5sta+jx8/jsjISNy6dQvOzs6WTIOImiB+DEdETcbvfvc7faEEAL169cKFCxfw8MMPIyAgAB07dsSkSZPw7rvv4s6dO0b7OXHiBIYPHw5/f3+4uLigf//+AICrV69aPAcianpYLBFRk+fs7Iwvv/wSW7ZsgY+PD15++WV06dIFN2/erLFtYWEhYmNj4erqinfffRfHjh3Djh07AAAlJSVWHjkRNQUsloioyTh69KjB+ueff46goCDY2NhAo9EgJiYGf/3rX3Hq1ClcuXIFBw4cqNHH+fPncePGDaSmpuKhhx5CSEgIT+4mIkU8wZuImoyrV69i1qxZePrpp/Hll1/i9ddfx7Jly5CWloZLly6hX79+8PDwwO7du6HT6RAcHFyjD39/f2i1Wrz++ut45plncObMGSxcuPA+ZENETQWPLBFRkzF58mT88ssv6NmzJ6ZPn46kpCQkJCTA3d0d27dvR3R0NEJDQ7FmzRps2bIFYWFhNfpo06YNNmzYgG3btqFTp05ITU3Fa6+9dh+yIaKmgt+GI6Im4V7faCMishQeWSIiIiJSwGKJiIiISAE/hiMiIiJSwCNLRERERApYLBEREREpYLFEREREpIDFEhEREZECFktEREREClgsERERESlgsURERESkgMUSERERkYL/D2o2GQ0yXb/QAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "prod_wells =['Prod {}'.format(i+1) for i in range(num_prod)]\n",
    "days = ['{}'.format((i+1)*100) for i in range(len(ACTIONS))]\n",
    "plt.imshow(Controls_opt[:, :num_prod].transpose())\n",
    "ax = plt.gca()\n",
    "ax.set_yticks(np.arange(num_prod))\n",
    "ax.set_xticks(np.arange(len(ACTIONS))+0.5)\n",
    "ax.set_yticklabels(labels=prod_wells)\n",
    "ax.set_xticklabels(labels=days, rotation=90, ha='left')\n",
    "# ax.grid(color='b', linestyle='-', linewidth=1)\n",
    "cbar = plt.colorbar(orientation=\"horizontal\")\n",
    "cbar.ax.set_xlabel('psia')\n",
    "# plt.colorbar()\n",
    "plt.clim(2200,2500)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "ae141f7f-d327-4d7b-badd-66502066d4e8",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAD1CAYAAACRBbVeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3wElEQVR4nO3deVxU9foH8M+ZYZVNQNmuLG79cIEkUcMVvaTmlmUZaqGV2aJX0V+5tGlWQtZVu9pP02talku3zDa1CDUvuC+YZZqWC6mAW4CgoMzz+0OYHGBYZpiZw/B5v17zOpzzPfM83+8cODxzzpkziogIiIiIiFRKY+sOEBEREVWFxQoRERGpGosVIiIiUjUWK0RERKRqLFaIiIhI1VisEBERkaqxWCEiIiJVc7B1B+qCTqfDuXPn4OHhAUVRbN0dIiIiqgERQX5+PoKCgqDRGD9+YhfFyrlz5xAcHGzrbhAREZEJMjMz0axZM6PtdlGseHh4AADSAbhbMM+dmG7B6ERERA1NEYD5+v/jxthFsVJ26scdQNXDNZezRaMTERE1RNVdwsELbImIiEjVWKwQERGRqrFYISIiIlVjsUJERESqxmKFiIiIVM3ixUpYWBgWLFhg6TRERERkp2pVrIwZMwZDhw6tVYK9e/di3LhxRtt//vlnDBs2DGFhYVAUhYUNERERGbD4kZWmTZuiUaNGRtsLCwvRokULJCcnIyAgwNLdISIionrGrGIlNjYWEydOxNSpU+Hj44OAgADMmjXLYJ3qTgN16tQJb731FuLj4+HsXLObrhUVFSEvL8/gQURERPbJ7CMrH3zwAdzc3LB7927MnTsXs2fPRkpKSl30zaikpCR4eXnpH/xeICIiIvtldrESGRmJmTNnonXr1khISEB0dDRSU1Prom9GzZgxA7m5ufpHZmamRfMRERGR7Zj93UCRkZEG84GBgcjJyTE3bJWcnZ1rfMqIiIiI6jezj6w4OjoazCuKAp1OZ25YIiIiIgC8KRwRERGpnNmngcxVXFyMI0eO6H8+e/YsMjIy4O7ujlatWtm4d0RERGRrNj+ycu7cOURFRSEqKgrnz5/H22+/jaioKIwdO9bWXSMiIiIVqNWRlZUrVxrMb9u2rcI6GzZsMJgvKiqCu7u70ZhhYWEQkdp0g4iIiBoQi50GKiwsRHp6OrKzs9GuXTtLpSEiIiI7Z7HTQEuXLkV8fDwSExMRExNjqTRERERk5xSxg3MweXl58PLywiEAHhbM0wIzLRidiIiooSkCkIzc3Fx4enoaXcvmF9gSERERVcXmH12uS2F9AE9Ljug7C8YmIiKiSvHIChEREakaixUiIiJSNRYrREREpGosVoiIiEjVWKwQERGRqrFYISIiIlVjsUJERESqxmKFiIiIVM3ixUpYWBgWLFhg6TRERERkp2pVrIwZMwZDhw6tVYK9e/di3LhxRtuXLVuGHj16wNvbG97e3oiLi8OePXtqlYOIiIjsl8WPrDRt2hSNGjUy2r5t2zaMGDECW7duxc6dOxEcHIy+ffvi7Nmzlu4aERER1QNmFSuxsbGYOHEipk6dCh8fHwQEBGDWrFkG61R3Gujjjz/Gs88+iw4dOiA8PBz//ve/odPpkJqaavQ5RUVFyMvLM3gQERGRfTL7yMoHH3wANzc37N69G3PnzsXs2bORkpJicrzCwkLcuHEDPj4+RtdJSkqCl5eX/hEcHGxyPiIiIlI3s4uVyMhIzJw5E61bt0ZCQgKio6OrPCpSnWnTpiEoKAhxcXFG15kxYwZyc3P1j8zMTJPzERERkbo5mBsgMjLSYD4wMBA5OTkmxUpOTsbatWuxbds2uLi4GF3P2dkZzs7OJuUgIiKi+sXsYsXR0dFgXlEU6HS6Wsd5++23kZycjO+//75CAUREREQNl9nFSl2YO3cu3njjDXz77beIjo62dXeIiIhIRWxerLz55pt45ZVXsHr1aoSFhSErKwsA4O7uDnd3dxv3joiIiGzN5rfbX7x4MYqLi/Hggw8iMDBQ/3j77bdt3TUiIiJSgVodWVm5cqXB/LZt2yqss2HDBoP5oqKiKo+QnDp1qjZdICIiogbGYqeBCgsLkZ6ejuzsbLRr185SaYiIiMjOWew00NKlSxEfH4/ExETExMRYKg0RERHZOUVExNadMFdeXh68vLyQ2wfwtOAlw8p3My0XnIiIqMEpApCM3NxceHp6Gl3L5hfYEhEREVXF5h9drktJWwDe15aIiMi+8MgKERERqRqLFSIiIlI1FitERESkaixWiIiISNVYrBAREZGqsVghIiIiVWOxQkRERKrGYoWIiIhUzeLFSlhYGBYsWGDpNERERGSnalWsjBkzBkOHDq1Vgr1792LcuHFG29evX4/o6Gg0btwYbm5u6NChA1atWlWrHERERGS/LH67/aZNm1bZ7uPjgxdffBHh4eFwcnLC119/jcceewx+fn7o16+fpbtHREREKmfWaaDY2FhMnDgRU6dOhY+PDwICAjBr1iyDdao7DRQbG4v7778fbdq0QcuWLTFp0iRERkYiLS3N6HOKioqQl5dn8CAiIiL7ZPY1Kx988AHc3Nywe/duzJ07F7Nnz0ZKSopJsUQEqampOHbsGHr27Gl0vaSkJHh5eekfwcHBpnafiIiIVM7sYiUyMhIzZ85E69atkZCQgOjoaKSmptYqRm5uLtzd3eHk5ISBAwdi4cKFuOeee4yuP2PGDOTm5uofmZmZ5g6DiIiIVMrsa1YiIyMN5gMDA5GTk1OrGB4eHsjIyMDVq1eRmpqKKVOmoEWLFoiNja10fWdnZzg7O5vaZSIiIqpHzC5WHB0dDeYVRYFOp6tVDI1Gg1atWgEAOnTogF9++QVJSUlGixUiIiJqOFR5UzidToeioiJbd4OIiIhUwOIfXa5OUlISoqOj0bJlSxQVFWHjxo1YtWoVFi9ebOuuERERkQrYvFgpKCjAs88+iz/++AOurq4IDw/HRx99hIcfftjWXSMiIiIVUERELJkgMDAQr732GsaOHWuxHHl5efDy8sJ0AJa87PZVzLRgdCIiooamCEAycnNz4enpaXQtix1ZKSwsRHp6OrKzs9GuXTtLpSEiIiI7Z7ELbJcuXYr4+HgkJiYiJibGUmmIiIjIzlnsyEpiYiISExMtFZ6IiIgaCFV+dJmIiIiojM0/DVSXUnK/gdbTzXIJlK2Wi01ERESV4pEVIiIiUjUWK0RERKRqLFaIiIhI1VisEBERkaqxWCEiIiJVY7FCREREqsZihYiIiFTN4sVKWFgYFixYYOk0REREZKdqVayMGTMGQ4cOrVWCvXv3Yty4cTVad+3atVAUpdY5iIiIyH5Z/A62TZs2rdF6p06dwnPPPYcePXpYuEdERERUn5h1Gig2NhYTJ07E1KlT4ePjg4CAAMyaNctgnZqcBiopKcGoUaPw6quvokWLFuZ0iYiIiOyM2desfPDBB3Bzc8Pu3bsxd+5czJ49GykpKbWKMXv2bPj5+eGJJ56o0fpFRUXIy8szeBAREZF9MrtYiYyMxMyZM9G6dWskJCQgOjoaqampNX5+Wloali9fjmXLltX4OUlJSfDy8tI/goODTek6ERER1QN1UqzcLjAwEDk5OTV6bn5+Ph599FEsW7YMTZo0qXHOGTNmIDc3V//IzMysVZ+JiIio/jD7AltHR0eDeUVRoNPpavTc3377DadOncLgwYP1y8qe6+DggGPHjqFly5YVnufs7AxnZ2czek1ERET1hcU/DVSV8PBwHD582GDZSy+9hPz8fLzzzjs8vUNERES2LVZcXFzQvn17g2WNGzcGgArLiYiIqGHi7faJiIhI1Wp1ZGXlypUG89u2bauwzoYNGwzmi4qK4O7ubnIOIiIiatgsdhqosLAQ6enpyM7ORrt27SyVhoiIiOycxU4DLV26FPHx8UhMTERMTIyl0hAREZGdU0REbN0Jc+Xl5cHLywsdc7+B1tPNYnn2KFstFpuIiKjhKQKQjNzcXHh6ehpdixfYEhERkarZ9KPLdaXs4FBJXqGFMxVZOD4REVFDcuv/anUneeyiWMnPzwcAZAQ/ZOOeEBERUW1dunQJXl5eRtvt4poVnU6Hc+fOwcPDA4qiWCRHXl4egoODkZmZWeV5NbXGZw515bCHMTCHeuIzh7py2MMYrJUjNzcXISEhuHLliv6msJWxiyMrGo0GzZo1s0ouT09Pi200a8RnDnXlsIcxMId64jOHunLYwxislUOjqfoSWl5gS0RERKrGYoWIiIhUjcVKDTk7O2PmzJlwdnaul/GZQ1057GEMzKGe+Myhrhz2MAa15bCLC2yJiIjIfvHIChEREakaixUiIiJSNRYrREREpGosVoiIiEjVWKwQERGRqrFYISIiIlWzi9vtq9nNmzfx888/IysrCwAQEBCAtm3bwtHR0cY9qx17GUdWVhZ2795tMI4uXbogICDAxj2rOWuMwR5eJyKyvaKiW9+qbO69Whp8sWKpnbJOp8Mrr7yCd999F7m5uQZtXl5emDBhAl599dVqvw+hJixZSNjLOAoKCvDUU09h7dq1UBQFPj4+AIDLly9DRDBixAi89957aNSoUYMeg7VepzK5ubkGr1VV37pqCmsV2ZYehzVy2MMY7GV7W3oclo6fkpKC+fPnY+fOncjLywNw6/uFYmJiMGXKFMTFxdU+qDRQV69elVGjRolWqxUHBwfx8/MTPz8/cXBwEK1WK4888ogUFBSYHP/555+Xpk2bypIlS+TkyZNSWFgohYWFcvLkSXnvvffEz89Ppk6datYYSkpK5MUXX5TGjRuLoigGj8aNG8tLL70kJSUlZuWwl3E88cQT0rp1a9m8ebPcvHlTv/zmzZvy7bffyh133CFjx45t0GOwVg4RkWXLlkmbNm1Eo9EYPNq0aSP//ve/zY5vje1hjXFYI4c9jMFetrelx2GN12nlypXi4OAg8fHxsmLFCtm4caNs3LhRVqxYISNGjBBHR0f58MMPax23wRYrlt4p+/v7y+bNm422b968Wfz8/EyOL2KdQsJextG4cWNJT0832p6WliaNGzc2Ob49jMFaOebOnSuNGjWS6dOny9atW+XIkSNy5MgR2bp1q8yYMUPc3NzkrbfeMiuHNbaHNcZh6Rz2MAYR+9nelh6HNV6n1q1by6JFi4y2v/vuu9KqVatax22wxYqld8qNGjWSH3/80Wj7oUOHxM3NzeT4ItYpJOxlHJ6enrJ3716j7Xv27BFPT0+T49vDGKyVIyQkRNatW2e0fe3atRIcHGxWDmtsD2uMw9I57GEMIvazvS09Dmu8Ts7OznL06FGj7UePHhUXF5dax22wnwbS6XRwcnIy2u7k5ASdTmdy/NjYWDz33HO4ePFihbaLFy9i2rRpiI2NNTk+AOTn5yMoKMhoe2BgIAoKCszKYS/jGDRoEMaNG4eDBw9WaDt48CCeeeYZDB482OT49jAGa+XIyclBRESE0faIiIhKf99qwxrbwxrjsHQOexgDYD/b29LjsMbr1K5dOyxfvtxo+/vvv4+2bdvWPrA5FVR9NnLkSImKipIDBw5UaDtw4IB07NhRRo0aZXL8M2fOSPv27cXBwUGioqKkf//+0r9/f4mKihIHBweJjIyUM2fOmDMEGTBggPTt21cuXLhQoe3ChQvSv39/GThwoFk57GUcly9flv79+4uiKOLj4yPh4eESHh4uPj4+otFo5N5775UrV66YHN8exmCtHD169JCEhAS5ceNGhbabN29KQkKC9OzZ06wc1tge1hiHpXPYwxhE7Gd7W3oc1nidtm7dKm5ubhIRESGTJ0+W5ORkSU5OlsmTJ0tkZKS4u7vLDz/8UOu4DfZbl69cuYKRI0fi22+/hbe3N/z8/ADcqp7//PNP9OvXD6tXr0bjxo1NzqHT6fDtt99i165dBlddx8TEoG/fvmZ/giYzMxMDBgzA0aNHERERAX9/fwBAdnY2Dh8+jLZt2+Lrr79GcHCwWXnsZRwA8Msvv1Q6jvDw8AY/Bmvl+PHHH9GvXz/cuHEDPXv2NHittm/fDicnJ3z33Xdo3769yTmssT2sMQ5L57CHMQD2s70tPQ5r7adOnTqFxYsXV7oPefrppxEWFlbrmA22WCljjR2/JVm6kLAWexiHPYzBWvLz8/HRRx9V+lqNHDkSnp6eZuewxvawxjgsncMexgDYz/a29Djq636qwRcrlrZnzx7s3LnT4Jeia9eu6NSpk417Vjv2MI7i4mJs2LCh0nHcd999VV7DpBbWGIM9vE5EZFvl7+USGBiINm3amHwvlwZdrFhyp5yTk4Nhw4YhPT0dISEhBofbzpw5g27duuGzzz7Tn34yhyULCXsZx4kTJ9CvXz+cO3cOXbp0MRjH7t270axZM2zatAmtWrVq0GOw1usEVLwhY2BgIDp37lynd8m1RpFtjXFYOoc9jAGwn+1t6XFYMr7FbiRq1pU09djx48elRYsW4uLiIr169ZLhw4fL8OHDpVevXuLi4iKtWrWS48ePmxx/2LBhEhMTU+lHuI4ePSpdu3aVBx980JwhSHZ2tnTv3l0URZHQ0FDp3LmzdO7cWUJDQ0VRFOnevbtkZ2eblcNexhEXFyf33Xef5ObmVmjLzc2V++67T/r27WtyfHsYg7VyWPqGjCLW2R7WGIelc9jDGETsZ3tbehzWeJ0sdS+XBlusWHqn7O7uXuknjcrs27dP3N3dTY4vYp1Cwl7G4erqKocPHzba/uOPP4qrq6vJ8e1hDNbKYY275Fpje9jDHYXtYQwi9rO9LT0Oa7xOlrqXS4MtViy9U/b19ZVt27YZbd+6dav4+vqaHF/EOoWEvYwjMDBQvvrqK6PtX375pQQGBpoc3x7GYK0c1rhLrjW2hz3cUdgexiBiP9vb0uOwxutkqRuJqvOyXyto3LgxTp06ZbT91KlTZn1s+eGHH8bo0aPx+eef67/ICQDy8vLw+eef47HHHsOIESNMjg/c+hbL22OXl5+fb/Y3XdrLOMaOHYuEhATMnz8fP/74I7Kzs5GdnY0ff/wR8+fPx5gxYzBu3DiT49vDGKyVw9I3ZASssz2sMQ5L57CHMQD2s70tPQ5rvE4Wu5GoORVUffbyyy+Lt7e3zJs3Tw4dOiRZWVmSlZUlhw4dknnz5omPj4/MnDnT5PjXr1+Xp59+WpycnESj0YiLi4u4uLiIRqMRJycneeaZZ+T69etmjeHZZ5+V0NBQWb9+vcHprNzcXFm/fr2EhYXJhAkTzMphbByKotSrcYiIJCcnS2BgoCiKov8CMkVRJDAwUN58802zYtvDGKyVw9I3ZBSxzvawxjgsncMexiBiP9vb0uOwxutkqRuJNthiRcQ6O/7c3FzZsmWLrF69WlavXi1btmyp9DoZU1ijICqTm5srqamp+nGkpqbWy3GIiPz++++yY8cO2bFjh/z+++91EtMexmCtHNa4S641toc93FHYHsYgYj/b29LjsNZ+qqSkRDZu3CivvPKKjBs3TsaNGyevvPKKbNq0yeRvdW7QH10uc/LkSYOPcDVv3tzGPaqdvLw87N+/32AMHTt2rJMbFBnj5OSEQ4cOoU2bNnUW0xbjqGt5eXnYt28fsrOzAdTPMZw/fx6LFy9GWloazp8/D41GgxYtWmDo0KEYM2YMtFptneSxxg0ZrfE7Vd/vKAwAR48erfBR1roegzVyWGN728M46uN+isWKEZmZmZg5cybef/99k2Ncu3YN+/fvh4+PT4Uvbrp+/To++eQTJCQkmNXPsp1Y2R/L0aNH8c4776CoqAiPPPII+vTpY1b8KVOmVLr8nXfewSOPPAJfX18AwLx588zKc7uCggJ88sknOHHiBIKCghAfH6/PY6oDBw7A29tbX4iuWrUKS5YswZkzZxAaGooJEyYgPj7e5Pj/+Mc/MHz4cPTo0cOsflZn0aJF2LNnDwYMGID4+HisWrUKSUlJ0Ol0eOCBBzB79mw4ODiYHH/fvn2Ii4tDq1at4Orqip07d2LkyJEoLi7Gt99+i7Zt22Lz5s3w8PCow1ERkb2p83u5mH28x05lZGSIRqMx+fnHjh3Tf3Zdo9FIz5495ezZs/r2rKwss+KLiGzatEmcnJzEx8dHXFxcZNOmTdK0aVOJi4uTPn36iFarldTUVLNyKIoiHTp0kNjYWIOHoijSqVMniY2Nld69e5uVo02bNnLp0iURuXW+MywsTLy8vKRTp07i4+Mjfn5+Zp+GiIyMlJSUFBERWbZsmbi6usrEiRNl8eLFkpiYKO7u7rJ8+XKT45dt59atW0tycrKcP3/erP5W5rXXXhMPDw8ZNmyYBAQESHJysvj6+srrr78uc+bMkaZNm8orr7xiVo5u3brJrFmz9POrVq2SLl26iMitw+AdOnSQiRMnmpVDRKSoqEjWrVsniYmJEh8fL/Hx8ZKYmCiffPKJFBUVmR2/OllZWfLqq6/WSazMzEzJz8+vsLy4uNikL2wr7+LFi7Jlyxb938iFCxckOTlZXn31VTly5IjZ8SvTvHlz+fXXXy0SW6fTyZYtW2Tp0qXy1VdfSXFxsdkxMzMzDb6cb/v27TJy5Ejp3r27jBo1Snbs2GF2jrfffltOnTpldpzqfPXVV/Lyyy9LWlqaiIikpqbKvffeK/369ZP33nvP7PiFhYWyfPlyeeyxx6R///4yYMAAmTBhgnz//fdmxxax3L1cGmyx8sUXX1T5mD9/vlnFxNChQ2XgwIFy4cIFOX78uAwcOFCaN28up0+fFpG6KVZiYmLkxRdfFBGRNWvWiLe3t7zwwgv69unTp8s999xjVo6kpCRp3rx5haLHwcFBfv75Z7Nil1EURf/LO2rUKOnatav8+eefIiKSn58vcXFxMmLECLNyuLq66nc0UVFRsnTpUoP2jz/+WNq2bWtyfEVR5Pvvv5dJkyZJkyZNxNHRUYYMGSJfffWVyedoy2vZsqV89tlnInKrmNZqtfLRRx/p29evXy+tWrUyK4erq6v89ttv+vmSkhJxdHSUrKwsERH57rvvJCgoyKwclr4hY02Y+2ZEROTcuXPSqVMn0Wg0otVq5dFHHzUoWurib3z37t3i5eUliqKIt7e37Nu3T5o3by6tW7eWli1biqurq+zfv9/k+O+8806lD61WKzNmzNDPm+Pee+/V/z1funRJunTpIoqiSNOmTUWj0Uh4eLjk5OSYlaNz5876j9xv2LBBNBqNDBkyRKZNmyb333+/ODo6VvmR/JpQFEW0Wq3ExcXJ2rVrLVJUL1myRBwcHKRjx47i6ekpq1atEg8PDxk7dqw89dRT4urqKgsWLDA5/vHjxyU0NFT8/PwkODhYFEWRgQMHSpcuXUSr1cpDDz1U6bdK14al7uXSYIuVsnfCiqIYfZizo/Hz8zP4rLlOp5Onn35aQkJC5LfffquTHZmnp6d+p15SUiIODg4GV6ofPnxY/P39zcohIrJnzx6544475H//93/174IsVay0aNFCvvvuO4P29PR0CQ4ONiuHr6+v7Nu3T0RubZuMjAyD9hMnTph1X53bx1BcXCzr1q2Tfv36iVarlaCgIHnhhRfM/gfs6uqqL3ZFRBwdHeWnn37Sz586dUoaNWpkVo7Q0FD9OzqRW/+QFUWRwsJCERE5efKkuLi4mJXDGnfJPXToUJWPdevWmf33l5CQIF26dJG9e/dKSkqKdOzYUaKjo+Xy5csicqtYURTFrBxxcXEyduxYycvLk7feekuaNWtmcOOxxx57TIYOHWpyfEVRpFmzZhIWFmbwUBRF/va3v0lYWJg0b97crDHc/rfxzDPPSNu2bfVHSjMzM6Vjx47y9NNPm5XDzc1NH7NLly6SnJxs0L5w4UKJiooyK4eiKLJixQq57777xNHRUXx9fWXSpElV3q+rttq2bat/I7VlyxZxcXGRd999V9++YsUKadOmjcnx7733XnnqqadEp9OJyK0Pmdx7770iIvLrr79KWFiYWZ+CFbHcvVwabLESFBQkGzZsMNp+8OBBs3ZmHh4elR6iHT9+vDRr1ky2b99eJ8XKiRMn9PPu7u4G74pPnTpl9j+WMvn5+ZKQkCCRkZFy+PBhcXR0rNNipeydVVBQUIU//roYxyOPPCJPPPGEiIg89NBD8tJLLxm0z5kzRyIiIkyOf/sO+XanT5+WmTNnSmhoqNnbu3nz5rJp0yYRubVj0Wg08sknn+jbv/nmGwkLCzMrx6RJk6R9+/ayadMm2bJli/Tu3VtiY2P17Zs3b5aWLVualcMad8mt6s1I2XJzt0dQUJDs3r1bP3/9+nUZPHiwdOjQQS5dulQnb0i8vb31+5Hi4mLRaDQGOffv3y9/+9vfTI7/1FNPSYcOHSrsqyz1ZuR//ud/5IsvvjBo//77780uiLy8vOTQoUMicuvNSNnPZU6cOGF2IX/7OLKzs+XNN9+U8PBw0Wg00qlTJ1m6dKnk5eWZlaOyNyS3/62cPHnSrHE0atTI4PReUVGRODo6ysWLF0Xk1lEpc/chlrqRaIMtVgYPHiwvv/yy0faMjAyz3hV16tRJPvzww0rbxo8fL40bNzZ7RxYZGan/5yVy60jK7Yfwtm/fbvZOoLw1a9aIv7+/aDSaOt2ZRURESFRUlLi7u8unn35q0P7DDz+YtUMWETl79qyEhYVJz549ZcqUKeLq6irdu3eXJ598Unr27ClOTk7yzTffmDWGqs7D6nS6CkeMauull16Spk2bytixY6V58+Yyffp0CQkJkcWLF8uSJUskODhYJk+ebFaO/Px8GT58uDg4OIiiKNK1a1eD64W+/fZbgwLJFNa4S66vr68sX75cTp06Venjm2++Mfvvz83NrcJ1HTdu3JChQ4dKZGSk/Pjjj3WS4+TJk/r58m9ITp8+bXYhv379egkODpaFCxfql9V1sVL2ZsTPz8/gaKDIrTcjzs7OZuUYMmSITJ8+XURE+vXrV+HU1bJly6R169Zm5TD2N759+3YZPXq0uLm5mXRn1tuVvZEVubXPUhTFYL+0bds2adasmcnxg4KCDE4bXrlyRRRF0RdZv//+u9nbwlL3cmmwxcr27dsN/tGXd/Xq1Sqrw+rMmTNHf3itMs8884zZh4gXL14sX3/9tdH2GTNm6I8m1KXMzEzZsGGDXL16tU7izZo1y+BR/nslnnvuOYmPjzc7z5UrV2TatGnStm1bcXFxEScnJwkNDZWRI0fK3r17zYodFhamf3diKSUlJfLGG2/IoEGDZM6cOaLT6WTNmjUSHBwsvr6+MmbMmDrbJteuXav0otG6YOkbMoqI9O3bV1577TWj7ea+GRERiYiIqFBYi/xVsISEhJhdrISHhxtcL/b111/rT8mJiOzatcusf15l/vjjD+nTp4/0799fzp8/X+fFyoABA+T+++8Xb2/vCoXqrl27zD5dfeTIEfH19ZWEhAR57bXXxN3dXR555BF54403JCEhQZydnWXFihVm5dBoNFW+IcnNza1wLVxtjR8/Xlq3bi2vv/66dO7cWUaPHi3h4eGyadMm2bx5s0RERMjjjz9ucvzRo0dLr1695JdffpHff/9dHn74YYPTY9u2bTP7lLul7uXSYIsVIrIdS9+Qcf369bJq1Sqj7ZcvX5aVK1ealWPq1KlGr625ceOGDBkyxOyCaNasWbJmzRqj7S+88II88MADZuUoo9PpZM6cORIQECBarbbOipUxY8YYPNatW2fQ/vzzz0u/fv3MznPixAmJj48XDw8P/Sk/R0dH6dq1q3z++edmx6/u6GlduHr1qjz55JPSvn17GTdunBQVFclbb70lTk5OoiiKxMbGmv2ty3fffbf+7y40NNTg+pL//Oc/8q9//asuhlLnN0TlfVaIyGbq8w0Zb968icLCQqM30rp58ybOnj2L0NBQi/WhsLAQWq3W7O9zud3+/fuRlpaGhIQEeHt711lcYwoKCqDVauHi4lIn8UQEOTk50Ol0aNKkCRwdHeskri1dv34dN27cqLP7Gx0/fhxFRUUIDw83675M1tRgv8iQiGyvefPmiImJQUxMjL5QyczMxOOPP27RvHWRw8HBoco7fp4/fx6vvvqqWTmqc+nSJTzzzDN1GrNjx46YNGkSvL29rbItLl++jGeffbbO4imKAn9/fwQGBuoLlfryO2WMi4sLPDw86ixH69at0b59+wqFSl3Fv3btGtLS0nDkyJEKbdevX8eHH35Y65g8skJEqnLo0CHcddddKCkpYQ4b57CHMTCHdeP/+uuv6Nu3L86cOQNFUdC9e3esWbMGQUFBAIDs7GwEBQXVOkf9OP5DRHbjyy+/rLL9999/Zw4r5bCHMTCHeuIDwLRp09C+fXvs27cPf/75JxITE9G9e3ds27YNISEhJsflkRUisiqNRgNFUVDVrkdRFLPe3TGHOuIzh7pyWGMM/v7++P777xEREQHg1jVEzz77LDZu3IitW7fCzc3NpCMrvGaFiKwqMDAQ69evh06nq/Rx4MAB5rBSDnsYA3OoJz5w63qV26+FURQFixcvxuDBg9GrVy/8+uuvJsVlsUJEVtWxY0fs37/faHt17/yYo+5y2MMYmEM98QEgPDwc+/btq7B80aJFuO+++zBkyBCT4vKaFSKyqueffx4FBQVG21u1aoWtW7cyhxVy2MMYmEM98QHg/vvvx5o1a/Doo49WaFu0aBF0Oh2WLFlS67i8ZoWIiIhUjaeBiIiISNVYrBAREZGqsVghIiIiVWOxQkRERKrGYoWIiIhUjcUKERERqRqLFSIiIlI1u7kp3PXr11FcXGzrbhAREVEtODk5wcXFpcp17KJYuX79OlxdmwK4auuuEBERUS0EBATg5MmTVRYsdlGs3DqichXAZABupUsdS6cO+GuYDpW0VTZ1rOV6VbXX9rnll5fNV0NbydPLQmgqabt9qjWy3Jz28m2ORtatrN+3L6+uvbI4te23sRymjN/oOqU3inYoMZgq2pu3ZkvntfqprnR603C5pnSKilOH0p81FdpKc0BX6XM1FdYrm9eVW698+1/Tsp+N98FwWj5GzftQcb3y65S1Gcv91/rVrVf563Z7H4z1s/rXrKZ9qPx106AE2tJvrdXeLM1RIqXzMJgqt1KhNFXFaUm5afnllc2XbyubGstVm9jVxavLWJU9r6p4dRmrln27Ufrc0s1dYf5mieFTbhgJZWx5Ve3l20yJYay9CMD8rCwUFxfbf7HyF2cAZYO9/T+PsQKgunlTn1dZsWJuDiOUKqZlP2tqOC3/z7v8tLp/1pUVSjUtBExtr6pYqem0pjnMyl11saI4ls47lJ/eateUzms0hv8wNbf9UzReINT0H6LxYqS69uqfa34Ow6lWv17FNk0NYyil62kMnvfXVCkXRyk3BRz0P5folxlOpTRG+eVl07K+isH6Fae60vX+6oO2pPTnm6X9NFasGPsHWd0/0KqeV9vn1GUfjK2rLfec8i922XJj+7/yccvvP2/fp5ata4xSbr78l9qIkeXG+ibAjdKCpnRz44ZSbr501eoKC1PmaxvToZp2bbn5muAFtkRERKRqLFaIiIhI1VisEBERkaqxWCEiIiJVY7FCREREqsZihYiIiFSNxQoRERGpGosVIiIiUjUWK0RERKRqLFaIiIhI1VisEBERkaqxWCEiIiJVY7FCREREqsZihYiIiFSNxQoRERGpGosVIiIiUjUWK0RERKRqLFaIiIhI1VisEBERkaqxWCEiIiJVY7FCREREquZg6w7UrSL8NaSS0qkDgJu3/YxK5m+Um3csN3UwMq1Je22fW3552bwRUsW07GddNVOldFpSbr5sCiPz5ftw+3z53Cg3XzbVlotRtrykXHv55eU3ofa29vLLjG0CbbmpsU1S0zgOVa1TOkCHEoOpaG8NRMrm9VNd6fRWu650uaIpnaLiVKd/ccpPb5auU/5FLM2ln94st1xn0K4rbf9r/q+prlysktKp1shUV25eo5+/We75uirXc0BJhXXK2ozl/mv96tYry2G8Dw5G+ll+XQejfajpVErjSGlcgbbk1s/am6VtpZu19FdKP1XKNvtNI9PKf12qni/fVjY1lqs2sauLV5exKnteVfFMjVXdtLp9tA64UboLKd3cFebLd+EGare8qvbybabEMNZehJqxi2JFRODu7o6rV+fbuiu2U/6Xn+qFshqt/E6AiKihCAgIgJOTU5Xr2EWxoigKrl69iszMTHh6etq6O2REXl4egoODuZ1UjNuofuB2qh+4nWrGyckJLi4uVa5jF8VKGU9PT/5C1APcTurHbVQ/cDvVD9xO5uMFtkRERKRqLFaIiIhI1eyiWHF2dsbMmTPh7Oxs665QFbid1I/bqH7gdqofuJ3qjiIi5T94SkRERKQadnFkhYiIiOwXixUiIiJSNRYrREREpGosVoiIiEjV6n2x8u677yIsLAwuLi7o0qUL9uzZY+suUTlJSUno1KkTPDw84Ofnh6FDh+LYsWO27hYZkZycDEVRkJiYaOuuUCXOnj2LRx55BL6+vnB1dUVERAT27dtn627RbUpKSvDyyy+jefPmcHV1RcuWLfHaa6+Bn2cxXb0uVtatW4cpU6Zg5syZOHDgAO68807069cPOTk5tu4a3eaHH37A+PHjsWvXLqSkpODGjRvo27cvCgoKbN01Kmfv3r147733EBkZaeuuUCWuXLmCbt26wdHREZs2bcKRI0fwz3/+E97e3rbuGt3mzTffxOLFi7Fo0SL88ssvePPNNzF37lwsXLjQ1l2rt+r1R5e7dOmCTp06YdGiRQAAnU6H4OBg/OMf/8D06dNt3Dsy5sKFC/Dz88MPP/yAnj172ro7VOrq1au466678H//9394/fXX0aFDByxYsMDW3aLbTJ8+Henp6fjvf/9r665QFQYNGgR/f38sX75cv2zYsGFwdXXFRx99ZMOe1V/19shKcXEx9u/fj7i4OP0yjUaDuLg47Ny504Y9o+rk5uYCAHx8fGzcE7rd+PHjMXDgQIO/KVKXL7/8EtHR0XjooYfg5+eHqKgoLFu2zNbdonK6du2K1NRU/PrrrwCAQ4cOIS0tDffee6+Ne1Z/1dsvMrx48SJKSkrg7+9vsNzf3x9Hjx61Ua+oOjqdDomJiejWrRvat29v6+5QqbVr1+LAgQPYu3evrbtCVfj999+xePFiTJkyBS+88AL27t2LiRMnwsnJCaNHj7Z196jU9OnTkZeXh/DwcGi1WpSUlOCNN97AqFGjbN21eqveFitUP40fPx4//fQT0tLSbN0VKpWZmYlJkyYhJSWl2q9pJ9vS6XSIjo7GnDlzAABRUVH46aefsGTJEhYrKvLJJ5/g448/xurVq9GuXTtkZGQgMTERQUFB3E4mqrfFSpMmTaDVapGdnW2wPDs7GwEBATbqFVVlwoQJ+Prrr7F9+3Y0a9bM1t2hUvv370dOTg7uuusu/bKSkhJs374dixYtQlFREbRarQ17SGUCAwPRtm1bg2Vt2rTBZ599ZqMeUWWef/55TJ8+HfHx8QCAiIgInD59GklJSSxWTFRvr1lxcnJCx44dkZqaql+m0+mQmpqKmJgYG/aMyhMRTJgwAZ9//jm2bNmC5s2b27pLdJu///3vOHz4MDIyMvSP6OhojBo1ChkZGSxUVKRbt24VPvb/66+/IjQ01EY9osoUFhZCozH896rVaqHT6WzUo/qv3h5ZAYApU6Zg9OjRiI6ORufOnbFgwQIUFBTgscces3XX6Dbjx4/H6tWr8cUXX8DDwwNZWVkAAC8vL7i6utq4d+Th4VHh+iE3Nzf4+vryuiKVmTx5Mrp27Yo5c+Zg+PDh2LNnD5YuXYqlS5faumt0m8GDB+ONN95ASEgI2rVrh4MHD2LevHl4/PHHbd21+kvquYULF0pISIg4OTlJ586dZdeuXbbuEpUDoNLHihUrbN01MqJXr14yadIkW3eDKvHVV19J+/btxdnZWcLDw2Xp0qW27hKVk5eXJ5MmTZKQkBBxcXGRFi1ayIsvvihFRUW27lq9Va/vs0JERET2r95es0JEREQNA4sVIiIiUjUWK0RERKRqLFaIiIhI1VisEBERkaqxWCEiIiJVY7FCREREqsZihYiIiFSNxQoR1ZiIYNy4cfDx8YGiKMjIyKhy/RMnTsDf3x+NGjVCenq6dTpJRMjMzERsbCzatm2LyMhI/Oc//7F1lxAWFobIyEh06NABvXv3rtVzWawQUY1t3rwZK1euxNdff43z588jKioKGzZsqHTdc+fO4Z577kH37t3xxBNPYNCgQTh8+HCF9WbNmoXw8HC4ubnB29sbcXFx2L17d4X1rl27Bjc3N5w4cQIrV66EoihQFAVarRbe3t7o0qULZs+ejdzc3LoeNlG94+DggAULFuDIkSP47rvvkJiYiIKCAlt3Czt27EBGRga2bt1aq+fV6y8yJCLr+u233xAYGIiuXbtWud6VK1fQr18/9OjRAytWrIBWq4W7uzv69euH9PR0g2/evuOOO7Bo0SK0aNEC165dw/z589G3b1+cOHECTZs21a+XkpKC0NBQtGrVCmlpafD09MSxY8cgIvjzzz+xY8cOJCUlYcWKFUhPT0dQUJDFXgcitQsMDERgYCAAICAgAE2aNMHly5fh5uZm456ZyLZfTURE9cXo0aMNvogyNDS0wryISEFBgcTExMhTTz0lOp3OIMbrr78urVq1kuzsbKN5cnNzBYB8//33Bssff/xxmTZtmoiIrFixQry8vCo8Nzs7W5o0aSKjRo3SL9u0aZN069ZNvLy8xMfHRwYOHCgnTpzQt/fu3VvGjx9vECcnJ0ccHR0r9IGoLvzwww8yaNAgCQwMFADy+eefV7reokWLJDQ0VJydnaVz586ye/duk/Lt27dP2rVrZ/P+hoWFyV133SXR0dHy0Ucf1aoPLFaIqEb+/PNPmT17tjRr1kzOnz8vOTk5+m/PLps3V1FRkbz11lvi5eUlFy5c0C8vKSkRPz8/2bFjh4gYL1ZERCZNmiQeHh5y8+ZNERH59NNP5bPPPpPjx4/LwYMHZfDgwRIRESElJSUiIvLxxx+Lt7e3XL9+XR9j3rx5EhYWVqHYIqpKWlqaFBcXV1j+888/S1ZWln5+48aN8uKLL8r69euN/vNfu3atODk5yfvvvy8///yzPPnkk9K4cWODQv/OO++Udu3aVXicPXtWv86lS5ekbdu2kp6eXmXfe/XqJStWrKi0ra76+8cff4iIyLlz56Rt27Zy6NChKvt0OxYrRFRj8+fP1x9BEZEq32XVxldffSVubm6iKIoEBQXJnj17DNrT09PFz89PX2BUVawsXrxYABg9enPhwgUBIIcPHxYRkWvXrom3t7esW7dOv05kZKTMmjXL7HFRw1FSUiJ33nmnPPjgg/pCWUTk6NGj4u/vL2+++WalzzP2N9S5c2eDI34lJSUSFBQkSUlJNe7T9evXpUePHvLhhx9Wu25VxYol+vvcc8/VKF8ZXmBLRDbXu3dvZGRkYMeOHejfvz+GDx+OnJwcffsXX3yBQYMGQaOpfpd1a38KKIoCADh+/DhGjBiBFi1awNPTE2FhYQCAM2fOAABcXFzw6KOP4v333wcAHDhwAD/99BPGjBlThyMke6fRaLBx40YcPHgQCQkJ0Ol0+O2339CnTx8MHToUU6dOrXGs4uJi7N+/H3FxcQbx4+LisHPnzhrFEBGMGTMGffr0waOPPlrr8dRGTfpbUFCA/Px8AMDVq1exZcsWtGvXrsY5WKwQkc25ubmhVatWuPvuu7F8+XI4ODhg+fLl+vYvv/wSQ4YMqVGsX375BZ6envD19QUADB48GJcvX8ayZcuwe/du/SeNiouL9c8ZO3YsUlJS8Mcff2DFihXo06cPQkND63CE1BAEBQVhy5YtSEtLw8iRI9GnTx/ExcVh8eLFtYpz8eJFlJSUwN/f32C5v78/srKyahQjPT0d69atw4YNG9ChQwd06NCh0k/j1YWa9Dc7Oxvdu3fHnXfeibvvvhsJCQno1KlTjXPw00BEZDJHR0eUlJTUeVydToeioiIAt46MnD59Gvfcc0+1z8vJycHq1asxdOhQaDQaXLp0CceOHcOyZcvQo0cPAEBaWlqF50VERCA6OhrLli3D6tWrsWjRorodEDUYISEhWLVqFXr16oUWLVpg+fLl+qN81tS9e3fodDqj7XPmzMGcOXP089euXcOuXbswYcIE/bIjR44gJCSkTvrTokULHDp0yOTns1ghIpOFhYUhNTUV3bp1g7OzM7y9vWv1/IKCArzxxhsYMmQIAgMDcfHiRbz77rs4e/YsHnroIQC3TgHFxcWhUaNGBs8VEWRlZek/urxz507MmTMHXl5eSE5OBgB4e3vD19cXS5cuRWBgIM6cOYPp06dX2pexY8diwoQJcHNzw/3332/Cq0F06wjCuHHjMHjwYOzduxeTJ0/GwoULaxWjSZMm0Gq1yM7OrhA7ICCgTvr59NNPY/jw4fr5UaNGYdiwYXjggQf0y2r68X9r9JengYjIZP/85z+RkpKC4OBgREVF1fr5Wq0WR48exbBhw3DHHXdg8ODBuHTpEv773//qz2d/8cUXlZ4CysvLQ2BgIP72t78hJiYG7733HkaPHo2DBw/q7y+h0Wiwdu1a7N+/H+3bt8fkyZPx1ltvVdqXESNGwMHBASNGjICLi0utx0J08eJF/P3vf0ebNm2wfv16pKamYt26dXjuuedqFcfJyQkdO3ZEamqqfplOp0NqaipiYmLqpK8+Pj5o1aqV/uHq6go/Pz+DZQ4ONTueYY3+8tNARKRaFy5cEAcHB4OPfVrKyZMnRaPRyP79+y2ei+xPSUmJREdHy4ABA6SoqEi/PCMjQ3x8fGTevHn6Zfn5+XLw4EE5ePCgAJB58+bJwYMH5fTp0/p11q5dK87OzrJy5Uo5cuSIjBs3Tho3bmyxv4WqPg2khv6yWCEi1Tp27Jj861//smiO4uJiOX/+vIwaNUq6du1q0Vxk37777ju5du1aheUHDhyQzMxM/fzWrVsNbqhY9hg9erTB8xYuXCghISHi5OQknTt3ll27dlms71UVK2roryJS+jk/IqIGaNu2bejduzfuuOMOfPrpp4iIiLB1l4ioHBYrREREpGq8wJaIiIhUjcUKERERqRqLFSIiIlI1FitERESkaixWiIiISNVYrBAREZGqsVghIiIiVWOxQkRERKrGYoWIiIhUjcUKERERqRqLFSIiIlK1/wdLrb8zK1JsqQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.ticker as ticker\n",
    "def fmt(x, pos):\n",
    "    a, b = '{:.2e}'.format(x).split('e')\n",
    "    b = int(b)\n",
    "    return r'${} \\times 10^{{{}}}$'.format(a, b)\n",
    "\n",
    "inj_wells =['Inj {}'.format(i+1) for i in range(num_inj)]\n",
    "days = ['{}'.format((i+1)*100) for i in range(len(ACTIONS))]\n",
    "plt.imshow(Controls_opt[:, num_prod:num_well].transpose())\n",
    "plt.yticks(np.arange(num_inj), labels=inj_wells)\n",
    "plt.xticks(np.arange(0.5,20.5,1), labels=days, rotation=90, ha='right')\n",
    "cbar = plt.colorbar(orientation=\"horizontal\")\n",
    "cbar.formatter.set_powerlimits((0, 0))\n",
    "cbar.formatter.set_useMathText(True)\n",
    "cbar.ax.set_xlabel('ft^3/Day')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "08045afe-3180-435a-86fa-4da3946395d4",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[100000.0859375, 100000.078125 , 100000.09375  , 100000.0390625],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ],\n",
       "       [100000.       , 100000.       , 100000.       , 100000.       ]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Controls_opt[:, num_prod:num_well]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "5c04017d-9f44-4c32-916a-aed024a5aa49",
   "metadata": {},
   "outputs": [],
   "source": [
    "hf_r = h5py.File(data_dir + yobs_data)\n",
    "Qrate_w = torch.tensor(np.array(hf_r.get('Qpro_w_norm_slt')).transpose((2,1,0)), dtype=torch.float32)\n",
    "Qrate_g = torch.tensor(np.array(hf_r.get('Qpro_g_norm_slt')).transpose((2,1,0)), dtype=torch.float32)\n",
    "BHP_inj = torch.tensor(np.array(hf_r.get('BHPinj_norm_slt')).transpose((2,1,0)), dtype=torch.float32)\n",
    "hf_r.close()\n",
    "yobs = torch.cat((Qrate_w,Qrate_g,BHP_inj),dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "6aa10bc4-f4c8-41f3-8ee3-9773dc2d46ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "hf_r = h5py.File(data_dir + ctrl_data, 'r')\n",
    "bhp0 = torch.tensor(np.array(hf_r.get('Pwf_norm_slt')).transpose((2,1,0)) , dtype=torch.float32)\n",
    "rate0 = torch.tensor(np.array(hf_r.get('Qinj_norm_slt')).transpose((2,1,0)) , dtype=torch.float32)\n",
    "hf_r.close()\n",
    "controls = torch.cat((bhp0,rate0),dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4af708c7-f820-4826-be0c-5e8faf66f56c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([600, 14, 20])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yobs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "f42fd2d8-f473-4d32-8a03-73739ee8478c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "npv_all =0\n",
    "for i_time in range(yobs.shape[2]):\n",
    "    alpha = np.power((1-0.05),i_time*100/365)\n",
    "    Qw = yobs[:,:num_prod, i_time]*Q_diff_w+Q_min\n",
    "    Qg = yobs[:,num_prod:num_prod*2, i_time]*Q_diff_g+Q_min\n",
    "    Qinj_g = controls[:, num_prod:-1, i_time]*Q_diff_g+Q_min\n",
    "    lf3toton =0.1167*4.536e-4 # convert lf^3 to ton \n",
    "    pv = np.sum((50-10)*lf3toton*Qinj_g.detach().cpu().numpy().squeeze())- np.sum(5*Qw.detach().cpu().numpy().squeeze()) - np.sum(50*lf3toton*Qg.detach().cpu().numpy().squeeze())\n",
    "    npv_all += pv*100*alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "4b0073c3-becb-4731-875c-d37d399033d0",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4329393110.151995"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(npv_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "0614e275-85f3-4b0b-9f35-ae3819eb0464",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "yobs_file_best = 'rate_norm_best2'\n",
    "yobs_data_best = yobs_file_best + '.mat'\n",
    "\n",
    "ctrl_file_best = 'controls_norm_best2'\n",
    "ctrl_data_best = ctrl_file_best + '.mat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "c635bddd-b7e3-47e9-9e32-25d401dd012c",
   "metadata": {},
   "outputs": [],
   "source": [
    "hf_r_best = h5py.File(data_dir + yobs_data_best, 'r')\n",
    "Qrate_w = torch.tensor(np.array(hf_r_best.get('Qpro_w_norm_best')).transpose((2,1,0)), dtype=torch.float32)\n",
    "Qrate_g = torch.tensor(np.array(hf_r_best.get('Qpro_g_norm_best')).transpose((2,1,0)), dtype=torch.float32)\n",
    "BHP_inj = torch.tensor(np.array(hf_r_best.get('BHPinj_norm_best')).transpose((2,1,0)), dtype=torch.float32)\n",
    "hf_r_best.close()\n",
    "yobs_best = torch.cat((Qrate_w,Qrate_g,BHP_inj),dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "e7d937fe-bfe7-4060-bd93-b1fca86464ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hf_r = h5py.File(data_dir + ctrl_data_best, 'r')\n",
    "bhp0 = torch.tensor(np.array(hf_r.get('Pwf_norm_best')).transpose((2,1,0)) , dtype=torch.float32)\n",
    "rate0 = torch.tensor(np.array(hf_r.get('Qinj_norm_best')).transpose((2,1,0)) , dtype=torch.float32)\n",
    "hf_r.close()\n",
    "controls_best = torch.cat((bhp0,rate0),dim=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "dc714704-dcf9-41f7-b574-6b5075697631",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "npv_all =0\n",
    "for i_time in range(yobs_best.shape[2]):\n",
    "    alpha = np.power((1-0.05),i_time*100/365)\n",
    "    Qw = yobs_best[:,:num_prod, i_time]*Q_diff_w+Q_min\n",
    "    Qg = yobs_best[:,num_prod:num_prod*2, i_time]*Q_diff_g+Q_min\n",
    "    Qinj_g = controls_best[:, num_prod:-1, i_time]*Q_diff_g+Q_min\n",
    "    lf3toton =0.1167*4.536e-4 # convert lf^3 to ton \n",
    "    pv = np.sum((50-10)*lf3toton*Qinj_g.detach().cpu().numpy().squeeze())- np.sum(5*Qw.detach().cpu().numpy().squeeze()) - np.sum(50*lf3toton*Qg.detach().cpu().numpy().squeeze())\n",
    "    npv_all += pv*100*alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "deef23cf-ecd7-4d4a-a468-25ae28de0ab6",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-10327732.66562316"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npv_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "4cc2b8cb-4fc1-4438-97b2-69fdb54ce1e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGsCAYAAAAPJKchAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABE5ElEQVR4nO3deXRU9f3/8dckJBMSkrBmITtLgiD7GhRwQRaphdoqUi2oiNVif1Jqa/H0W2u1jS3i0paKVgGtRRQVbHGhkVUksqOgkgVCEiALWzJJINvM/f0RMhhJQibbzSTPxzlzjnPzuTPv681kXtz7WSyGYRgCAAAwiYfZBQAAgPaNMAIAAExFGAEAAKYijAAAAFMRRgAAgKkIIwAAwFSEEQAAYCrCCAAAMBVhBAAAmIowAgAATOVWYWTbtm265ZZb1LNnT1ksFq1bt86l/X//+9/LYrFc9vDz82ueggEAwBW5VRgpLi7W4MGDtXTp0gbt/8gjjyg7O7vao3///rrtttuauFIAAFBfbhVGpk6dqqeeeko/+MEPavx5aWmpHnnkEYWFhcnPz0+jR4/Wli1bnD/v1KmTQkJCnI/c3Fx9/fXXmjt3bgsdAQAA+C63CiNX8tBDDykpKUmrV6/Wl19+qdtuu01TpkxRampqje1feeUVxcbGaty4cS1cKQAAqNJmwkhmZqZWrFihNWvWaNy4cerdu7ceeeQRXXvttVqxYsVl7UtKSvTvf/+bqyIAAJisg9kFNJWDBw/KbrcrNja22vbS0lJ169btsvZr165VYWGh5syZ01IlAgCAGrSZMFJUVCRPT0/t3btXnp6e1X7WqVOny9q/8sor+t73vqfg4OCWKhEAANSgzYSRoUOHym63Ky8v74p9QNLT07V582b95z//aaHqAABAbdwqjBQVFSktLc35PD09XQcOHFDXrl0VGxurO++8U7Nnz9aSJUs0dOhQnTp1Shs3btSgQYM0bdo0537Lly9XaGiopk6dasZhAACAb7EYhmGYXUR9bdmyRddff/1l2+fMmaOVK1eqvLxcTz31lF5//XWdOHFC3bt315gxY/TEE09o4MCBkiSHw6GoqCjNnj1bf/zjH1v6EAAAwHe4VRgBAABtT6OG9j799NOyWCxasGBBne3WrFmjfv36ycfHRwMHDtSHH37YmLcFAABtSIPDyO7du/XSSy9p0KBBdbbbsWOHZs2apblz52r//v2aMWOGZsyYoUOHDjX0rQEAQBvSoNs0RUVFGjZsmP7xj3/oqaee0pAhQ/T888/X2HbmzJkqLi7W+vXrndvGjBmjIUOGaNmyZfV6P4fDoZMnT8rf318Wi8XVcgEAgAkMw1BhYaF69uwpD4/ar380aDTN/PnzNW3aNE2cOFFPPfVUnW2TkpK0cOHCatsmT55c54q7paWlKi0tdT4/ceKE+vfv35BSAQCAybKyshQeHl7rz10OI6tXr9a+ffu0e/fuerXPycm5bGKx4OBg5eTk1LpPQkKCnnjiicu2Z2VlKSAgwLWCAQCAKWw2myIiIuTv719nO5fCSFZWlh5++GElJibKx8enUQXWZdGiRdWuplQdTEBAAGEEAAA3c6UuFi6Fkb179yovL0/Dhg1zbrPb7dq2bZv+/ve/q7S09LKp2ENCQpSbm1ttW25urkJCQmp9H6vVKqvV6kppAADATbk0mubGG2/UwYMHdeDAAedjxIgRuvPOO3XgwIHLgogkxcfHa+PGjdW2JSYmKj4+vnGVAwCANsGlKyP+/v66+uqrq23z8/NTt27dnNtnz56tsLAwJSQkSJIefvhhTZgwQUuWLNG0adO0evVq7dmzRy+//HITHQIAAHBnjZr0rCaZmZnKzs52Ph87dqxWrVqll19+WYMHD9Y777yjdevWXRZqAABA++QW08HbbDYFBgaqoKCADqwAALiJ+n5/N/mVEQAAAFcQRgAAgKkIIwAAwFSEEQAAYCrCCAAAMBVhBAAAmIowAgAATEUYAQCgnbI7DL1/4ITue22P7A7zph1zaTp4AADg/gzD0IavcvRsYopScoskSf/94qRmDA0zpR7CCAAA7YRhGNqcnKcl/0vRVydtkqQAnw66f3wvTewfbFpdhBEAANo4wzC0Pe20lvwvRQey8iVJft6emnttjOaO66XAjl6m1kcYAQCgDdt59IyWJKZoV/pZSZKPl4fmjI3WT8f3Vlc/b5Orq0QYAQCgDdqfeU7PJqbo09TTkiRvTw/9eHSkfnZ9bwX5+5hcXXWEEQAA2pBDJwr0XGKKNh7OkyR18LBo5sgIPXRDH4UGdjS5upoRRgAAaANScgv1XGKKPjqUI0nysEi3DgvXwzf2VURXX5OrqxthBAAAN5Z+uljPf5Ki/3xxUoYhWSzSLYN66uGJfdW7Ryezy6sXwggAAG4o6+x5/XVjqt7bf8I5YdnUq0O0YGKs4kL8Ta7ONYQRAADcSHbBBf19U5re3pOlcntlCLmhX5AW3hSrq8MCTa6uYQgjAAC4gVOFpfrHljT9e2emyiockqRr+3TXwkmxGhbZxeTqGocwAgBAK1ZwvlwvbTuiFZ8d04VyuyRpVHRXLZwUqzG9uplcXdMgjAAA0ApdKLNrxY50LdtyRLaSCknS4IjOemRSrK7t010Wi8XkCpsOYQQAgFakrMKht3Zn6q+b0nSqsFSSFBvcSY9MitNN/YPbVAipQhgBAKAVcDgM/eeLk3o2MUWZZ89LksK7dNTCm2I1fUiYPD3aXgipQhgBAMBEhmFo4zd5euZ/yTqcUyhJ6t7Jqv93Yx/dMTJS3h08TK6w+RFGAAAwyedHz2jxhmTtzTgnSfL36aAHJvTWPddEy9e7/XxFt58jBQCglTh0okCLNyRra8opSZUr6d49NkYPTOilzr6tYyXdlkQYAQCghRw9VaQliSn64MtsSZWL2N0xKkI/v6GvggNa10q6LYkwAgBAM8suuKAXPknVmr3HZXcYslik7w/uqYU3xSqqm5/Z5ZmOMAIAQDM5V1ymf2xJ02tJGc5ZU2/sF6RHJsfpqtAAk6trPQgjAAA0saLSCr36abr++elRFZVWTlg2Krqrfj0lTiOiu5pcXetDGAEAoImUVtj1788ztXRzms4Ul0mS+ocG6FdT4nRdbI82OWFZUyCMAADQSHaHoXX7T+jZxBSdyL8gSYrp7qeFN8Vq2sBQebThCcuaAmEEAIAGMgxDm5Pz9OePkpWcWzlhWXCAVQsmxupHw8Pl5dn2JyxrCi79X3rxxRc1aNAgBQQEKCAgQPHx8froo49qbb9y5UpZLJZqDx+f9jt0CQDQduzLPKeZL3+ue1fuUXJuoQJ8Oug3U/tp66+u16xRkQQRF7h0ZSQ8PFxPP/20+vbtK8Mw9Nprr2n69Onav3+/BgwYUOM+AQEBSk5Odj7nfhkAwJ2l5RVp8YbD2vBVriTJu4OH7rkmWj+b0EeBvl4mV+eeXAojt9xyS7Xnf/zjH/Xiiy/q888/rzWMWCwWhYSENLxCAABagZyCEj3/SYre3pMlhyF5WKQfDQ/Xgomx6tm5o9nlubUG9xmx2+1as2aNiouLFR8fX2u7oqIiRUVFyeFwaNiwYfrTn/5Ua3CpUlpaqtLSUudzm83W0DIBAGiUgvPl+sfWNK387JhKL84VclP/YP16cpz6BvubXF3b4HIYOXjwoOLj41VSUqJOnTpp7dq16t+/f41t4+LitHz5cg0aNEgFBQV65plnNHbsWH311VcKDw+v9T0SEhL0xBNPuFoaAABNpqTcrtd2HNPSzWmylVTOFTIyuot+M7WfhkcxV0hTshiGYbiyQ1lZmTIzM1VQUKB33nlHr7zyirZu3VprIPm28vJyXXXVVZo1a5aefPLJWtvVdGUkIiJCBQUFCghgxjoAQPOpsDv03r4Teu6TFGUXlEiSYoM76dEp/XRDvyD6PrrAZrMpMDDwit/fLl8Z8fb2Vp8+fSRJw4cP1+7du/XCCy/opZdeuuK+Xl5eGjp0qNLS0upsZ7VaZbVaXS0NAIAGMwxDiV/n6i8bkpWWVyRJ6hnoo4WT4vSDoWHyZK6QZtPoeUYcDke1qxh1sdvtOnjwoG6++ebGvi0AAE1mV/pZ/fnjw9qbcU6S1NnXSw9d30d3jYmSj5enydW1fS6FkUWLFmnq1KmKjIxUYWGhVq1apS1btmjDhg2SpNmzZyssLEwJCQmSpD/84Q8aM2aM+vTpo/z8fC1evFgZGRm67777mv5IAABwUXJOof7y8WFtPJwnSfLx8tDca2N0//jeCuzIMN2W4lIYycvL0+zZs5Wdna3AwEANGjRIGzZs0E033SRJyszMlIfHpUlezp07p3nz5iknJ0ddunTR8OHDtWPHjnr1LwEAoLmcyL+gZ/+Xovf2H5dhSJ4eFs0cGaGHb+yr4AAm52xpLndgNUN9O8AAAFCXc8Vl+seWNL2WlKGyi8N0p14dokcmx6l3j04mV9f2NFsHVgAA3M2FMrtW7EjXi1uOqPDiMN0xvbrq0Sn9NDSyi8nVgTACAGizKuwOvbvvuJ5LTFWOrXKYbr8Qf/1maj9NiO3BMN1WgjACAGhzahqmG9a5o345KVYzhoTJg2G6rQphBADQpuw5dlZPf3RYexim6zYIIwCANiE1t1B/2ZCsxK8rV9P18fLQvdfE6IHreivAh2G6rRlhBADg1nIKSvRcYorW7L20mm7lMN1YhQQyTNcdEEYAAG6p4EK5lm09ouXb052r6U4eEKxfTY5TnyBW03UnhBEAgFspKbfrX0kZ+vvmNBVcKJdUtZruVRoexTBdd0QYAQC4BbvD0Lr9J/RsYopO5F+QJPUNqlxN98arWE3XnRFGAACtmmEY2pJ8Sn/++LAO5xRKkkIDffSLm2L1w2HhrKbbBhBGAACt1oGsfD390Tf6/OhZSVKATwf97Po+untsNMN02xDCCACg1Tl6qkjP/C9ZHx7MkSR5d/DQPWOj9eB1vdXZ19vk6tDUCCMAgFYjz1aiFzamavXuLNkdhiwW6UfDwvWLm2LVs3NHs8tDMyGMAABMV1hSrn9uO6p/fpquC+V2SdKN/YL06yn9FBfCMN22jjACADBNaYVdq3Zm6m+b0nS2uEySNDSys34zpZ9G9+pmcnVoKYQRAECLczgM/eeLk1qSmKyss5XDdHv18NOvJ/fT5AHBDNNtZwgjAIAWYxiGPk09rac/Oqyvs22SpCB/q35xU6xuGx6uDp4eJlcIMxBGAAAt4svj+frzx4f1WdoZSZK/tYMeuK637r0mRh29GabbnhFGAADN6tjpYj3zv2St/zJbkuTt6aHZ8VGaf30fdfFjmC4IIwCAZnKqsFR/25SqVTszVXFxmO4PhoRp4aRYhXfxNbs8tCKEEQBAkyoqrbg4TPeozpdVDtO9Lq6Hfj25n/r3DDC5OrRGhBEAQJMoq3DozV2Z+uvGVJ25OEx3cHigHp3aT2N7dze5OrRmhBEAQKM4HIbWH8zWMxuSlXn2vCQppruffjU5TlOvDmGYLq6IMAIAaLDtqaf1548P6+CJAklS905WLZjYVzNHRsiLYbqoJ8IIAMBlB48X6M8fH9b2tNOSpE7WDvrp+F6699oY+Vn5aoFr+I0BANRb+sVhuh9cHKbr5WnRXWOi9ND1fdStk9Xk6uCuCCMAgCuqWk33rd1Z1Ybp/uKmWEV0ZZguGocwAgCola2kXC9vPapXt19aTff6uB769ZR+uiqUYbpoGoQRAMBlSsrteuPzDC3dnKZz58slsZoumg9hBADgZHcYWrv/hJ5LTNGJ/MrVdPsEddKvJsdpUn9W00XzIIwAAGQYhjZ+k6fFG5KVnFsoSQoJ8NHCm2J167AwVtNFsyKMAEA7t+fYWT390WHtyTgnSQrs6KWfXddbc8ZGy8eL1XTR/AgjANBOJecUavGGw/rkmzxJkrWDh+69NkYPjO+tQF8vk6tDe+LSdbcXX3xRgwYNUkBAgAICAhQfH6+PPvqozn3WrFmjfv36ycfHRwMHDtSHH37YqIIBAI1zIv+CHlnzhaa8sE2ffJMnTw+LZo2K0NZfXa9Hp/QjiKDFuXRlJDw8XE8//bT69u0rwzD02muvafr06dq/f78GDBhwWfsdO3Zo1qxZSkhI0Pe+9z2tWrVKM2bM0L59+3T11Vc32UEAAK7sXHGZlm5O0+ufZ6iswiFJmnp1iB6ZHKfePTqZXB3aM4thGEZjXqBr165avHix5s6de9nPZs6cqeLiYq1fv965bcyYMRoyZIiWLVtW7/ew2WwKDAxUQUGBAgIY1w4ArigurdDKHce0bMsRFZZWSJLie3XTo1P7aUhEZ3OLQ5tW3+/vBvcZsdvtWrNmjYqLixUfH19jm6SkJC1cuLDatsmTJ2vdunV1vnZpaalKS0udz202W0PLBIB2q7TCrlU7M7V0c5pOF5VJkvqHBujRqf00vm93humi1XA5jBw8eFDx8fEqKSlRp06dtHbtWvXv37/Gtjk5OQoODq62LTg4WDk5OXW+R0JCgp544glXSwMASKqwO/TevhN6YWOqc66QqG6+WnhTrG4Z1FMeHoQQtC4uh5G4uDgdOHBABQUFeueddzRnzhxt3bq11kDSEIsWLap2RcVmsykiIqLJXh8A2iKHw9BHh3K0JDFZR08VS6qcK+T/3dhXt40IlxdzhaCVcjmMeHt7q0+fPpKk4cOHa/fu3XrhhRf00ksvXdY2JCREubm51bbl5uYqJCSkzvewWq2yWln9EQDqwzAMbUk5pWc2JOurk5W3tbv4emn+9X1015go5gpBq9foeUYcDke1/h3fFh8fr40bN2rBggXObYmJibX2MQEAuGZX+lkt3nBYu49VTljWydpB88b10r3XRsvfhyG6cA8uhZFFixZp6tSpioyMVGFhoVatWqUtW7Zow4YNkqTZs2crLCxMCQkJkqSHH35YEyZM0JIlSzRt2jStXr1ae/bs0csvv9z0RwIA7cjB4wV65n/J2ppySlLlhGV3j43WAxN6q4uft8nVAa5xKYzk5eVp9uzZys7OVmBgoAYNGqQNGzbopptukiRlZmbKw+PSPcmxY8dq1apV+u1vf6vHHntMffv21bp165hjBAAaKC2vUM8mpujDg5UDATp4WHTHqAj9/Ia+Cg7wMbk6oGEaPc9IS2CeEQDtXdbZ83phY6re23dcDkOyWKQfDAnTgomxiuzma3Z5QI2afZ4RAEDzyyss0dJNaVq1K1Pl9sp/O07qH6xfTopTXIi/ydUBTYMwAgCtUP75Mr207ahWfJaukvLKqduv7dNdj0yOY9ZUtDmEEQBoRYpLK7Tis3S9tO2oCksqp24fGtlZv5oUp7F9uptcHdA8CCMA0AqUlNv1752Z+sfmNJ0prpy6vV+Ivx6ZFKcbrwpi6na0aYQRADBRaYVdq3dlaenmNOUVVs7ZFN3NV79g6na0I4QRADBBWYVDa/Zm6e+b0pRdUCJJCuvcUQ/d0Ec/Gs7U7WhfCCMA0IKqFrH766ZUHT9XuYhdSICP5t/QR7ePCJe1A1O3o/0hjABAC7A7DL1/4IT+ujFVx86clyR172TV/Ot7a9aoSNaPQbtGGAGAZuRwGFp/MFsvfJKiIxdX0u3q560HJ/TWXWOi1NGbEAIQRgCgGTgchjZ8laPnPklRSm6RJKmzr5fuH99Lc+Kj5Wflzy9QhU8DADQhwzD0yTd5ei4xRV9n2yRJ/j6VK+necw0r6QI1IYwAQBMwDENbUk7pucQUfXm8QJLUydpB914TrbnjeimwIyEEqA1hBAAawTAMfZZ2Rs8mJmtfZr4kqaOXp+6+Jlr3j+ulLn7e5hYIuAHCCAA00OdHz+jZxBTtSj8rSbJ28NDs+Cj9dEJvde9kNbk6wH0QRgDARXszzurZxBR9lnZGkuTdwUM/HhWpn13XW0EBPiZXB7gfwggA1NPOo2f0102pzhDi5WnRzJERmn99H4UGdjS5OsB9EUYAoA6GYSjpyBm9sDFVOy/ejungYdGPhofroRv6KLyLr8kVAu6PMAIANTAMQ9tST+uvG1O1N+OcJMnb00O3jQjXg9f1JoQATYgwAgDfYhiGNh3O0183puqLi0N0q/qE/HRCL27HAM2AMAIAqpwx9X9f5+pvm1L11cnKycp8vDx01+go3T++Fx1TgWZEGAHQrtkdhj46lK2/b0rT4ZxCSZKvt6d+Eh+leeN6MUQXaAGEEQDtkt1haP2XJ/W3TWlKy6tcO6aTtYPuHhute6+NUVcmKwNaDGEEQLtSYXdo3YGTWro5TemnK1fRDfDpoHuvjdE9Y2MU6Mu07UBLI4wAaBfKKhx6b99x/WPLEWWePS+pchXdeeN66SfxUQpgATvANIQRAG1aaYVdb+85rmVbjuhE/gVJUjc/b80b30t3jYlSJyt/BgGz8SkE0CaVlNu1elemlm09qhxbiSSph79VPx3fSz8eHSlfb/78Aa0Fn0YAbcq54jL96/MMvbbjmM4Ul0mSQgJ89MCEXrpjVKR8vDxNrhDAdxFGALQJWWfP65VPj+rtPcd1odwuSQrr3FEPXtdbt40Il7UDIQRorQgjANzaweMFemnbEX14MFsOo3LbgJ4Bun98L00bGKoOnh7mFgjgiggjANyOYRjamnJKL287qh1Hzji3j+vbXT8d31vX9Okmi8ViYoUAXEEYAeA2yu0O/feLk3p521HnbKmeHhbdMihU94/vrf49A0yuEEBDEEYAtHpFpRVavStTr25PV3ZB5cgYX29PzRoVqXuvjVFYZxavA9wZYQRAq5VnK9Hyz47p3zszVFhSIalyeO7dY6N11+goZksF2giXenYlJCRo5MiR8vf3V1BQkGbMmKHk5OQ691m5cqUsFku1h48Pq18CqF1aXqF+/c4XuvbPm7Vs6xEVllSoVw8//fmHA7X90es1//o+BBGgDXHpysjWrVs1f/58jRw5UhUVFXrsscc0adIkff311/Lz86t1v4CAgGqhhY5lAL7LMAztPnZOL287ok++yXNuHxndRfeP760b+wXJw4O/HUBb5FIY+fjjj6s9X7lypYKCgrR3716NHz++1v0sFotCQkIaViGANs3uMJT4dY6WbT2qA1n5kiSLRZrUP1j3j++t4VFdzC0QQLNrVJ+RgoICSVLXrl3rbFdUVKSoqCg5HA4NGzZMf/rTnzRgwIBa25eWlqq0tNT53GazNaZMAK1QcWmF3t13XMu3p+vYmcqF67w7eOhHw8N137Ux6tWjk8kVAmgpFsMwjIbs6HA49P3vf1/5+fnavn17re2SkpKUmpqqQYMGqaCgQM8884y2bdumr776SuHh4TXu8/vf/15PPPHEZdsLCgoUEMDQPcCdnci/oNd3HNObuzJlu9gpNbCjl2bHR2l2fLR6+FtNrhBAU7HZbAoMDLzi93eDw8iDDz6ojz76SNu3b681VNSkvLxcV111lWbNmqUnn3yyxjY1XRmJiIggjABubG/GOS3/LF0fH8qR/eJUqTHd/XTvNdG6dVi4/Fg9F2hz6htGGvTpf+ihh7R+/Xpt27bNpSAiSV5eXho6dKjS0tJqbWO1WmW18q8jwN2V2x366FCOlm9Pd/YHkaRr+nTT3GtjdF0snVIBuBhGDMPQz3/+c61du1ZbtmxRTEyMy29ot9t18OBB3XzzzS7vC8A95J8v05u7svR60jHnJGXeHTw0Y0hP3XttjPqFcIUTwCUuhZH58+dr1apVev/99+Xv76+cnBxJUmBgoDp2rJwBcfbs2QoLC1NCQoIk6Q9/+IPGjBmjPn36KD8/X4sXL1ZGRobuu+++Jj4UAGZLyyvSyh3penfvCefKud07WfWTMVG6c0ykunfiiieAy7kURl588UVJ0nXXXVdt+4oVK3T33XdLkjIzM+XhcWkutXPnzmnevHnKyclRly5dNHz4cO3YsUP9+/dvXOUAWgXDMLQ97bRe3Z6uLcmnnNuvCg3Q3GtjdMvgUFk7eJpYIYDWrsEdWFtSfTvAAGg5JeV2rdt/Qss/S1dKbpGkyvlBbuwXrLnXxmhMr65McAi0c83agRVA+5VnK9G/Ps/Qv3dm6mxxmSTJz9tTt42I0N1joxXdvfbZmAGgJoQRAPVy6ESBlm9P13+/PKlye+UF1bDOHXX32GjdPjJCgR1ZKwZAwxBGANTK4TC08XCe/vnpUe1KP+vcPiKqi+69NkaT+gerg6dL620CwGUIIwAuc6HMrnf2Zmn5Z8eUfrpYktTBw6Jpg0J17zUxGhzR2dwCAbQphBEATnm2Er2elKE3dmYo/3y5JCnAp4N+PDpKc8ZGKTSwo8kVAmiLCCMA9E22Ta9uT9d/DpxUmd0hSYrs6qt7r4nWbSMimKodQLPiLwzQThmGoa0pp/Tq9nR9mnrauX14VBfNGxejm/qHyJOp2gG0AMII0M6UlNv1/oETeuXTdKXmVc4P4mGRpl4dqrnjYjQssovJFQJobwgjQDtxpqhUb3yeqX99fkyniy7NDzJzZKTuuSZaEV19Ta4QQHtFGAHauLS8Ir26PV3v7Tuu0orK/iChgT6655po3TEqUgE+zA8CwFyEEaANMgxDSUfP6JVP07XpcJ5z+8CwQN03LkY3DwyVF/ODAGglCCNAG1JW4dAHB0/qlU/T9dVJm6RL68XMGxejUTGsFwOg9SGMAG1AwflyrdqVqdd2HFOOrUSS5OPloduGR+iea6LVq0cnkysEgNoRRgA3lnX2vF7dnq6392TpfJldktTD36q7x0brx6Mi1cXP2+QKAeDKCCOAG9qXeU6vfHpUHx/KkaNyzTrFBfvrvnEx+v6QnrJ28DS3QABwAWEEcBN2h6HEr3P0z0/TtTfjnHP7uL7dNW9cL43r253+IADcEmEEaOXOl1VozZ7jWv5ZujLOnJckeXlaNH1ImO4bF6N+IQEmVwgAjUMYAVqpXFuJXttxTP/emamCC5WL1gV29NJdYyI1Jz5aQQE+JlcIAE2DMAK0Mt9k2/TKp+n6zxcnVG6v7BAS1c1Xc6+N0Y+Gh8vXm48tgLaFv2pAK2AYhralntYrnx6ttmjdyOguum9cL028KphF6wC0WYQRwESlFXa9v/+kXtl+VCm531q0bmCo5o3rpSERnc0tEABaAGEEMMG54jK98XmGXkvK0OmiUkksWgeg/SKMAC0oJbdQr+04pnf3HVdJeeWidSEBlxatC+zIonUA2h/CCNDMKuwOffJNrl7bkaGko2ec2wf0DNC8cb00bRCL1gFo3wgjQDM5U1Sq1buz9O/PM3SyoHK9GA+LNKl/iGaPjVJ8r25MUgYAIowATe7L4/laueOY1n+RrTJ75a2Yrn7emjUqQj8eHaWwzh1NrhAAWhfCCNAESivs+vBgtl7bkaEDWfnO7YPCAzUnPlrTBoXKx4v1YgCgJoQRoBGyCy5o1c5MvbkrU6eLyiRJ3p4emjYoVHPGRjM0FwDqgTACuMgwDO1KP6vXko5pw1e5sl9cNjckwEd3jYnUzJGR6uFvNblKAHAfhBGgns6XVWjd/pN6PemYDucUOrePjumqOWOjdVP/YEbFAEADEEaAK8g4U6x/JWXo7T1ZspVUSJI6ennqB8PCNDs+ilVzAaCRCCNADRwOQ9tST+n1pAxtTs6TUXknRlHdfPWTMVG6bXiEAn2ZoAwAmgJhBPiWXFuJ1uzJ0urdWTp+7oJz+3VxPTRnbLQm9O0hDxasA4Am5dIN7oSEBI0cOVL+/v4KCgrSjBkzlJycfMX91qxZo379+snHx0cDBw7Uhx9+2OCCgaZmdxjadDhX817fo7FPb9Iz/0vR8XMXFODTQXOvjdGWR67TyntG6fq4IIIIADQDl66MbN26VfPnz9fIkSNVUVGhxx57TJMmTdLXX38tPz+/GvfZsWOHZs2apYSEBH3ve9/TqlWrNGPGDO3bt09XX311kxwE0BAn8y/ord1ZWrMnyzlDqiSNjO6iWaMidfNA5gYBgJZgMYyqu+GuO3XqlIKCgrR161aNHz++xjYzZ85UcXGx1q9f79w2ZswYDRkyRMuWLavX+9hsNgUGBqqgoEABAXQWRMNV2B3adDhPb+7K1NaUU7o4KlddfL1067BwzRoVoT5B/uYWCQBtRH2/vxvVZ6SgoECS1LVr11rbJCUlaeHChdW2TZ48WevWrat1n9LSUpWWljqf22y2xpQJKOvseb21O0tv78lSXuGl3634Xt10x6gITR4QwlUQADBJg8OIw+HQggULdM0119R5uyUnJ0fBwcHVtgUHBysnJ6fWfRISEvTEE080tDRAklRWUbla7pu7MrU97bRzREw3P2/9aES47hgZqZjuNd9eBAC0nAaHkfnz5+vQoUPavn17U9YjSVq0aFG1qyk2m00RERFN/j5om9JPF2v17ky9u/e4c4p2SRrXt7tmjYrUxKuC5d2ByckAoLVoUBh56KGHtH79em3btk3h4eF1tg0JCVFubm61bbm5uQoJCal1H6vVKquV6bRRf6UVdm34Kldv7sxU0tEzzu09/K26fUS4Zo6IVGQ3XxMrBADUxqUwYhiGfv7zn2vt2rXasmWLYmJirrhPfHy8Nm7cqAULFji3JSYmKj4+3uVige9KyyvU6l1ZenffcZ07Xy5Jslik62J76I5RkbqhXxBTtANAK+dSGJk/f75WrVql999/X/7+/s5+H4GBgerYsaMkafbs2QoLC1NCQoIk6eGHH9aECRO0ZMkSTZs2TatXr9aePXv08ssvN/GhoL0oKq3QB1+e1Fu7s7QvM9+5PTTQR7ePiNDtIyMU1rmjeQUCAFziUhh58cUXJUnXXXddte0rVqzQ3XffLUnKzMyUh8elf4mOHTtWq1at0m9/+1s99thj6tu3r9atW8ccI3CJYRjak3FOb+/O0gcHs3W+zC5J8vSw6Pq4IP14dIQmxAbJk0nJAMDtNGqekZbCPCPtV15hid7bd0Jv78nS0VPFzu29uvvp9pERunVomIICfEysEABQmxaZZwRoDuV2hzYfztPbe45rc3Ke7BdnJuvo5anvDQrV7SMjNCKqiywWroIAQFtAGEGrkZZXpDV7s/Tu3hM6XXRpYrJhkZ01c2SEpg3qqU5WfmUBoK3hLztMVVxaoQ++zNbbe7K0J+Occ3v3Tt66dVi4bhserr7BTM8OAG0ZYQQtzjAM7cs8p7d2Z2n9l5c6o3pYpOvjgnT7yAiG5AJAO0IYQYs5VViq9/Yd19t7snTkW51RY7r76bYR4frhsHAF0xkVANodwgiald1haFvKKa3alanNh/NU8a3OqDcPDNXMkREaGU1nVABozwgjaBani0r19p4srdqZqePnLji3D4mo7Iz6vUGh8vfxMrFCAEBrQRhBkzEMQ7vSz+rfOzP10aFsldsrr4IE+HTQj4ZH6I5REYqlMyoA4DsII2g0W0m51u47oX/vzFBKbpFz++CIzrpzdKRuGdRTHb09TawQANCaEUbQYIdOFOiNzzP0/oGTulBeOSKmo5enpg/pqTtHR2lgeKDJFQIA3AFhBC4pKbfrv1+c1Bs7M/VFVr5ze9+gTrprTJR+MCxMAfQFAQC4gDCCejlyqkirdmbqnb3HVXChXJLk5WnRlKtDddfoSI2K6cqIGABAgxBGUKtyu0OJX+fqjc8ztOPIGef28C4d9ePRkbpteIR6+FtNrBAA0BYQRnCZk/kXtHpXplbvzlJeYeUaMRaLdENckO4aE6XxsT3k6cFVEABA0yCMQFLlsNxPU0/r9aQMbTqcq4tzk6l7J6vuGFk5LDe8i6+5RQIA2iTCSDtXVuHQf744qX9uO6rk3ELn9jG9uuquMVGa1D9E3h1YIwYA0HwII+1UwYVyvbkrUys+S1eurfJWjJ+3p24bEaG7xkSqTxCTkwEAWgZhpJ05kX9By7ena/WuTBVfXC03yN+qe66J0Y9HRyqwI8NyAQAtizDSThw6UaB/fnpU67/Mlv1ih5DY4E6aN66Xpg8J41YMAMA0hJE2zDAMbU05pX9+elSfpV0amju2dzfdP76XJsT2YG4QAIDpCCNtUE2dUj09LJo2MFT3j++lq8OYph0A0HoQRtqQ2jql3jEqUvdcE83QXABAq0QYaQNO5F/Qiu3pWr07S0WlFZK+1Sl1VKQCfemUCgBovQgjbqyuTqnfH9JT1g6eJlcIAMCVEUbcjGEY2pZ6Wi9vO3JZp9R543vpOjqlAgDcDGHETZwuKtW7e4/rrd1ZOnq6WNKlTqnzxvXSwHA6pQIA3BNhpBVzOAxtTzut1bszlfh1rsrtlbdi/Lw9NXNkZafUiK50SgUAuDfCSCuUayvRmj1ZWr07S8fPXXBuHxweqDtGReqWwT3VycqpAwC0DXyjtRIVdoe2ppzSm7uytDk5z9kh1d+ng34wNEx3jIxU/54BJlcJAEDTI4yY7Pi583p7d5be3nNcObYS5/aR0V10x8hI3TwwVB29GRUDAGi7CCMmKLc7tPGbXL25K0vbUk/JqLwIoi6+Xrp1WLjuGBmhvsGsmgsAaB8IIy3o2Olird6dpXf2HtfpolLn9rG9u+mOUZGaPCCYuUEAAO0OYaSZlVbY9fGhHK3elaWko5fmBeneyarbRoRr5ogIRXf3M7FCAADM5XIY2bZtmxYvXqy9e/cqOztba9eu1YwZM2ptv2XLFl1//fWXbc/OzlZISIirb+82jpwq0qqdmXpv33GdO18uSbJYpPF9e2jWqAjdeFWwvDw9TK4SAADzuRxGiouLNXjwYN1777269dZb671fcnKyAgIujQYJCgpy9a3dRuaZ85ry/DbnvCChgT66bUSEbh8RzmJ1AAB8h8thZOrUqZo6darLbxQUFKTOnTu7vJ872pt5VuV2Q+FdOuoP0wdoQmyQPD2Yoh0AgJq02H2CIUOGKDQ0VDfddJM+++yzOtuWlpbKZrNVe7iTlNwiSdKE2B66oV8wQQQAgDo0exgJDQ3VsmXL9O677+rdd99VRESErrvuOu3bt6/WfRISEhQYGOh8RERENHeZTSo1t1CSFMvwXAAArqjZR9PExcUpLi7O+Xzs2LE6cuSInnvuOf3rX/+qcZ9FixZp4cKFzuc2m82tAknVlRHCCAAAV2bK0N5Ro0Zp+/bttf7carXKarW2YEVN53xZhTLPnpckxQZ3MrkaAABaP1PGlh44cEChoaFmvHWzS8urvCrSzc9b3Tq5Z6ACAKAluXxlpKioSGlpac7n6enpOnDggLp27arIyEgtWrRIJ06c0Ouvvy5Jev755xUTE6MBAwaopKREr7zyijZt2qT//e9/TXcUrUjVLZq+XBUBAKBeXA4je/bsqTaJWVXfjjlz5mjlypXKzs5WZmam8+dlZWX65S9/qRMnTsjX11eDBg3SJ598UuNEaG0BnVcBAHCNxTCqlmlrvWw2mwIDA1VQUFBt4rTW6J4Vu7Q5+ZSemnG17hoTZXY5AACYpr7f38xH3sQYSQMAgGsII02oqLRCJ/IvSGIkDQAA9UUYaUJV/UV6+FvV2dfb5GoAAHAPhJEmlOq8RcNVEQAA6osw0oRSLl4Z6RtEfxEAAOqLMNKEUi5OeBYXQhgBAKC+CCNNKCWnao4RbtMAAFBfhJEmUnChXDm2EklSH27TAABQb4SRJpKWV3lVJCTAR4EdvUyuBgAA90EYaSKsSQMAQMMQRppI1UiaOGZeBQDAJYSRJpLCAnkAADQIYaSJcJsGAICGIYw0gfzzZTpVWCpJ6suVEQAAXEIYaQJVV0XCOndUJ2sHk6sBAMC9EEaawKX+ItyiAQDAVYSRJkDnVQAAGo4w0gScC+QRRgAAcBlhpAmkXuwzwm0aAABcRxhppDNFpTpTXCZJ6hNEGAEAwFWEkUaqGkkT2dVXvt6MpAEAwFWEkUZiJA0AAI1DGGkkOq8CANA4hJFGovMqAACNQxhpBMMwlJJ38cpIEFdGAABoCMJII5wqKlX++XJ5WBhJAwBAQxFGGqHqFk1UNz/5eHmaXA0AAO6JMNIIyTlVt2i4KgIAQEMRRhohNY81aQAAaCzCSCNUTXjWl5E0AAA0GGGkgQzDYLVeAACaAGGkgXJtpSosqZCnh0W9eviZXQ4AAG6LMNJAyRevikR385W1AyNpAABoKJfDyLZt23TLLbeoZ8+eslgsWrdu3RX32bJli4YNGyar1ao+ffpo5cqVDSi1dUnlFg0AAE3C5TBSXFyswYMHa+nSpfVqn56ermnTpun666/XgQMHtGDBAt13333asGGDy8W2JqxJAwBA03B5zfupU6dq6tSp9W6/bNkyxcTEaMmSJZKkq666Stu3b9dzzz2nyZMnu/r2rUYKa9IAANAkmr3PSFJSkiZOnFht2+TJk5WUlFTrPqWlpbLZbNUerYlhGM7bNHFcGQEAoFGaPYzk5OQoODi42rbg4GDZbDZduHChxn0SEhIUGBjofERERDR3mS45kX9BxWV2eXlaFN2dkTQAADRGqxxNs2jRIhUUFDgfWVlZZpdUTdWaNDHd/eTl2Sr/FwIA4DZc7jPiqpCQEOXm5lbblpubq4CAAHXs2LHGfaxWq6xWa3OX1mB0XgUAoOk0+z/r4+PjtXHjxmrbEhMTFR8f39xv3WycnVeDCCMAADSWy2GkqKhIBw4c0IEDByRVDt09cOCAMjMzJVXeYpk9e7az/QMPPKCjR4/q17/+tQ4fPqx//OMfevvtt/WLX/yiaY7ABJcWyGMkDQAAjeVyGNmzZ4+GDh2qoUOHSpIWLlyooUOH6ne/+50kKTs72xlMJCkmJkYffPCBEhMTNXjwYC1ZskSvvPKK2w7rdTgMZ5+R2BCujAAA0FgWwzAMs4u4EpvNpsDAQBUUFCggIMDUWjLPnNf4xZvl7emhr/8wWR3owAoAQI3q+/3NN6mLqjqv9urhRxABAKAJ8G3qopQ81qQBAKApEUZclMo08AAANCnCiItSWK0XAIAmRRhxgd1hKC2v6soIYQQAgKZAGHFB5tnzKq1wyNrBQxFdfc0uBwCANoEw4oKqWzR9gjrJ08NicjUAALQNhBEXpNJfBACAJkcYcYFzTRrCCAAATYYw4oJLI2kY1gsAQFMhjNRThd2ho6eKJXFlBACApkQYqadjZ86rzO5QRy9PhXXuaHY5AAC0GYSReqrqvNo3uJM8GEkDAECTIYzUE51XAQBoHoSRerq0QB6dVwEAaEqEkXpKyam6TcOVEQAAmhJhpB7KKhxKP81IGgAAmgNhpB6OnSlWhcNQJ2sH9Qz0MbscAADaFMJIPXx7TRqLhZE0AAA0JcJIPVSNpInjFg0AAE2OMFIPlzqvMpIGAICmRhiph0vDerkyAgBAUyOMXEFphV0ZZ85LIowAANAcCCNXcPRUsewOQ/4+HRQcYDW7HAAA2hzCyBVUjaSJC/ZnJA0AAM2AMHIFKbnMvAoAQHMijFzBpQXyGEkDAEBzIIxcQWouI2kAAGhOhJE6lJTblXG2ciQNc4wAANA8CCN1SMsrkmFIXXy91KMTI2kAAGgOhJE6fLvzKiNpAABoHoSROtB5FQCA5kcYqQOdVwEAaH4NCiNLly5VdHS0fHx8NHr0aO3atavWtitXrpTFYqn28PHxaXDBLalqTZq+QYQRAACai8th5K233tLChQv1+OOPa9++fRo8eLAmT56svLy8WvcJCAhQdna285GRkdGoolvC+bIKZZ29IInbNAAANCeXw8izzz6refPm6Z577lH//v21bNky+fr6avny5bXuY7FYFBIS4nwEBwc3quiWkJZX2V+keydvdWMkDQAAzcalMFJWVqa9e/dq4sSJl17Aw0MTJ05UUlJSrfsVFRUpKipKERERmj59ur766qs636e0tFQ2m63ao6Ul53CLBgCAluBSGDl9+rTsdvtlVzaCg4OVk5NT4z5xcXFavny53n//fb3xxhtyOBwaO3asjh8/Xuv7JCQkKDAw0PmIiIhwpcwmkZrHSBoAAFpCs4+miY+P1+zZszVkyBBNmDBB7733nnr06KGXXnqp1n0WLVqkgoIC5yMrK6u5y7wMC+QBANAyOrjSuHv37vL09FRubm617bm5uQoJCanXa3h5eWno0KFKS0urtY3VapXVam4/jVTnHCOEEQAAmpNLV0a8vb01fPhwbdy40bnN4XBo48aNio+Pr9dr2O12HTx4UKGhoa5V2oIKS8p1Ip+RNAAAtASXroxI0sKFCzVnzhyNGDFCo0aN0vPPP6/i4mLdc889kqTZs2crLCxMCQkJkqQ//OEPGjNmjPr06aP8/HwtXrxYGRkZuu+++5r2SJpQVX+RIH+rOvt6m1wNAABtm8thZObMmTp16pR+97vfKScnR0OGDNHHH3/s7NSamZkpD49LF1zOnTunefPmKScnR126dNHw4cO1Y8cO9e/fv+mOookx8yoAAC3HYhiGYXYRV2Kz2RQYGKiCggIFBAQ0+/s9uf5rvbo9XfdcE63HbxnQ7O8HAEBbVN/vb9amqUEKV0YAAGgxhJEaMJIGAICWQxj5joIL5cqxlUiS+jKSBgCAZkcY+Y6qzquhgT4K8PEyuRoAANo+wsh3pFy8RcPMqwAAtAzCyHc4O68GcYsGAICWQBj5jtQ8RtIAANCSCCPfUXWbJjaEMAIAQEsgjHzLueIynSoslST15TYNAAAtgjDyLVX9RcI6d5Sf1eWZ8gEAQAMQRr4lJa9qsjOuigAA0FIII9/CAnkAALQ8wsi3sCYNAAAtjzDyLSmsSQMAQIsjjFx0uqhUZ4vLZLFIfRhJAwBAiyGMXFR1iyaii686enuaXA0AAO0HYeSi1FxG0gAAYAbCyEV0XgUAwByEkYsIIwAAmIMwIskwDOdImr7cpgEAoEURRiSdKixVwYVyeVik3j0IIwAAtCTCiC7NLxLVzU8+XoykAQCgJRFG9O3+IlwVAQCgpRFGROdVAADMRBjRpTDSlzACAECLa/dhxDAMJjwDAMBE7T6M5NhKVFhaIU8Pi2K6+5ldDgAA7U67DyNVI2miu/nK2oGRNAAAtLR2H0ZSL/YXiQuhvwgAAGZo92EkOedi59UgwggAAGZo92EkJa+q8yphBAAAM7TrMGIYhtKY8AwAAFM1KIwsXbpU0dHR8vHx0ejRo7Vr1646269Zs0b9+vWTj4+PBg4cqA8//LBBxTa1E/kXVFxml5enRdGMpAEAwBQuh5G33npLCxcu1OOPP659+/Zp8ODBmjx5svLy8mpsv2PHDs2aNUtz587V/v37NWPGDM2YMUOHDh1qdPGNVTW/SK/uneTl2a4vEgEAYBqXv4GfffZZzZs3T/fcc4/69++vZcuWydfXV8uXL6+x/QsvvKApU6boV7/6la666io9+eSTGjZsmP7+9783uvjGSnbOvMotGgAAzOJSGCkrK9PevXs1ceLESy/g4aGJEycqKSmpxn2SkpKqtZekyZMn19pekkpLS2Wz2ao9mgNr0gAAYD6Xwsjp06dlt9sVHBxcbXtwcLBycnJq3CcnJ8el9pKUkJCgwMBA5yMiIsKVMuuNaeABADBfq+wosWjRIhUUFDgfWVlZzfI+s+OjNCc+SleHBTbL6wMAgCvr4Erj7t27y9PTU7m5udW25+bmKiQkpMZ9QkJCXGovSVarVVar1ZXSGuS2ERG6bUTzXHUBAAD149KVEW9vbw0fPlwbN250bnM4HNq4caPi4+Nr3Cc+Pr5ae0lKTEystT0AAGhfXLoyIkkLFy7UnDlzNGLECI0aNUrPP/+8iouLdc8990iSZs+erbCwMCUkJEiSHn74YU2YMEFLlizRtGnTtHr1au3Zs0cvv/xy0x4JAABwSy6HkZkzZ+rUqVP63e9+p5ycHA0ZMkQff/yxs5NqZmamPDwuXXAZO3asVq1apd/+9rd67LHH1LdvX61bt05XX3110x0FAABwWxbDMAyzi7gSm82mwMBAFRQUKCAgwOxyAABAPdT3+7tVjqYBAADtB2EEAACYijACAABMRRgBAACmIowAAABTEUYAAICpCCMAAMBUhBEAAGAqwggAADCVy9PBm6FqklibzWZyJQAAoL6qvrevNNm7W4SRwsJCSVJERITJlQAAAFcVFhYqMDCw1p+7xdo0DodDJ0+elL+/vywWS5O9rs1mU0REhLKystrsmjdt/Rg5PvfX1o+R43N/bf0Ym/P4DMNQYWGhevbsWW0R3e9yiysjHh4eCg8Pb7bXDwgIaJO/YN/W1o+R43N/bf0YOT7319aPsbmOr64rIlXowAoAAExFGAEAAKZq12HEarXq8ccfl9VqNbuUZtPWj5Hjc39t/Rg5PvfX1o+xNRyfW3RgBQAAbVe7vjICAADMRxgBAACmIowAAABTEUYAAICp2nwYWbp0qaKjo+Xj46PRo0dr165ddbZfs2aN+vXrJx8fHw0cOFAffvhhC1XquoSEBI0cOVL+/v4KCgrSjBkzlJycXOc+K1eulMViqfbw8fFpoYpd8/vf//6yWvv161fnPu50/iQpOjr6smO0WCyaP39+je1b+/nbtm2bbrnlFvXs2VMWi0Xr1q2r9nPDMPS73/1OoaGh6tixoyZOnKjU1NQrvq6rn+PmUtfxlZeX69FHH9XAgQPl5+ennj17avbs2Tp58mSdr9mQ3/PmdKVzePfdd19W75QpU674uu5wDiXV+Hm0WCxavHhxra/Zms5hfb4XSkpKNH/+fHXr1k2dOnXSD3/4Q+Xm5tb5ug397NZXmw4jb731lhYuXKjHH39c+/bt0+DBgzV58mTl5eXV2H7Hjh2aNWuW5s6dq/3792vGjBmaMWOGDh061MKV18/WrVs1f/58ff7550pMTFR5ebkmTZqk4uLiOvcLCAhQdna285GRkdFCFbtuwIAB1Wrdvn17rW3d7fxJ0u7du6sdX2JioiTptttuq3Wf1nz+iouLNXjwYC1durTGn//lL3/RX//6Vy1btkw7d+6Un5+fJk+erJKSklpf09XPcXOq6/jOnz+vffv26f/+7/+0b98+vffee0pOTtb3v//9K76uK7/nze1K51CSpkyZUq3eN998s87XdJdzKKnacWVnZ2v58uWyWCz64Q9/WOfrtpZzWJ/vhV/84hf673//qzVr1mjr1q06efKkbr311jpftyGfXZcYbdioUaOM+fPnO5/b7XajZ8+eRkJCQo3tb7/9dmPatGnVto0ePdr46U9/2qx1NpW8vDxDkrF169Za26xYscIIDAxsuaIa4fHHHzcGDx5c7/bufv4MwzAefvhho3fv3obD4ajx5+50/iQZa9eudT53OBxGSEiIsXjxYue2/Px8w2q1Gm+++Watr+Pq57ilfPf4arJr1y5DkpGRkVFrG1d/z1tSTcc4Z84cY/r06S69jjufw+nTpxs33HBDnW1a8zn87vdCfn6+4eXlZaxZs8bZ5ptvvjEkGUlJSTW+RkM/u65os1dGysrKtHfvXk2cONG5zcPDQxMnTlRSUlKN+yQlJVVrL0mTJ0+utX1rU1BQIEnq2rVrne2KiooUFRWliIgITZ8+XV999VVLlNcgqamp6tmzp3r16qU777xTmZmZtbZ19/NXVlamN954Q/fee2+dC0K60/n7tvT0dOXk5FQ7R4GBgRo9enSt56ghn+PWpKCgQBaLRZ07d66znSu/563Bli1bFBQUpLi4OD344IM6c+ZMrW3d+Rzm5ubqgw8+0Ny5c6/YtrWew+9+L+zdu1fl5eXVzke/fv0UGRlZ6/loyGfXVW02jJw+fVp2u13BwcHVtgcHBysnJ6fGfXJyclxq35o4HA4tWLBA11xzja6++upa28XFxWn58uV6//339cYbb8jhcGjs2LE6fvx4C1ZbP6NHj9bKlSv18ccf68UXX1R6errGjRunwsLCGtu78/mTpHXr1ik/P1933313rW3c6fx9V9V5cOUcNeRz3FqUlJTo0Ucf1axZs+pcfMzV33OzTZkyRa+//ro2btyoP//5z9q6daumTp0qu91eY3t3Poevvfaa/P39r3gLo7Wew5q+F3JycuTt7X1ZQL7Sd2NVm/ru4yq3WLUXVzZ//nwdOnToivcp4+PjFR8f73w+duxYXXXVVXrppZf05JNPNneZLpk6darzvwcNGqTRo0crKipKb7/9dr3+peJuXn31VU2dOlU9e/astY07nb/2rLy8XLfffrsMw9CLL75YZ1t3+z2/4447nP89cOBADRo0SL1799aWLVt04403mlhZ01u+fLnuvPPOK3YSb63nsL7fC61Bm70y0r17d3l6el7WQzg3N1chISE17hMSEuJS+9bioYce0vr167V582aFh4e7tK+Xl5eGDh2qtLS0Zqqu6XTu3FmxsbG11uqu50+SMjIy9Mknn+i+++5zaT93On9V58GVc9SQz7HZqoJIRkaGEhMTXV6S/Uq/561Nr1691L1791rrdcdzKEmffvqpkpOTXf5MSq3jHNb2vRASEqKysjLl5+dXa3+l78aqNvXdx1VtNox4e3tr+PDh2rhxo3Obw+HQxo0bq/3L8tvi4+OrtZekxMTEWtubzTAMPfTQQ1q7dq02bdqkmJgYl1/Dbrfr4MGDCg0NbYYKm1ZRUZGOHDlSa63udv6+bcWKFQoKCtK0adNc2s+dzl9MTIxCQkKqnSObzaadO3fWeo4a8jk2U1UQSU1N1SeffKJu3bq5/BpX+j1vbY4fP64zZ87UWq+7ncMqr776qoYPH67Bgwe7vK+Z5/BK3wvDhw+Xl5dXtfORnJyszMzMWs9HQz67DSm8zVq9erVhtVqNlStXGl9//bVx//33G507dzZycnIMwzCMn/zkJ8ZvfvMbZ/vPPvvM6NChg/HMM88Y33zzjfH4448bXl5exsGDB806hDo9+OCDRmBgoLFlyxYjOzvb+Th//ryzzXeP8YknnjA2bNhgHDlyxNi7d69xxx13GD4+PsZXX31lxiHU6Ze//KWxZcsWIz093fjss8+MiRMnGt27dzfy8vIMw3D/81fFbrcbkZGRxqOPPnrZz9zt/BUWFhr79+839u/fb0gynn32WWP//v3O0SRPP/200blzZ+P99983vvzyS2P69OlGTEyMceHCBedr3HDDDcbf/vY35/MrfY5by/GVlZUZ3//+943w8HDjwIED1T6TpaWltR7flX7PW1pdx1hYWGg88sgjRlJSkpGenm588sknxrBhw4y+ffsaJSUlztdw13NYpaCgwPD19TVefPHFGl+jNZ/D+nwvPPDAA0ZkZKSxadMmY8+ePUZ8fLwRHx9f7XXi4uKM9957z/m8Pp/dxmjTYcQwDONvf/ubERkZaXh7exujRo0yPv/8c+fPJkyYYMyZM6da+7ffftuIjY01vL29jQEDBhgffPBBC1dcf5JqfKxYscLZ5rvHuGDBAuf/j+DgYOPmm2829u3b1/LF18PMmTON0NBQw9vb2wgLCzNmzpxppKWlOX/u7uevyoYNGwxJRnJy8mU/c7fzt3nz5hp/J6uOweFwGP/3f/9nBAcHG1ar1bjxxhsvO+6oqCjj8ccfr7atrs9xS6rr+NLT02v9TG7evNn5Gt89viv9nre0uo7x/PnzxqRJk4wePXoYXl5eRlRUlDFv3rzLQoW7nsMqL730ktGxY0cjPz+/xtdozeewPt8LFy5cMH72s58ZXbp0MXx9fY0f/OAHRnZ29mWv8+196vPZbQzLxTcFAAAwRZvtMwIAANwDYQQAAJiKMAIAAExFGAEAAKYijAAAAFMRRgAAgKkIIwAAwFSEEQAAYCrCCAAAMBVhBAAAmIowAgAATEUYAQAApvr/z+fnHf52Q44AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "NPV_all =np.zeros((yobs.shape[0], yobs.shape[2]+1))\n",
    "NPV = np.zeros((yobs.shape[2]+1))\n",
    "for i_time in range(yobs.shape[2]):\n",
    "    alpha = np.power((1-0.05),i_time*100/365)\n",
    "    Qpro_w = Yobs_pred[i_time][:,:num_prod]*Q_diff_w+Q_min\n",
    "    Qpro_g = Yobs_pred[i_time][:,num_prod:num_prod*2]*Q_diff_g+Q_min\n",
    "    Qinj_g = ACTIONS[i_time][:, num_prod:]*Q_diff_g+Q_min\n",
    "    lf3toton =0.1167*4.536e-4 # convert lf^3 to ton \n",
    "    pv = np.sum((50-10)*lf3toton*Qinj_g.detach().cpu().numpy().squeeze())- np.sum(10*Qpro_w.detach().cpu().numpy().squeeze()) - np.sum(50*lf3toton*Qpro_g.detach().cpu().numpy().squeeze())\n",
    "    npv += pv*100*alpha\n",
    "    NPV[i_time+1] = npv\n",
    "    \n",
    "    Qw = yobs[:,:num_prod, i_time]*Q_diff_w+Q_min\n",
    "    Qg = yobs[:,num_prod:num_prod*2, i_time]*Q_diff_g+Q_min\n",
    "    Qinj_g = controls[:, num_prod:-1, i_time]*Q_diff_g+Q_min\n",
    "    pv_all = np.sum((50-10)*lf3toton*Qinj_g.detach().cpu().numpy().squeeze())- np.sum(10*Qw.detach().cpu().numpy().squeeze()) - np.sum(50*lf3toton*Qg.detach().cpu().numpy().squeeze())\n",
    "    npv_all += pv_all*100*alpha\n",
    "    NPV_all[:, i_time+1] = npv_all\n",
    "\n",
    "plt.plot(NPV)\n",
    "# for i in range(NPV_all.shape[0]):\n",
    "#     plt.plot(NPV_all[i,:],'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d47144b7-d18a-4596-9797-51a13d788b61",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "NPV.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a646896-7e82-4bc4-867e-a266bb149397",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.plot(np.arange(yobs.shape[2]+1)*100, NPV)\n",
    "# for i in range(NPV_all.shape[0]):\n",
    "#     plt.plot(NPV_all[i,:],'b')\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.xlabel('Days', fontsize=14)\n",
    "plt.ylabel('USD', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "246783f6-34e2-42c2-9a09-fd27b7852bb1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}