socnet_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b6c82652-b78a-494c-9f5e-38b68cbddc12",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random as rd\n",
    "import sys\n",
    "import copy\n",
    "import pickle\n",
    "import pathlib\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cf616090-c8a5-4f6d-81c1-8e1902e5aab4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ffd00329-f04d-4a46-9303-90c1995356d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow.keras.layers import Dense, Activation, Flatten\n",
    "from tensorflow.keras.layers import Conv1D, concatenate\n",
    "from tensorflow.keras.layers import ZeroPadding1D, Reshape, Input, Dropout, PReLU\n",
    "from tensorflow.keras.models import Sequential, Model\n",
    "import code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e2e13c29-8570-4e3f-ab8e-28527f36db1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir   = pathlib.Path('/project/informatics_lab/data_science_cop/socrates_emulation/')\n",
    "output_dir = pathlib.Path(os.environ['SCRATCH']) / 'ml_weather_tutorial'\n",
    "\n",
    "if not output_dir.is_dir():\n",
    "    output_dir.mkdir()\n",
    "    print(f'creating directory {output_dir}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c0213d-221d-4395-ad4f-6f6081326bb9",
   "metadata": {},
   "source": [
    "## Define inputs\n",
    "Specify the hyperparametersfor the pipeline and the locationof the input data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "231016d6-25cc-417a-8e02-1db4df4b1aeb",
   "metadata": {},
   "source": [
    "Set up the hyperparameters for out training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e36f5608-da97-47a5-b67e-269b1a79ffb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "wl='sw'\n",
    "target='nflx'\n",
    "nsamps = '50.0K'\n",
    "scale_data = True\n",
    "if wl=='sw':\n",
    "  model = 'sw_260'\n",
    "  model_ref = 'sw_ga7'\n",
    "elif wl=='lw':\n",
    "  model = 'lw_300'\n",
    "  model_ref = 'lw_ga7'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9237f4a-b90f-4c14-b7af-e9fc8cb6c73d",
   "metadata": {},
   "source": [
    "Construct the paths to file names that contain the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "94df4384-9f1e-4836-a69a-5defff906e5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "root dir: /project/informatics_lab/data_science_cop/socrates_emulation\n"
     ]
    }
   ],
   "source": [
    "fnext='train'\n",
    "fn_meta = model+'_meta_'+nsamps+'_'+fnext+'.npz'\n",
    "fn_dat_levs = model+'_dat_levs_'+nsamps+'_'+fnext+'.npz'\n",
    "fn_dat_lays = model+'_dat_lays_'+nsamps+'_'+fnext+'.npz'\n",
    "fn_dat_surf = model+'_dat_surf_'+nsamps+'_'+fnext+'.npz'\n",
    "if target=='nflx':\n",
    "  fn_trg = model+'_trg_levs_'+nsamps+'_'+fnext+'.npz'\n",
    "if target=='ndiv':\n",
    "  fn_trg = model+'_trg_lays_'+nsamps+'_'+fnext+'.npz'\n",
    " \n",
    "fnext='test'\n",
    "fn_meta_test = model+'_meta_'+nsamps+'_'+fnext+'.npz'\n",
    "fn_dat_levs_test = model+'_dat_levs_'+nsamps+'_'+fnext+'.npz'\n",
    "fn_dat_lays_test = model+'_dat_lays_'+nsamps+'_'+fnext+'.npz'\n",
    "fn_dat_surf_test = model+'_dat_surf_'+nsamps+'_'+fnext+'.npz'\n",
    "if target=='nflx':\n",
    "  fn_trg_test = model+'_trg_levs_'+nsamps+'_'+fnext+'.npz'\n",
    "  fn_trg_ref = model_ref+'_trg_levs_'+nsamps+'_'+fnext+'.npz'\n",
    "if target=='ndiv':\n",
    "  fn_trg_test = model+'_trg_lays_'+nsamps+'_'+fnext+'.npz'\n",
    "  fn_trg_ref = model_ref+'_trg_lays_'+nsamps+'_'+fnext+'.npz'\n",
    "\n",
    "print('root dir:',data_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5e1f5ff-3ebf-4db4-836c-aedf4b2fe328",
   "metadata": {},
   "source": [
    "### Loading and preparing the  training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9cf89ea2-84be-4400-91a9-1511305fcdee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading sw_260_dat_lays_50.0K_train.npz\n"
     ]
    }
   ],
   "source": [
    "print('loading',fn_dat_lays)\n",
    "with np.load(data_dir / fn_dat_lays) as npzfile:\n",
    "    dat_lays = npzfile['dat_lays']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "c8587209-3ae2-4919-85a2-4adb6d9cfbdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading sw_260_dat_surf_50.0K_train.npz\n"
     ]
    }
   ],
   "source": [
    "print('loading',fn_dat_surf)\n",
    "with np.load(data_dir / fn_dat_surf) as npzfile:\n",
    "    dat_surf = npzfile['dat_surf']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "c411a38d-f1d3-493c-89c3-39001d68b671",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading sw_260_trg_levs_50.0K_train.npz\n"
     ]
    }
   ],
   "source": [
    "print('loading',fn_trg)\n",
    "with np.load(data_dir / fn_trg) as npzfile:\n",
    "    if target=='nflx':\n",
    "         trg = npzfile['trg_levs']\n",
    "    elif target=='ndiv':\n",
    "         trg = npzfile['trg_lays']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "64fe3f4f-2e31-485c-97d8-ba62eb36f1c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "nsamps = trg.shape[0]\n",
    "nlays = dat_lays.shape[1]\n",
    "nlay_feats = dat_lays.shape[2]\n",
    "nsurf_feats = dat_surf.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf1ac304-f48f-495b-bbc7-17119771990a",
   "metadata": {},
   "source": [
    "Next we normalise the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "858047d0-1d84-4a1a-b71c-0542f5a02345",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[7.48596189e-02, 2.07117128e-02, 1.28316524e-02, ...,\n",
       "         0.00000000e+00, 9.53195985e-02, 0.00000000e+00],\n",
       "        [7.50777509e-02, 2.07628967e-02, 1.28662172e-02, ...,\n",
       "         0.00000000e+00, 2.39585537e-01, 0.00000000e+00],\n",
       "        [7.52075666e-02, 2.07914892e-02, 1.28858323e-02, ...,\n",
       "         0.00000000e+00, 3.97299512e-01, 0.00000000e+00],\n",
       "        ...,\n",
       "        [2.59644903e-03, 9.05737910e-04, 3.99721354e-03, ...,\n",
       "         7.23902259e-01, 9.38158696e-01, 0.00000000e+00],\n",
       "        [1.61318348e-03, 5.63421798e-04, 2.31647339e-03, ...,\n",
       "         7.23902259e-01, 9.41771698e-01, 0.00000000e+00],\n",
       "        [5.87507001e-04, 2.05322775e-04, 5.60418666e-04, ...,\n",
       "         7.23902259e-01, 9.44402078e-01, 0.00000000e+00]],\n",
       "\n",
       "       [[1.26202793e-04, 3.66322592e-06, 1.57547697e-03, ...,\n",
       "         0.00000000e+00, 4.29455379e-02, 0.00000000e+00],\n",
       "        [1.28851030e-04, 3.75739613e-06, 1.57675368e-03, ...,\n",
       "         0.00000000e+00, 1.72811216e-01, 0.00000000e+00],\n",
       "        [1.42091854e-04, 4.38634354e-06, 1.57934227e-03, ...,\n",
       "         0.00000000e+00, 3.61180706e-01, 0.00000000e+00],\n",
       "        ...,\n",
       "        [3.56097778e-03, 1.24615131e-03, 5.65986030e-03, ...,\n",
       "         2.00382283e-01, 6.01805922e-01, 0.00000000e+00],\n",
       "        [2.79017377e-03, 9.76690095e-04, 4.33783035e-03, ...,\n",
       "         2.00382283e-01, 6.05677306e-01, 0.00000000e+00],\n",
       "        [1.22915369e-03, 4.30274772e-04, 1.66176715e-03, ...,\n",
       "         2.00382283e-01, 6.09110093e-01, 0.00000000e+00]],\n",
       "\n",
       "       [[4.06717973e-02, 2.62244614e-03, 3.78832342e-03, ...,\n",
       "         0.00000000e+00, 6.46323584e-02, 0.00000000e+00],\n",
       "        [4.06380466e-02, 2.62819719e-03, 3.80200744e-03, ...,\n",
       "         0.00000000e+00, 2.13664850e-01, 0.00000000e+00],\n",
       "        [4.05636330e-02, 2.63300742e-03, 3.81482172e-03, ...,\n",
       "         0.00000000e+00, 3.54187027e-01, 0.00000000e+00],\n",
       "        ...,\n",
       "        [3.39010516e-03, 1.18568665e-03, 5.36066482e-03, ...,\n",
       "         6.17225572e-01, 8.74296899e-01, 0.00000000e+00],\n",
       "        [2.26714713e-03, 7.92685944e-04, 3.43576206e-03, ...,\n",
       "         6.17225572e-01, 8.77988308e-01, 0.00000000e+00],\n",
       "        [8.94574833e-04, 3.12762047e-04, 1.08630962e-03, ...,\n",
       "         6.17225572e-01, 8.80861147e-01, 0.00000000e+00]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[8.90418255e-02, 6.84193111e-02, 4.64770676e-02, ...,\n",
       "         0.00000000e+00, 3.58349724e-02, 0.00000000e+00],\n",
       "        [8.98498831e-02, 6.90489234e-02, 4.68789632e-02, ...,\n",
       "         0.00000000e+00, 2.15075214e-01, 0.00000000e+00],\n",
       "        [9.06251868e-02, 6.99099747e-02, 4.73874082e-02, ...,\n",
       "         0.00000000e+00, 3.93528972e-01, 0.00000000e+00],\n",
       "        ...,\n",
       "        [2.64375277e-03, 9.21061639e-04, 4.07203019e-03, ...,\n",
       "         7.20056039e-01, 9.34961854e-01, 0.00000000e+00],\n",
       "        [1.62101484e-03, 5.65673077e-04, 2.32830185e-03, ...,\n",
       "         7.20056039e-01, 9.38546389e-01, 0.00000000e+00],\n",
       "        [5.79207395e-04, 2.02461328e-04, 5.46636043e-04, ...,\n",
       "         7.20056039e-01, 9.41129724e-01, 0.00000000e+00]],\n",
       "\n",
       "       [[1.16998933e-02, 6.83098835e-03, 4.84883525e-02, ...,\n",
       "         0.00000000e+00, 6.19977568e-02, 0.00000000e+00],\n",
       "        [1.22530887e-02, 7.21981659e-03, 4.86410554e-02, ...,\n",
       "         0.00000000e+00, 2.31925277e-01, 0.00000000e+00],\n",
       "        [1.31090457e-02, 8.62433946e-03, 4.95706984e-02, ...,\n",
       "         0.00000000e+00, 3.23150936e-01, 0.00000000e+00],\n",
       "        ...,\n",
       "        [1.34099989e-03, 4.68035714e-04, 1.85051877e-03, ...,\n",
       "         2.09834244e-01, 7.37400872e-01, 0.00000000e+00],\n",
       "        [9.99488045e-04, 3.49430487e-04, 1.27065805e-03, ...,\n",
       "         2.09834244e-01, 7.40680463e-01, 0.00000000e+00],\n",
       "        [4.24838375e-04, 1.48772731e-04, 2.85306984e-04, ...,\n",
       "         2.09834244e-01, 7.43663555e-01, 0.00000000e+00]],\n",
       "\n",
       "       [[1.88133975e-02, 1.53180230e-04, 1.56832894e-03, ...,\n",
       "         0.00000000e+00, 8.68917926e-02, 0.00000000e+00],\n",
       "        [1.86371925e-02, 1.51315885e-04, 1.57021199e-03, ...,\n",
       "         0.00000000e+00, 1.88251397e-01, 0.00000000e+00],\n",
       "        [1.81890123e-02, 1.46603953e-04, 1.57086960e-03, ...,\n",
       "         0.00000000e+00, 3.57811731e-01, 0.00000000e+00],\n",
       "        ...,\n",
       "        [3.56657801e-03, 1.24783482e-03, 5.66784742e-03, ...,\n",
       "         2.90813620e-01, 7.76430273e-01, 0.00000000e+00],\n",
       "        [2.61987641e-03, 9.16881395e-04, 4.04464564e-03, ...,\n",
       "         2.90813620e-01, 7.79778622e-01, 0.00000000e+00],\n",
       "        [1.07235600e-03, 3.75419398e-04, 1.39306943e-03, ...,\n",
       "         2.90813620e-01, 7.82297060e-01, 0.00000000e+00]]])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat_lays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "049c3139-56e2-45fc-a41d-ee5adddaa68e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50000, 70, 29)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dat_lays.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f0cfd63a-2526-4ef0-915a-e1b5d6b12ff7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "removing constant layer features: 6\n"
     ]
    }
   ],
   "source": [
    "if scale_data: # normalize by range\n",
    "  scaler_lays = []\n",
    "  use_lays = []\n",
    "  for ic in range(nlay_feats):\n",
    "    min0 = np.min(dat_lays[:,:,ic])\n",
    "    range0 = np.max(dat_lays[:,:,ic]) - min0\n",
    "    if range0 > 0.:\n",
    "      dat_lays[:,:,ic] = (dat_lays[:,:,ic] - min0)/range0\n",
    "      scaler_lays.append([min0, range0])\n",
    "      use_lays.append(ic)\n",
    "  if len(use_lays)<nlay_feats:\n",
    "    print('removing constant layer features:', nlay_feats-len(use_lays))\n",
    "    dat_lays = dat_lays[:,:,use_lays]\n",
    "    nlay_feats = len(use_lays)\n",
    "      \n",
    "\n",
    "  scaler_surf = []\n",
    "  use_surf = []\n",
    "  for ic in range(nsurf_feats):\n",
    "    min0 = np.min(dat_surf[:,ic])\n",
    "    range0 = np.max(dat_surf[:,ic]) - min0\n",
    "    if range0 > 0.:\n",
    "      dat_surf[:,ic] = (dat_surf[:,ic] - min0)/range0\n",
    "      scaler_surf.append([min0, range0])\n",
    "      use_surf.append(ic)\n",
    "  if len(use_surf)<nsurf_feats:\n",
    "    print('removing constant surf features:', nsurf_feats-len(use_surf))\n",
    "    dat_surf = dat_surf[:,:,use_surf]\n",
    "    nsurf_feats = len(use_surf)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "95e28b0f-98cb-4a67-8223-8baf862eeae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "ntrg_samps = trg.shape[0]\n",
    "\n",
    "if target=='nflx' or target=='ndiv':\n",
    "  nouts=1\n",
    "  ntrg_levs = trg.shape[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "4dd67ee4-3bb4-41d2-9489-99aade718ade",
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_model_mlp(nlays, nlay_feats):\n",
    "    profile_input = Input(shape=(nlays, nlay_feats), name='profile_input')\n",
    "    surf_input = Input(shape=(nsurf_feats,), name='surf_input')\n",
    "    flat_profs = Flatten()(profile_input)\n",
    "    raw_in = tf.keras.layers.concatenate([flat_profs, surf_input])\n",
    "    raw_size = (nlays*nlay_feats)+nsurf_feats\n",
    "    prof_size = nlays*nlay_feats\n",
    "\n",
    "    x = Dense(512, use_bias=False, activation='relu')(raw_in)\n",
    "    x = Dense(512, use_bias=False, activation='relu')(x)\n",
    "    x = Dense(256, use_bias=False, activation='relu')(x)\n",
    "    x = Dense(256, use_bias=False, activation='relu')(x)\n",
    "    x = Dense(128, use_bias=False, activation='relu')(x)\n",
    "    x = Dense(128, use_bias=False, activation='relu')(x)\n",
    "\n",
    "    main_output = Dense(ntrg_levs, use_bias=True, activation='linear', name='main_output')(x)\n",
    "    model = Model(inputs=[profile_input, surf_input], outputs=[main_output])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "7add4bed-8d1f-4df7-865f-7b22304637fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_model_cnn(nlays, nlay_feats):\n",
    "    profile_input = Input(shape=(nlays, nlay_feats), name='profile_input')\n",
    "    surf_input = Input(shape=(nsurf_feats,), name='surf_input')\n",
    "    flat_profs = Flatten()(profile_input)\n",
    "    raw_in = tf.keras.layers.concatenate([flat_profs, surf_input])\n",
    "    raw_size = (nlays*nlay_feats)+nsurf_feats\n",
    "    prof_size = nlays*nlay_feats\n",
    "\n",
    "    out = ZeroPadding1D(padding=1)(profile_input)\n",
    "    out = Conv1D(32, 3, strides=1, activation='relu', use_bias=False, kernel_initializer='glorot_uniform', bias_initializer='zeros')(out)\n",
    "    ident = out\n",
    "    out = ZeroPadding1D(padding=1)(out)\n",
    "    out = Conv1D(32, 3, strides=1, activation='relu', use_bias=False, kernel_initializer='glorot_uniform', bias_initializer='zeros')(out)\n",
    "    out = ZeroPadding1D(padding=1)(out)\n",
    "    out = Conv1D(32, 3, strides=1, activation='relu', use_bias=False, kernel_initializer='glorot_uniform', bias_initializer='zeros')(out)\n",
    "    x = tf.keras.layers.add([out, ident])\n",
    "    out = Flatten()(x)\n",
    "    out = Dense(prof_size, use_bias=False, activation='relu')(out)\n",
    "\n",
    "    out = tf.keras.layers.concatenate([out, surf_input])\n",
    "    x = tf.keras.layers.add([out, raw_in])\n",
    "    x = Dense(1024, use_bias=False, activation='relu')(x)\n",
    "    x = Dense(1024, use_bias=False, activation='relu')(x)\n",
    "\n",
    "    main_output = Dense(ntrg_levs, use_bias=True, activation='linear', name='main_output')(x)\n",
    "    model = Model(inputs=[profile_input, surf_input], outputs=[main_output])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "3f9f6a79-f2dd-469a-bd46-98ef45651290",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_dict = {'mlp': {'build_func': build_model_mlp,},\n",
    "              'cnn_1d': {'build_func': build_model_cnn,},\n",
    "             }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "2eb36ffd-b068-4acd-8f2e-534fccc85bf5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building and training model mlp\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-07-08 15:56:53.774845: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  SSE4.1 SSE4.2 AVX AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building and training model cnn_1d\n",
      "CPU times: user 4min 19s, sys: 28.9 s, total: 4min 48s\n",
      "Wall time: 52.6 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for model_name in model_dict.keys():\n",
    "    print(f'building and training model {model_name}')\n",
    "    model_dict[model_name]['model_object'] = model_dict[model_name]['build_func'](nlays=nlays, nlay_feats=nlay_feats)\n",
    "    model_dict[model_name]['model_object'].compile(loss='mean_absolute_error',\n",
    "                                                   optimizer='adam')\n",
    "    model_dict[model_name]['model_object'].fit([dat_lays, dat_surf], trg, epochs=1, batch_size=32, verbose=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45d1d8ff-fe98-4c19-b177-a96ec53fd3ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "del dat_lays, dat_surf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af4b3b67-d72f-49ea-81e2-d4ce7de5cabb",
   "metadata": {},
   "source": [
    "### Testing the output\n",
    "\n",
    "Next we load in the test data and do inferenceto check the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "347d87fe-a2a7-489f-8d50-b02ac54c681b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading sw_260_dat_lays_50.0K_test.npz\n",
      "loading sw_260_dat_surf_50.0K_test.npz\n",
      "loading sw_260_trg_levs_50.0K_test.npz\n",
      "loading sw_ga7_trg_levs_50.0K_test.npz\n"
     ]
    }
   ],
   "source": [
    "print('loading',fn_dat_lays_test)\n",
    "with np.load(data_dir / fn_dat_lays_test) as npzfile:\n",
    "    dat_lays_test = npzfile['dat_lays']\n",
    "\n",
    "print('loading',fn_dat_surf_test)\n",
    "with np.load(data_dir / fn_dat_surf_test) as npzfile:\n",
    "    dat_surf_test = npzfile['dat_surf']\n",
    "\n",
    "print('loading',fn_trg_test)\n",
    "with np.load(data_dir / fn_trg_test)as npzfile:\n",
    "    if target=='nflx':\n",
    "        trg_test = npzfile['trg_levs']\n",
    "    elif target=='ndiv':\n",
    "         trg_test = npzfile['trg_lays']\n",
    "\n",
    "print('loading',fn_trg_ref)\n",
    "with np.load(data_dir / fn_trg_ref) as npzfile:\n",
    "    if target=='nflx':\n",
    "         trg_ref = npzfile['trg_levs_ref']\n",
    "    elif target=='ndiv':\n",
    "         trg_ref = npzfile['trg_lays_ref']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "0a53579e-00bb-4129-989a-8bbc881c696f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# scale test data\n",
    "if scale_data: # normalize by range\n",
    "  dat_lays_test = dat_lays_test[:,:,use_lays]\n",
    "  for ic in range(nlay_feats):\n",
    "    dat_lays_test[:,:,ic] = (dat_lays_test[:,:,ic] - scaler_lays[ic][0])/scaler_lays[ic][1]\n",
    "\n",
    "  dat_surf_test = dat_surf_test[:,use_surf]\n",
    "  for ic in range(nsurf_feats):\n",
    "    dat_surf_test[:,ic] = (dat_surf_test[:,ic] - scaler_surf[ic][0])/scaler_surf[ic][1]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6429339-8eeb-43a1-8cf4-78d2cad1df74",
   "metadata": {},
   "source": [
    "### Do inference on test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "5831bb99-bd5e-4b7a-9ccf-cfa7aafd71dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "8d3920db-faf8-4c36-8085-57ac6380d0c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "for model_name, selected_model in model_dict.items():\n",
    "  predictions[model_name] = selected_model['model_object'].predict([dat_lays_test, dat_surf_test])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24af6c78-19c6-4052-8957-104c1fb47ba7",
   "metadata": {},
   "source": [
    "### Calaculate metrics for test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "ec68fea0-9cfd-419e-bdc8-3ec2bba88b05",
   "metadata": {},
   "outputs": [],
   "source": [
    "metrics_dict = {}\n",
    "for model_name in model_dict.keys():\n",
    "    metrics_dict['me_p'] = np.zeros(ntrg_levs)\n",
    "    metrics_dict['me_ctl'] = np.zeros(ntrg_levs)\n",
    "    metrics_dict['mae_p'] = np.zeros(ntrg_levs)\n",
    "    metrics_dict['mae_ctl'] = np.zeros(ntrg_levs)\n",
    "    for ilev in range(ntrg_levs):\n",
    "      metrics_dict['me_p'][ilev] = np.mean(predictions[model_name][:,ilev] - trg_test[:,ilev])\n",
    "      metrics_dict['me_ctl'][ilev] = np.mean(trg_ref[:,ilev] - trg_test[:,ilev])\n",
    "      metrics_dict['mae_p'][ilev] = np.mean(np.abs(predictions[model_name][:,ilev] - trg_test[:,ilev]))\n",
    "      metrics_dict['mae_ctl'][ilev] = np.mean(np.abs(trg_ref[:,ilev] - trg_test[:,ilev]))    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21a0ebc4-2f21-4de9-ba5a-86f78e7f4752",
   "metadata": {},
   "source": [
    "### Visualise metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "80f046ee-3343-4f8f-8bce-9ae3a3d627a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "yax = np.arange(1,len(me_p[1:])+1)[::-1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "d013ba7a-0f0f-40f7-866f-c849c26d21ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6wAAAGDCAYAAAAvXp2OAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABVrklEQVR4nO3de5xcdX34/9ebkBBCuCdAmgCBb7kHCJigIF8aBUFpAq0g4A8pKVWw8QIWWsWiUKFesaZVaE2/+kUNkiDwRUClXGqqUQomEBCINBQWEkQugQghpBjy+f0xk3Wz7Myemd2Z89nM6/l4nMfs3M68d9gXm7NzzkyklJAkSZIkKTeblT2AJEmSJEl9cYNVkiRJkpQlN1glSZIkSVlyg1WSJEmSlCU3WCVJkiRJWXKDVZIkSZKUJTdYO0BELIiI9w/Sui6LiOcj4jeDsb4GH3taRKxo9+NKrWajUt5sVMqbjW7a3GDtMBExMyIWNnnfXYHzgf1TSrsM7mSSwEal3NmolDcb3fS4wZqJiNi87BkK2B1YmVJ6ttE7DpHvT6ppiPwM26g61hD5GbZRdawh8jNsoxlyg7VEEdEVER+PiAeAVyJi84h4S0T8PCJWRcT9ETGtx+1nRsRjEfFyRDweEadXL78kIub2uN3EiEi9w4mI/YB/AQ6PiNURsap6+fER8XB1vU9FxAV9zHoMcDvwB9X7XlW9/ISIeKg674LqY9T8/vpYb4qIWRGxrPr4l0bE/4qIuyLipYi4NiJG1Hn+LqzO/mJE/N+IGFn0+Zf6Y6M2qrzZqI0qbzZqo4MipeRS0gJ0AUuAXYEtgfHASuB4Kn9MeEf1/FhgK+AlYJ/qfccBB1S/vgSY22O9E4EEbF49vwB4f/XrmcDCXnM8Dfzv6tfbA4fWmHcasKLH+b2BV6pzDgf+BngUGNHX91djnQm4CdgGOAD4H+BOYE9gW+Bh4Mwaj98FPFhd/w7Az4DLyv7v6rLpLDZqoy55LzZqoy55LzZqo4Ox+Apr+f4ppbQ8pfQq8D7ghymlH6aU1qeUbgcWUYkaYD0wKSK2TCk9nVJ6aJBm+B2wf0Rsk1J6MaV0b8H7nQr8IKV0e0rpd8DlVP5ndESP2/T8/mr5Qkrpper38yBwW0rpsZTSb4EfAYfUue/Xqut/Afh74L0FZ5eKslEbVd5s1EaVNxu10QFxg7V8y3t8vTvwnuouB6uquzEcCYxLKb1CJZoPAk9HxA8iYt9BmuEkKv+jeCIi/iMiDi94vz8AnthwJqW0nsr3M77HbZb3vlMfnunx9at9nB9d57491/9EdSZpMNmojSpvNmqjypuN2uiAuMFavtTj6+XAd1JK2/VYtkopfR4gpfRvKaV3UNlF4lfAv1bv9wowqsd66r2rWXrDBSn9IqV0IrATcCNwbcHZf03lfzwARERQ2WXhqXqPN8h27fH1btWZpMFkowNjo2o1Gx0YG1Wr2ejAdHyjbrDmZS4wIyKOi4hhETEyKp/HNCEidq4e9L0VlX3fVwOvV++3BDgqInaLiG2BC+s8xjPAhA0Hd0fEiIg4PSK2re7q8FKP9fbnWuCPI+LoiBhO5W3A/wf4eYPf90B8qPr87AB8EpjfxsdW57HRxtmo2slGG2ejaicbbVzHN+oGa0ZSSsuBE6n8MD5H5a9Qf03lv9NmVCL5NfAC8EfArOr9bqfyw/sAsBi4pc7D/DvwEPCbiHi+etkZQFdEvERlN4z3FZz3keptvwo8D8wAZqSUXiv2HQ+K7wK3AY9Vl8va+NjqMDbaFBtV29hoU2xUbWOjTen4RiOlVr+KLbVGRHRReUe4O8qeRdIb2aiUNxuV8majFb7CKkmSJEnKkhuskiRJkqQsuUuwJEmSJClLvsIqSZIkScqSG6ySJEmSpCxtXvYARYwZMyZNnDix7DGkUi1evPj5lNLYsufoi41KNirlzkalvNVqdEhssE6cOJFFixaVPYZUqoh4ouwZarFRyUal3NmolLdajbpLsCRJkiQpS26wSpIkSZKy5AarJEmSJClLQ+IYVpXsvPMqp7NnN72K3/3ud6xYsYK1a9cOykibspEjRzJhwgSGDx9e9ijaxNjh4LBRNeq8ZcsAmL3XXoOyPluuz0bVDnbYvEYbdYNV/VuyZMCrWLFiBVtvvTUTJ04kIgY+0yYqpcTKlStZsWIFe+yxR9njaKiYNq1yumBB3ZvZ4cDZqJqxZPXqQV2fLddmo2rGtPvuA2DBIYcUvo8dNqeZRt0lWG2xdu1adtxxR4PuR0Sw4447+tc6tYQdDpyNKge2XJuNql3ssDnNNOoGq9rGoIvxeVIr+fM1cD6HyoE/h7X53Khd/FlrTqPPmxus6hgRwRlnnNF9ft26dYwdO5bp06cDcNVVVzF27FgmT57cvTz88MNljSttkvrrcIMTTzyRww8/fKPLLrnkEsaPH79Ro6tWrWrH2JJ6Ket36sc+9rHu9e29995st912A16nNFTl+Dv11FNP7V7fxIkTmTx58oDX6TGs6hhbbbUVDz74IK+++ipbbrklt99+O+PHj9/oNqeeeipf+9rXSppQ2vQV6XDVqlXce++9jB49mscff3yjY1w+9rGPccEFF7R7bEm9lPU79Stf+Ur311/96le5r3rsodSJcvydOn/+/O6vzz//fLbddtsBr9NXWNVR3vWud/GDH/wAgGuuuYb3vve9Da9j7ty5HHbYYUyePJlzzjmH119/HYDRo0fz8Y9/nDe96U0cc8wx3HPPPUybNo0999yTm266Caj8xfnDH/5w97qmT5/Ogn7eKEfq1ymnVJYhor8Or7/+embMmMFpp53GvHnzGl7/l770JaZOncpBBx3ExRdfDEBXVxf77rsv73//+5k0aRKnn346d9xxB29961vZa6+9uOeee4DKX5wvv/zy7nVNmjSJrq6uJr9TadM2GL9TL730Uvbdd1/e8Y538N73vre7v3/9139l6tSpHHzwwZx00kmsWbPmDfdt9jGl3k7ZaSdO2WmnssdoSit/p65fv55Zs2ZxwAEHMH36dI4//niuu+46AD7zmc8wdepUJk2axNlnn01KaaP7ppS49tprB6VRX2FV+5133qC88/BGJk8u9LE7p512Gp/5zGeYPn06DzzwAGeddRY//elPu6+fP38+Cxcu7D5/1113seWWW3afX7p0KfPnz+dnP/sZw4cPZ9asWVx99dX82Z/9Ga+88grTpk3jC1/4An/6p3/KRRddxO23387DDz/MmWeeyQknnDCY37H0e7NmNXe/De8u3NMpp1TWt2YNHH/8G6+fObOyPP88nHzyxtcV/ONLfx1ec801XHzxxey8886cfPLJXHjhhd3XfeUrX2Hu3LkAbL/99vz4xz/eaN233XYby5Yt45577iGlxAknnMBPfvITdtttNx599FG+973vMWfOHKZOncp3v/tdFi5cyE033cRnP/tZbrzxxkLzSzk5b9myQX8X4smjRxf6CJ6B/k5dtGgR119/Pffddx/r1q3j0EMP5U1vehMA7373u/nABz4AwEUXXcQ3vvENPvKRj3Tf94knnuDxxx/n7W9/+4C/X2lWr1clmzGtj1f7T9lpJ2aNH8+a11/n+AceeMP1M3fZhZnjxvH8a69x8kMPbXRd0XcsbuXv1BtuuIGuri5++ctf8uyzz7Lffvtx1llnAfDhD3+YT3/60wCcccYZ3HLLLcyYMaP7vj/96U/Zeeed2WsQPs7LDVZ1lIMOOoiuri6uueYaju/jH+P97b505513snjxYqZOnQrAq6++yk7Vv8iNGDGCd77znQAceOCBbLHFFgwfPpwDDzzQV2jUWhteeRg1qtw5CqrX4TPPPMOjjz7KkUceSUSw+eab8+CDDzJp0iSg/92XbrvtNm677TYOqf6iX716NcuWLWO33XZjjz324MADDwTggAMO4OijjyYibFRq0kB/py5cuJATTzyxeyO25z92H3zwQS666CJWrVrF6tWrOe644za677x58zj55JMZNmzYIH036mRrqnvLjRqCP0+t/J26cOFC3vOe97DZZpuxyy678La3va37uh//+Md88YtfZM2aNbzwwgsccMABGzU8mHtAFN5gjYitUkqvDMqjqrMVeCW0lU444QQuuOACFixYwMqVKxu6b0qJM888k8997nNvuG748OHd73q22WabscUWW3R/vW7dOgA233xz1q9f332fwXzbfRvtYBt+QTW6e3m9248aVf/6MWMaf7weanU4f/58Xnzxxe5jbF566SXmzZvHZZddVmi9KSUuvPBCzjnnnI0u7+rq6m4SbFSbjiKvhLbSQH+n1jJz5kxuvPFGDj74YK666qo3HD4zb948rrjiimZGrslGO9eGVz8b+RzW3urdd9SwYXWvHzNixIAeu5W/U/uydu1aZs2axaJFi9h111255JJLNvp9uW7dOm644QYWL17c9PfUU7/HsEbEERHxMLC0ev7giLhyUB5dKsFZZ53Fpz/96e5XWhpx9NFHc9111/Hss88C8MILL/DEE08Uvv/EiRNZsmQJ69evZ/ny5d3HzQ2EjWooqtXhNddcw6233kpXVxddXV0sXry4oWNujjvuOL75zW+yurqL5FNPPdXdaxETJ07k3nvvBeDee+/l8ccfL3zfWmxUm7KB/E498sgjufnmm1m7di2rV6/uPg4P4OWXX2bcuHH87ne/4+qrr97ofo888ggvvvjiG971tFk2qqGuVb9TjzzySK6//nrWr1/PM8880/2How0bp2PGjGH16tXdx7VucMcdd7DvvvsyYcKEgX1jVUVeYf0KcBxwE0BK6f6IOGpQHl0qwYQJEzj33HP7vK738TZXXnklRxxxRPf5/fffn8suu4xjjz2W9evXM3z4cK644gp23333Qo/91re+tXu3xEmTJnHooYcO7JupsFENOX112NXVxZNPPslb3vKW7sv22GMPttlmG+6++25g4+NtAG688UYmTpzYff7YY49l6dKl3f+QHT16NHPnzi282+BJJ53Et7/9bSZPnszUqVPZe++9m/0We7JRbbIG8jt16tSpnHDCCRx88MHsvvvuTJkypfsdRS+99FLe/OY3s/vuu3PggQfy8ssvd9/vmmuu4bTTThvMz8C0UQ1prfqdetJJJ3HnnXcyadIk9t57b9785jez7bbbst122/GBD3yAAw88kIkTJ3YfKrfBvHnzBvUN0aLe7hgAEXF3SunNEXFfSumQ6mX3p5QOHrQp+jFlypS0aNGidj2cetvwxiwD2P1v6dKl7LfffoMyTifo6/mKiMUppSm9b2ujKtqoHQ4eG1UjNrwZy0B2+etpU2p59erVjB49mjVr1nDUUUcxZ86cQfljro2qEc00uil1WM+GRleuXMlhhx3Gz372M3bZZZcBr7eRRou8wro8Io4AUkSMAD5KdZcJSVmwUSlvNirVcPbZZ/Pwww+zdu1azjzzzMHa86hRNirVMH36dFatWsVrr73Gpz71qUHZWG1UkQ3WDwL/CIwHVgC3AR9q5VCSGmKjnW7mzLInUH02KtXw3e9+t+wRwEY73swSNsKGit5veFaGfjdYU0rPA6e3YRZJTbBRucGaNxuV8majmjluXNkjqI4i7xL8rYjYrsf57SPimy2dSpuk/o6XVkWjz5ONiuefrywF2OHA2ahyYMu12aga9fxrr/H8a681fD87bE6jz1u/G6zAQSmlVT0e4EWg0BHJEbFdRFwXEb+KiKURcXhE7BARt0fEsurp9g1NrCFp5MiRrFy50rD7kVJi5cqVjBw5spG72WinO/nkytIPOxw4G1UObLk2G1UzTn7oIU5+6KGG7mOHzWmm0SLHsG4WEdtX4yUidih4P6gcD3BrSunk6kHso4BPAnemlD4fEZ8APgF8vPDEGpImTJjAihUreO6558oeJXsjR45s9HOrbFSF2OHgsFGVzZbrs1G1gx02r9FGi8T4ZeDnEbHhE2HfA/x9f3eKiG2Ao4CZACml14DXIuJEYFr1Zt8CFmDEm7zhw4ezxx57lD3GpspGVYgdlsZGNahsedDZqBpmh+1T5E2Xvh0Ri4G3AQG8O6X0cIF17wk8B/zfiDgYWAycC+ycUnq6uu6nI2KnpqeXZKNS5mxUypuNSnkrcgwrwK+AG4DvA6sjYrcC99kcOBT45+qHML9CZZeIQiLi7IhYFBGLfKld6peNSnmzUSlvNiplqt9XWCPiI8DFwDPA61T+8pSAg/q56wpgRUrp7ur566hE/ExEjKv+xWkc8Gxfd04pzQHmAEyZMsWjmaUabFT85V+WPYHqsFEpbzaqvxw/vuwRVEeRY1jPBfZJKa1sZMUppd9ExPKI2Cel9AhwNPBwdTkT+Hz19PsNzixpYzba6U49tewJVJ+NSnmz0Q536k7utZ2zIhusy4HfNrn+jwBXV9817THgz6nshnxtRPwF8CSVA9slNc9GO93y5ZXTXXctdw7VYqNS3my0wy1fuxaAXRv7OCS1SZEN1seABRHxA+B/NlyYUvqH/u6YUloCTOnjqqOLDiipXzba6c44o3K6YEGpY6gmG5XyZqMd7oylSwFYcEihj99VmxXZYH2yuoyoLpLyYqNS3mxUypuNShkr8rE2fwcQEVullF5p/UiSGmGjUt5sVMqbjUp56/djbSLi8Ih4GFhaPX9wRFzZ8skkFWKjUt5sVMqbjUp5K/I5rLOB44CVACml+4GjWjiTpMbMxkalnM3GRqWczcZGpWwVOYaVlNLyiOh50eutGUdSM2y0w51/ftkTqB82KuXNRjvb+b7LftYKfaxNRBwBpOpbdn+U6i4TkrJgo51uxoyyJ1B9NirlzUY73IwxY8oeQXUU2SX4g8CHgPHACmBy9bykPNhop3vkkcqiXNmolDcb7XCPrFnDI2vWlD2Gaqj7CmtEDANmp5ROb9M8khpgowLgnHMqp34Oa3ZsVMqbjQrgnOofff0c1jzVfYU1pfQ6MLa6e4SkzNiolDcblfJmo1L+ihzD2gX8LCJuAro/myql9A+tGkpSQ7qwUSlnXdiolLMubFTKVpEN1l9Xl82ArVs7jqQm2KiUNxuV8majUsb63WBNKf0dQERslVJ6pb/bS2ovG5XyZqNS3mxUylu/7xIcEYdHxMNU3947Ig6OiCtbPpmkQmxUXHRRZVGWbFTKm43qot1356Lddy97DNVQ5GNtZgPHASsBUkr3A0e1cCZJjZmNjXa2Y46pLMrVbGxUytlsbLSjHbPDDhyzww5lj6EaimywklJa3uui11swi6Qm2WiHW7KksihbNirlzUY725KXX2bJyy+XPYZqKPKmS8sj4gggVd/y+6NUd5mQlAUb7XTnnVc59XNYc2WjUt5stMOd9+ijgJ/Dmqsir7B+EPgQMB5YAUyunpeUBxuV8majUt5sVMpYzVdYI+ILKaWPA29LKZ3expkkFWCjUt5sVMqbjUpDQ71XWI+PiOHAhe0aRlJDbFTKm41KebNRaQiodwzrrcDzwFYR8RIQQNpwmlLapg3zSarNRqW82aiUNxuVhoB6G6wXpZT+OiK+n1I6sW0TSSrKRlXx2c+WPYH6ZqNS3mxUAHx2zz3LHkF11NtgvQs4FHipTbNIaoyNquKII8qeQH2zUSlvNioAjth227JHUB31NlhHRMSZwBER8e7eV6aUbmjdWJIKsFFV/PznlVM3XHNjo1LebFQA/Py3vwXccM1VvQ3WDwKnA9sBM3pdlwAjlsplo6r45Ccrp34Oa25sVMqbjQqATz72GODnsOaq5gZrSmkhsDAiFqWUvtHGmSQVYKNS3mxUypuNSkNDvc9hfXtK6d+BF91NQsqPjUp5s1EpbzYqDQ31dgn+I+DfeeMuEuBuElIObFTKm41KebNRaQiot0vwxdXTP2/fOJKKslEpbzYq5c1GpaGh3i7Bf1Xvjimlfxj8cSQVZaPqNnt22ROoDzYq5c1GtcHsP/zDskdQHfV2Cd66eroPMBW4qXp+BvCTVg4lqRAbVcXkyWVPoL7ZqJQ3GxUAk7feuv8bqTT1dgn+O4CIuA04NKX0cvX8JcD32jKdpJpsVN3uuKNyeswx5c6hjdiolDcb1QZ3vPACAMfssEPJk6gvmxW4zW7Aaz3OvwZMLPoAETEsIu6LiFuq53eIiNsjYln1dPuGJpbUm412ussuqyzKVdON2qfUFjba4S574gkue+KJssdQDUU2WL8D3BMRl0TExcDdwLcaeIxzgaU9zn8CuDOltBdwZ/W8pObZqJS3gTRqn1Lr2aiUsX43WFNKfw/8OfAisAr485TS54qsPCImAH8M/J8eF5/I7/8n8C3gT4qPK6k3G5Xy1myj9im1h41Keav3pkvdUkr3Avc2sf7ZwN/w+4PaAXZOKT1dXe/TEbFTX3eMiLOBswF22223Jh5a6hw2KuWtyUZn02SfYKNSI2xUyleRXYKbEhHTgWdTSoubuX9KaU5KaUpKacrYsWMHeTpJNirla6B9go1KrWSjUvsUeoW1SW8FToiI44GRwDYRMRd4JiLGVf/qNA54toUzSKrNRjcVX/962RNo8NmnlDcb3YR8fZ99yh5BdbTsFdaU0oUppQkppYnAacC/p5TeR+Uzrs6s3uxM4PutmkFSbTa6Cdlnn8qiTYZ9Snmz0U3LPqNGsc+oUWWPoRqa2mCNiDkDeMzPA++IiGXAO6rnJQ0iG+0wN99cWTRkDKBR+5TawEY7y83PP8/Nzz9f9hiqod9dgiNiWErp9V4XN7T/WUppAbCg+vVK4OhG7i+pNhsVX/5y5XTGjHLnUJ8G2qh9Sq1lo/ry8uUAzBgzpuRJ1Jcir7A+HhFzIuLoiAiAgRxgLmnQ2aiUNxuV8majUsaKbLDuA9wBfIhK0F+LiCNbO5akBtiolDcblfJmo1LG+t1gTSm9mlK6NqX0buAQYBvgP1o+maRCbFTKm41KebNRKW+F3nQpIv4oIq6k8oHKI4FTWjqVpIbYqJQ3G5XyZqNSvoq86dLjwBLgWuCvU0qvtHooScXZqPjOd8qeQHXYqJQ3G9V39tuv7BFUR78brMDBKaWXWj6JpGbZaKfbddeyJ1B9NirlzUY73K4jR5Y9guqoucEaEX+TUvoicFn1DdM2klL6aCsHk1Sfjarb/PmV01NPLXcObcRGpbzZqDaY/+yzAJy6004lT6K+1HuFdWn11Lf1lvJko6r453+unLrBmhsblfJmowLgn596CnCDNVc1N1hTSjdXT7/VvnEkFWWjUt5sVMqbjUpDQ71dgm8GUq3rU0ontGQiSYXYqJQ3G5XyZqPS0FBvl+DLq6fvBnYB5lbPvxfoauFMkoqxUSlvNirlzUalIaDeLsH/ARARl6aUjupx1c0R8ZOWTyapLhuV8majUt5sVBoainyszdiI2DOl9BhAROwBjG3tWJIaYKOd7rrryp5A9dmolDcb7XDXHXBA2SOojiIbrB8DFkTEY9XzE4GzWzaRpEbZaKcbM6bsCVSfjUp5s9EON2bEiLJHUB39brCmlG6NiL2AfasX/Sql9D+tHUtSUTYqrrqqcjpzZplTqAYblfJmo7rq6acBmDluXMmTqC9FXmGlGu39LZ5FUpNstMO5wZo9G5XyZqOd7arf/AZwgzVXm5U9gCRJkiRJfXGDVZIkSZKUpX43WCPiL3qdHxYRF7duJEmNsFEpbzYq5c1GpbwVeYX16Ij4YUSMi4hJwH8CW7d4LknF2aiUNxuV8majUsaKvEvw/xcRpwK/BNYA700p/azlk0kqxEbFD39Y9gSqw0alvNmofnjQQWWPoDqK7BK8F3AucD3QBZwREaNaPJekgmxUjBpVWZQlG5XyZqMaNWwYo4YNK3sM1VBkl+CbgU+nlM4B/ghYBvyipVNJaoSNdrorr6wsypWNSnmz0Q535VNPceVTT5U9hmoo8jmsh6WUXgJIKSXgyxFxU2vHktQAG+10115bOZ01q9w5VIuNSnmz0Q537bPPAjBr/PiSJ1Ffimyw/klE9HX5skGeRVJzbFTKm41KebNRKWNFNlin9vh6JHA0cC/w7ZZMJKlRNirlzUalvNmolLEi7xL8kZ7nI2Jb4Dstm0hSQ2xUypuNSnmzUSlvRd50qbc1wF6DPYikQWOjUt5sVMqbjUoZ6fcV1oi4GUjVs5sB+wPXtnIoScXZqFiwoOwJVIeNSnmzUS045JCyR1AdRY5hvbzH1+uAJ1JKK1o0j6TG2aiUNxuV8majUsaKHMP6H82sOCJ2pXKw+i7AemBOSukfI2IHYD4wkcqHM5+SUnqxmceQZKMCLq/+W+uCC8qdQ32yUSlvNqrLn3wSgAt2263kSdSXmsewRsTLEfFSj9OXep4vsO51wPkppf2AtwAfioj9gU8Ad6aU9gLurJ6X1CAbVbdbbqksyoqNSnmzUW1wy8qV3LJyZdljqIZ6r7AenFJ6rNkVp5SeBp6ufv1yRCwFxgMnAtOqN/sWsAD4eLOPI3UwG5XyZqNS3mxUGgLqvUvw9wAi4s6BPkhETAQOAe4Gdq4GviH0nWrc5+yIWBQRi5577rmBjiBtimxUypuNSnmzUWkIqPcK62YRcTGwd0T8Ve8rU0r/UOQBImI0cD1wXkrppYgoNFhKaQ4wB2DKlCmpn5tLnchGpbzZqJQ3G5WGgHqvsJ4GrKWyUbt1H0u/ImI4lYCvTindUL34mYgYV71+HPBsc6NLHc9GVbHllpVFubFRKW82KgC2HDaMLYcNK3sM1VDzFdaU0iPAFyLigZTSjxpdcVT+vPQNYGmvv1DdBJwJfL56+v1G1y3JRtXDjxr+z682sFEpbzaqDX500EFlj6A6inysTbP/EnorcAbwy4hYUr3sk1TivTYi/gJ4EnhPk+uXhI1KubNRKW82KuWt3w3WZqWUFgK1duI/ulWPK6kYG92EXHpp5fRTnyp3Dg0qG5XyZqObjku7ugD41MSJpc6hvtX7HNZx7RxEUmNsVN3uvLOyKCs2KuXNRrXBnS++yJ0vvlj2GKqh3ius34yI7al8dtStwMKU0rq2TCWpCBuV8majUt5sVBoC6r3p0rsiYiSVDz7+U+DyiHiSStC3ppSebM+Ikvpio1LebFTKm41KQ0PdY1hTSmupRgsQEXsA7wK+FhG7pJQOa/2IkmqxUSlvNirlzUal/DX0pksppceBK4ErI2JEa0aS1Cwb7VA77lj2BCrIRqW82Whn2nH48LJHUB1Nv0twSum1wRxE0uCy0Q5y/fVlT6Am2KiUNxvtHNdPmlT2CKqj5rsES5IkSZJUpn43WCNipz4u26c140hqlI2KCy+sLMqSjUp5s1Fd+NhjXPjYY2WPoRqKvML604g4ZcOZiDgf+H+tG0lSg2y00911V2VRrmxUypuNdri7fvtb7vrtb8seQzUUOYZ1GjAnIt4D7AwsBXzHNCkf07BRKWfTsFEpZ9OwUSlb/b7CmlJ6mspbfR8OTAS+nVJa3eK5JBVko1LebFTKm41Keev3FdaIuB14GpgETAC+GRE/SSld0OrhJPXPRqW82aiUNxuV8lbkGNYrUkp/llJalVJ6EDgCcCdvKR822ukmTKgsypWNSnmz0Q43YYstmLDFFmWPoRr6fYU1pXRjr/PrgEtbNZCkxtiomDu37AlUh41KebNRzd1//7JHUB1Fdgl+GUjVsyOA4cDqlNK2rRxMUjE2KuXNRqW82aiUtyKvsG7d83xE/Am+c5qUDRsV551XOZ09u8wpVIONSnmzUZ23bBkAs/faq+RJ1Jcix7BupLrbxNsHfxRJg8FGO9CSJZVFQ4KNSnmz0c6zZPVqlqz2jaFzVWSX4Hf3OLsZMIXf7zYhqWQ2KuXNRqW82aiUt343WIEZPb5eB3QBJ7ZkGknNsFEpbzYq5c1GpYwVOYb1z9sxiKTm2KiUNxuV8majUt5qbrBGxFepsztESumjLZlIUiE2qm577132BOqDjUp5s1FtsPeoUWWPoDrqvcK6qG1TSGqGjapizpyyJ1DfbFTKm40KgDn77FP2CKqj3gbrMSmlMyLi3JTSP7ZtIklF2aiUNxuV8maj0hBQ72Nt3hQRuwNnRcT2EbFDz6VdA0qqyUZVcfbZlUW5sVEpbzYqAM5+5BHOfuSRssdQDfVeYf0X4FZgT2AxED2uS9XLJZXHRlXxX/9V9gTqm41KebNRAfBfa9aUPYLqqPkKa0rpn1JK+wHfTCntmVLao8diwFLJbFTKm41KebNRaWiot0swACmlv2zHIJKaY6NS3mxUypuNSnnrd4NVkiRJkqQy1DuGVZI0FEyeXPYEkiQNWZNHjy57BNXhBqskDXWzZ5c9gSRJQ9bsvfYqewTVUcouwRHxzoh4JCIejYhPlDGDpNpsVMqbjUp5s1Fp8LT9FdaIGAZcAbwDWAH8IiJuSik93O5ZVJC7G3YUG5XyZqNDj7sbdhYblQZXGbsEHwY8mlJ6DCAi5gEnAkacK3c37DQ2KuXNRocYdzfsODYqDaIydgkeDyzvcX5F9bKNRMTZEbEoIhY999xzbRtOko1KmbNRKW82Kg2iMjZYo4/L0hsuSGlOSmlKSmnK2LFj2zCWpCoblfJmo1LebFQaRGVssK4Adu1xfgLw6xLmkNQ3G5XyZqNS3mxUGkRlbLD+AtgrIvaIiBHAacBNJcwhqW82KuXNRqW82ag0iNr+pksppXUR8WHg34BhwDdTSg+1ew5JfbNRKW82KuXNRqXBFSm9YZf67ETEc8ATbXioMcDzbXiconKax1lqa9c8u6eUsjzIxUaz4Cx9a+csNprXf3vIax5n6ZuN0rGN5jQL5DVPp87SZ6NDYoO1XSJiUUppStlzbJDTPM5SW27zbMpye65zmsdZ+pbTLJ0gt+c7p3mcpW85zdIJcnq+c5oF8prHWTZWxjGskiRJkiT1yw1WSZIkSVKW3GDd2JyyB+glp3mcpbbc5tmU5fZc5zSPs/Qtp1k6QW7Pd07zOEvfcpqlE+T0fOc0C+Q1j7P04DGskiRJkqQs+QqrJEmSJClLHbfBGhE7RMTtEbGserp9jdu9MyIeiYhHI+ITfVx/QUSkiBhT1iwR8aWI+FVEPBAR/y8itmtihv6+z4iIf6pe/0BEHFr0vs1odp6I2DUifhwRSyPioYg4t6xZelw/LCLui4hbBjpLJ7HRYuvucb2N2mhb2Wixdfe43kZttK1stNi6e1zftkZz6nMg8/S4vj2NppQ6agG+CHyi+vUngC/0cZthwH8DewIjgPuB/XtcvyuVD4N+AhhT1izAscDm1a+/0Nf9+3n8ut9n9TbHAz8CAngLcHfR+zbxfAxknnHAodWvtwb+ayDzDGSWHtf/FfBd4Jayf+6H0mKjg9aEjdpoSxYbHbQmbNRGW7LY6KA1MaiN5tTnQOfpcX1bGu24V1iBE4FvVb/+FvAnfdzmMODRlNJjKaXXgHnV+23wFeBvgIEeADygWVJKt6WU1lVv95/AhAYfv7/vc8OM304V/wlsFxHjCt63UU3Pk1J6OqV0L0BK6WVgKTC+jFkAImIC8MfA/xnADJ3KRgusu9eMNmqj7WSjBdbda0YbtdF2stEC6+41YzsazanPAc0D7W20EzdYd04pPQ1QPd2pj9uMB5b3OL+iehkRcQLwVErp/rJn6eUsKn8BaUSRdde6TdG52jVPt4iYCBwC3F3iLLOp/I9+/QBm6FQ22ti6bbS5WWZjo82y0cbWbaPNzTIbG22WjTa27nY1mlOfgzHPbNrU6OatfoAyRMQdwC59XPW3RVfRx2UpIkZV13Fs2bP0eoy/BdYBVxedq+i669ymyH0bNZB5KldGjAauB85LKb1UxiwRMR14NqW0OCKmDWCGTZaNFmajLZjFRvtno4XZaAtmsdH+2WhhOTWaU58DmqfdjW6SG6wppWNqXRcRz2x4ab36kvazfdxsBZV99zeYAPwa+F/AHsD9EbHh8nsj4rCU0m/aPMuGdZwJTAeOTik1GlHddfdzmxEF7tuogcxDRAynEvHVKaUbSpzlZOCEiDgeGAlsExFzU0rvG+BMmwwbLcxGWzOLjfbDRguz0dbMYqP9sNHCcmo0pz4HOk97G00tPEA2xwX4Ehsf/P3FPm6zOfAYlWA3HIR8QB+362JgB6IPaBbgncDDwNgmH7/f75PKvuk9D7a+p5HnqI3zBPBtYPYg/Zw0PUuv20zDN4to9Lm30UH4ObRRG23VYqOD83NoozbaqsVGB+fncLAbzanPgc7T6zYtb7T0qNq9ADsCdwLLqqc7VC//A+CHPW53PJV34Ppv4G9rrGugEQ9oFuBRKvuVL6ku/9LEDG9YN/BB4IPVrwO4onr9L4EpjTxH7ZoHOJLKbgwP9Hg+ji9jll7raHnEm9pio4P3c2ijNtqKxUYH7+fQRm20FYuNDt7P4WA3mlOfA31ueqyj5Y1G9YE0BEXEAmBuSmnA784VEZdR+QFdl1Lq6ziETUZEzATen1I6suxZtGmz0ebYqNrFRptjo2oXG23OptZoJ75L8CYpImZGxMIm77srcD6Vz14acgFHxIiIuC4iuqLyAdfTyp5J6s1GbVR5s1EbVd5stHMbdYO1xSJiKLyx1e7AypRSXwfC15XR97cQeB/Q5xsCSLVk9DNcj42qY2X0M1yPjapjZfQzXI+NDmFusLZA9a8fH4+IB4BXImLziHhLRPw8IlZFxP09/zJS/YvRYxHxckQ8HhGnVy+/JCLm9rjdxOpfVTbv9Xj7Af8CHB4RqyNiVfXy4yPi4ep6n4qIC/qY9RjgduAPqve9qnr5CRHxUHXeBdXHqPn99bHeAyLi9oh4ofoOcZ/s8T1dGxHfrs71UERM6bXuCyLigYj4bUTMj4iR9Z7vlNJrKaXZKaWFwOt9zLJjRNwUES9FxD1U3gFPHcxGbVR5s1EbVd5s1EbbqpUHyHbqQuUA9SVU3gZ6SyofsLuSyoHNmwHvqJ4fC2wFvATsU73vOH7/zmiXUNlvf8N6J1I54Hrz6vkFVPZPB5gJLOw1x9PA/65+vT1waI15pwErepzfG3ilOudwKh8K/Cgwoq/vr4/1bV197POpvNX11sCbe3xPa6vPxTDgc8B/9nru7qFyMP4OwFKqB34XfO5XANN6XTYPuLb6XE8Cnur9XLl01mKjNuqS92KjNuqS92KjNtrOxVdYW+efUkrLU0qvUnn5/ocppR+mlNanlG4HFlH5QQZYD0yKiC1TSk+nlB4apBl+B+wfEduklF5MKd1b8H6nAj9IKd2eUvodcDmV/xkd0eM2Pb+/3qYDv0kpfTmltDal9HJK6e4e1y+sPhevA98BDu51/39KKf06pfQCcDMwueDcbxARw4CTgE+nlF5JKT0IfKvZ9WmTYqM2qrzZqI0qbzZqo23hBmvrLO/x9e7Ae6q7HKyq7sZwJDAupfQKlWg+CDwdET+IiH0HaYaTqPyP4omI+I+IOLznlRGxusb9/gB4YsOZlNL66vczvsdtlve+U4917krl7a9r6bnv/RpgZK9dLXpfP7rOuvozlsrnTPWc94kat1VnsdHabFQ5sNHabFQ5sNHabHQQucHaOj0/L2g58J2U0nY9lq1SSp8HSCn9W0rpHVR2kfgV8K/V+70CjOqxnnrvavaGzydKKf0ipXQisBNwI5VdBYr4NZX/8QAQEUElzKfqPV4Py8ln3/nngHVU5t9gt5JmUV5sNA82qlpsNA82qlpsNA+bfKNusLbHXGBGRBwXEcMiYmRETIuICRGxc/Wg762A/wFW8/uDqZcAR0XEbhGxLXBhncd4BpgQESOg++2vT4+Ibau7OrxEHwdp9zA6In5RPbh8Z+CPI2JeRHyYyv75/wMcHRHnV2//ng23j4i/67WuW4BdIuK8iNgiIraOiDcXfK6aUn2cDQesj6g+x1HdFeMG4JKIGBUR+wNntnIWDUk2aqPKm43aqPJmozbaMm6wtkFKaTlwIvBJKn8FWQ78NZXnfzMqkfwaeAH4I2BW9X63A/OBB4DFVOKo5d+Bh4DfRMTz1cvOALoi4iUqu2G8r8Z9p1DZleAwKvvQTwQ+Wz0/G5hRXU4Cvkfl4PIJPW7/pog4qsf3+zKVg9hnUNnlYRnwtjqzD4ZHgFep7Mrxb9WvN/zl7MNUdrX4DXAV8H9bPIuGGBu1UeXNRm1UebNRG22lSKneq93alEXE6pTS6Ii4HDgZWFW9ajTwuZTSNyJiKXA0lf3jr0wpvbWf269OKQ1kP3xJVTYq5c1GpbzZ6KYhlw/CVbmCSoRf7+O666gEuwuVt8zu7/aSBp+NSnmzUSlvNjqEuUuwoLJbwVkRMRogIsZHxE7V6+YBp1EJ+boCt2+JiPhkVD7suffyo1Y+rpQJG5XyZqNS3mx0CHOX4A7Wc5eGiDgXeH/1qtXA+1JK/1297pfA8ymlt/W4b5+3dzcJafDYqJQ3G5XyZqObBjdYJUmSJElZGhLHsI4ZMyZNnDix7DGkUi1evPj5lNLYsufoi41KNirlzkalvNVqdEhssE6cOJFFixaVPYZUqoh4ouwZarFRyUal3NmolLdajfqmS5IkSZKkLLnBKkmSJEnKkhuskiRJkqQsDYljWFWy886rnM6e3fQqfve737FixQrWrl07KCNtykaOHMmECRMYPnx42aNoE2OHg8NG1ajzli0DYPZeew3K+my5PhtVO9hh8xpt1A1W9W/JkgGvYsWKFWy99dZMnDiRiBj4TJuolBIrV65kxYoV7LHHHmWPo6Fi2rTK6YIFdW9mhwNno2rGktWrB3V9tlybjaoZ0+67D4AFhxxS+D522JxmGnWXYLXF2rVr2XHHHQ26HxHBjjvu6F/r1BJ2OHA2qhzYcm02qnaxw+Y006gbrGobgy7G50mt5M/XwPkcKgf+HNbmc6N28WetOY0+b26wqmNEBGeccUb3+XXr1jF27FimT58OwFVXXcXYsWOZPHly9/Lwww+XNa60Seqvww1OPPFEDj/88I0uu+SSSxg/fvxGja5ataodY0vqpazfqR/72Me617f33nuz3XbbDXid0lCV4+/UU089tXt9EydOZPLkyQNep8ewqmNstdVWPPjgg7z66qtsueWW3H777YwfP36j25x66ql87WtfK2lCadNXpMNVq1Zx7733Mnr0aB5//PGNjnH52Mc+xgUXXNDusSX1Utbv1K985SvdX3/1q1/lvuqxh1InyvF36vz587u/Pv/889l2220HvE5fYVVHede73sUPfvADAK655hre+973NryOuXPncthhhzF58mTOOeccXn/9dQBGjx7Nxz/+cd70pjdxzDHHcM899zBt2jT23HNPbrrpJqDyF+cPf/jD3euaPn06C/p5oxypX6ecUlmGiP46vP7665kxYwannXYa8+bNa3j9X/rSl5g6dSoHHXQQF198MQBdXV3su+++vP/972fSpEmcfvrp3HHHHbz1rW9lr7324p577gEqf3G+/PLLu9c1adIkurq6mvxOpU3bYPxOvfTSS9l33315xzvewXvf+97u/v71X/+VqVOncvDBB3PSSSexZs2aN9y32ceUejtlp504Zaedyh6jKa38nbp+/XpmzZrFAQccwPTp0zn++OO57rrrAPjMZz7D1KlTmTRpEmeffTYppY3um1Li2muvHZRGfYVV7XfeeYPyzsMbmTy50MfunHbaaXzmM59h+vTpPPDAA5x11ln89Kc/7b5+/vz5LFy4sPv8XXfdxZZbbtl9funSpcyfP5+f/exnDB8+nFmzZnH11VfzZ3/2Z7zyyitMmzaNL3zhC/zpn/4pF110EbfffjsPP/wwZ555JieccMJgfsfS782a1dz9Nry7cE+nnFJZ35o1cPzxb7x+5szK8vzzcPLJG19X8I8v/XV4zTXXcPHFF7Pzzjtz8sknc+GFF3Zf95WvfIW5c+cCsP322/PjH/94o3XfdtttLFu2jHvuuYeUEieccAI/+clP2G233Xj00Uf53ve+x5w5c5g6dSrf/e53WbhwITfddBOf/exnufHGGwvNL+XkvGXLBv1diCePHl3oI3gG+jt10aJFXH/99dx3332sW7eOQw89lDe96U0AvPvd7+YDH/gAABdddBHf+MY3+MhHPtJ93yeeeILHH3+ct7/97QP+fqVZvV6VbMa0Pl7tP2WnnZg1fjxrXn+d4x944A3Xz9xlF2aOG8fzr73GyQ89tNF1Rd+xuJW/U2+44Qa6urr45S9/ybPPPst+++3HWWedBcCHP/xhPv3pTwNwxhlncMsttzBjxozu+/70pz9l5513Zq9B+DgvN1jVUQ466CC6urq45pprOL6Pf4z3t/vSnXfeyeLFi5k6dSoAr776KjtV/yI3YsQI3vnOdwJw4IEHssUWWzB8+HAOPPBAX6FRa2145WHUqHLnKKheh8888wyPPvooRx55JBHB5ptvzoMPPsikSZOA/ndfuu2227jttts4pPqLfvXq1SxbtozddtuNPfbYgwMPPBCAAw44gKOPPpqIsFGpSQP9nbpw4UJOPPHE7o3Ynv/YffDBB7noootYtWoVq1ev5rjjjtvovvPmzePkk09m2LBhg/TdqJOtqe4tN2oI/jy18nfqwoULec973sNmm23GLrvswtve9rbu63784x/zxS9+kTVr1vDCCy9wwAEHbNTwYO4BUXiDNSK2Sim9MiiPqs5W4JXQVjrhhBO44IILWLBgAStXrmzoviklzjzzTD73uc+94brhw4d3v+vZZpttxhZbbNH99bp16wDYfPPNWb9+ffd9BvNt9220g234BdXo7uX1bj9qVP3rx4xp/PF6qNXh/PnzefHFF7uPsXnppZeYN28el112WaH1ppS48MILOeeccza6vKurq7tJsFFtOoq8EtpKA/2dWsvMmTO58cYbOfjgg7nqqqvecPjMvHnzuOKKK5oZuSYb7VwbXv1s5HNYe6t331HDhtW9fsyIEQN67Fb+Tu3L2rVrmTVrFosWLWLXXXflkksu2ej35bp167jhhhtYvHhx099TT/0ewxoRR0TEw8DS6vmDI+LKQXl0qQRnnXUWn/70p7tfaWnE0UcfzXXXXcezzz4LwAsvvMATTzxR+P4TJ05kyZIlrF+/nuXLl3cfNzcQNqqhqFaH11xzDbfeeitdXV10dXWxePHiho65Oe644/jmN7/J6uoukk899VR3r0VMnDiRe++9F4B7772Xxx9/vPB9a7FRbcoG8jv1yCOP5Oabb2bt2rWsXr26+zg8gJdffplx48bxu9/9jquvvnqj+z3yyCO8+OKLb3jX02bZqIa6Vv1OPfLII7n++utZv349zzzzTPcfjjZsnI4ZM4bVq1d3H9e6wR133MG+++7LhAkTBvaNVRV5hfUrwHHATQAppfsj4qhBeXSpBBMmTODcc8/t87rex9tceeWVHHHEEd3n999/fy677DKOPfZY1q9fz/Dhw7niiivYfffdCz32W9/61u7dEidNmsShhx46sG+mwkY15PTVYVdXF08++SRvectbui/bY4892Gabbbj77ruBjY+3AbjxxhuZOHFi9/ljjz2WpUuXdv9DdvTo0cydO7fwboMnnXQS3/72t5k8eTJTp05l7733bvZb7MlGtckayO/UqVOncsIJJ3DwwQez++67M2XKlO53FL300kt585vfzO67786BBx7Iyy+/3H2/a665htNOO20wPwPTRjWktep36kknncSdd97JpEmT2HvvvXnzm9/Mtttuy3bbbccHPvABDjzwQCZOnNh9qNwG8+bNG9Q3RIt6u2MARMTdKaU3R8R9KaVDqpfdn1I6eNCm6MeUKVPSokWL2vVw6m3DG7MMYPe/pUuXst9++w3KOJ2gr+crIhanlKb0vq2Nqmijdjh4bFSN2PBmLAPZ5a+nTanl1atXM3r0aNasWcNRRx3FnDlzBuWPuTaqRjTT6KbUYT0bGl25ciWHHXYYP/vZz9hll10GvN5GGi3yCuvyiDgCSBExAvgo1V0mJGXBRqW82ahUw9lnn83DDz/M2rVrOfPMMwdrz6NG2ahUw/Tp01m1ahWvvfYan/rUpwZlY7VRRTZYPwj8IzAeWAHcBnyolUNJaoiNdrqZM8ueQPXZqFTDd7/73bJHABvteDNL2AgbKnq/4VkZ+t1gTSk9D5zehlkkNcFG5QZr3mxUypuNaua4cWWPoDqKvEvwtyJiux7nt4+Ib7Z0Km2S+jteWhWNPk82Kp5/vrIUYIcDZ6PKgS3XZqNq1POvvcbzr73W8P3ssDmNPm/9brACB6WUVvV4gBeBQkckR8R2EXFdRPwqIpZGxOERsUNE3B4Ry6qn2zc0sYakkSNHsnLlSsPuR0qJlStXMnLkyEbuZqOd7uSTK0s/7HDgbFQ5sOXabFTNOPmhhzj5oYcauo8dNqeZRoscw7pZRGxfjZeI2KHg/aByPMCtKaWTqwexjwI+CdyZUvp8RHwC+ATw8cITa0iaMGECK1as4Lnnnit7lOyNHDmy0c+tslEVYoeDw0ZVNluuz0bVDnbYvEYbLRLjl4GfR8SGT4R9D/D3/d0pIrYBjgJmAqSUXgNei4gTgWnVm30LWIARb/KGDx/OHnvsUfYYmyobVSF2WBob1aCy5UFno2qYHbZPkTdd+nZELAbeBgTw7pTSwwXWvSfwHPB/I+JgYDFwLrBzSunp6rqfjoidmp5eko1KmbNRKW82KuWtyDGsAL8CbgC+D6yOiN0K3Gdz4FDgn6sfwvwKlV0iComIsyNiUUQs8qV2qV82KuXNRqW82aiUqX5fYY2IjwAXA88Ar1P5y1MCDurnriuAFSmlu6vnr6MS8TMRMa76F6dxwLN93TmlNAeYAzBlyhSPZpZqsFHxl39Z9gSqw0alvNmo/nL8+LJHUB1FjmE9F9gnpbSykRWnlH4TEcsjYp+U0iPA0cDD1eVM4PPV0+83OLOkjdlopzv11LInUH02KuXNRjvcqTu513bOimywLgd+2+T6PwJcXX3XtMeAP6eyG/K1EfEXwJNUDmyX1Dwb7XTLl1dOd9213DlUi41KebPRDrd87VoAdm3s45DUJkU2WB8DFkTED4D/2XBhSukf+rtjSmkJMKWPq44uOqCkftlopzvjjMrpggWljqGabFTKm412uDOWLgVgwSGFPn5XbVZkg/XJ6jKiukjKi41KebNRKW82KmWsyMfa/B1ARGyVUnql9SNJaoSNSnmzUSlvNirlrd+PtYmIwyPiYWBp9fzBEXFlyyeTVIiNSnmzUSlvNirlrcjnsM4GjgNWAqSU7geOauFMkhozGxuVcjYbG5VyNhsblbJV5BhWUkrLI6LnRa+3ZhxJzbDRDnf++WVPoH7YqJQ3G+1s5/su+1kr9LE2EXEEkKpv2f1RqrtMSMqCjXa6GTPKnkD12aiUNxvtcDPGjCl7BNVRZJfgDwIfAsYDK4DJ1fOS8mCjne6RRyqLcmWjUt5stMM9smYNj6xZU/YYqqHuK6wRMQyYnVI6vU3zSGqAjQqAc86pnPo5rNmxUSlvNiqAc6p/9PVzWPNU9xXWlNLrwNjq7hGSMmOjUt5sVMqbjUr5K3IMaxfws4i4Cej+bKqU0j+0aihJDenCRqWcdWGjUs66sFEpW0U2WH9dXTYDtm7tOJKaYKNS3mxUypuNShnrd4M1pfR3ABGxVUrplf5uL6m9bFTKm41KebNRKW/9vktwRBweEQ9TfXvviDg4Iq5s+WSSCrFRcdFFlUVZslEpbzaqi3bfnYt2373sMVRDkY+1mQ0cB6wESCndDxzVwpkkNWY2NtrZjjmmsihXs7FRKWezsdGOdswOO3DMDjuUPYZqKLLBSkppea+LXm/BLJKaZKMdbsmSyqJs2aiUNxvtbEtefpklL79c9hiqocibLi2PiCOAVH3L749S3WVCUhZstNOdd17l1M9hzZWNSnmz0Q533qOPAn4Oa66KvML6QeBDwHhgBTC5el5SHmxUypuNSnmzUSljNV9hjYgvpJQ+DrwtpXR6G2eSVICNSnmzUSlvNioNDfVeYT0+IoYDF7ZrGEkNsVEpbzYq5c1GpSGg3jGstwLPA1tFxEtAAGnDaUppmzbMJ6k2G5XyZqNS3mxUGgLqbbBelFL664j4fkrpxLZNJKkoG1XFZz9b9gTqm41KebNRAfDZPfcsewTVUW+D9S7gUOClNs0iqTE2qoojjih7AvXNRqW82agAOGLbbcseQXXU22AdERFnAkdExLt7X5lSuqF1Y0kqwEZV8fOfV07dcM2NjUp5s1EB8PPf/hZwwzVX9TZYPwicDmwHzOh1XQKMWCqXjarik5+snPo5rLmxUSlvNioAPvnYY4Cfw5qrmhusKaWFwMKIWJRS+kYbZ5JUgI1KebNRKW82Kg0N9T6H9e0ppX8HXnQ3CSk/NirlzUalvNmoNDTU2yX4j4B/5427SIC7SUg5sFEpbzYq5c1GpSGg3i7BF1dP/7x940gqykalvNmolDcblYaGersE/1W9O6aU/mHwx5FUlI2q2+zZZU+gPtiolDcb1Qaz//APyx5BddTbJXjr6uk+wFTgpur5GcBPWjmUpEJsVBWTJ5c9gfpmo1LebFQATN566/5vpNLU2yX47wAi4jbg0JTSy9XzlwDfK/oAETEMWAQ8lVKaHhE7APOBiUAXcEpK6cUm55c6lo2q2x13VE6POabcObSRwWjUPqXWsVFtcMcLLwBwzA47lDyJ+rJZgdvsBrzW4/xrVCIs6lxgaY/znwDuTCntBdxZPS+peTba6S67rLIoVwNp1D6l1rPRDnfZE09w2RNPlD2Gaiiywfod4J6IuCQiLgbuBr5VZOURMQH4Y+D/9Lj4xB73/xbwJ4WnldQXG5Xy1lSj9im1jY1KGat3DCsAKaW/j4gfAf+7etGfp5TuK7j+2cDf8PtjBAB2Tik9XV330xGxUwPzSurFRqW8DaDR2din1HI2KuWt3w1WgJTSvcC9jaw4IqYDz6aUFkfEtEYHi4izgbMBdtttt0bvLnUUG5Xy1mijA+2zug4blQqyUSlfRXYJbtZbgRMioguYB7w9IuYCz0TEOIDq6bN93TmlNCelNCWlNGXs2LEtHFPqWDYq5WtAfYKNSi1mo1KbtGyDNaV0YUppQkppInAa8O8ppfdRecvwM6s3OxP4fqtmkFSbjW5Cvv71yqJNhn1KebPRTcvX99mHr++zT9ljqIZCuwQPss8D10bEXwBPAu8pYQZJtdnoUOMv2U5in1LebHQI2mfUqLJHUB1NbbBGxJyU0tlFb59SWgAsqH69Eji6mceVVIyNdpibb66czphR7hwqrJFG7VNqPxvtLDc//zwAM8aMKXkS9aXfDdaIGJZSer3Xxe57JmXCRsWXv1w5dYM1SzYq5c1G9eXlywE3WHNV5BjWxyNiTkQcHREBkFJa3OK5JBVno1LebFTKm41KGSuywboPcAfwISpBfy0ijmztWJIaYKNS3mxUypuNShnrd4M1pfRqSunalNK7gUOAbYD/aPlkkgqxUSlvNirlzUalvBX6WJuI+KOIuJLKByqPBE5p6VSSGmKjUt5sVMqbjUr5KvKmS48DS4Brgb9OKb3S6qEkFWej4jvfKXsC1WGjUt5sVN/Zb7+yR1AdRT7W5uCU0kstn0RSs2y00+26a9kTqD4blfJmox1u15Ejyx5BddTcYI2Iv0kpfRG4rPqGaRtJKX20lYNJqs9G1W3+/MrpqaeWO4c2YqNS3mxUG8x/9lkATt1pp5InUV/qvcK6tHrq23pLebJRVfzzP1dO3WDNjY1KebNRAfDPTz0FuMGaq5obrCmlm6un32rfOJKKslEpbzYq5c1GpaGh3i7BNwOp1vUppRNaMpGkQmxUypuNSnmzUWloqLdL8OXV03cDuwBzq+ffC3S1cCZJxdiolDcblfJmo9IQUG+X4P8AiIhLU0pH9bjq5oj4Scsnk1SXjUp5s1EpbzYqDQ1FPtZmbETsmVJ6DCAi9gDGtnYsSQ2w0U533XVlT6D6bFTKm412uOsOOKDsEVRHkQ3WjwELIuKx6vmJwNktm0hSo2y0040ZU/YEqs9GpbzZaIcbM2JE2SOojn43WFNKt0bEXsC+1Yt+lVL6n9aOJakoGxVXXVU5nTmzzClUg41KebNRXfX00wDMHDeu5EnUlyKvsFKN9v4WzyKpSTba4dxgzZ6NSnmz0c521W9+A7jBmqvNyh5AkiRJkqS+uMEqSZIkScpSvxusEfEXvc4Pi4iLWzeSpEbYqJQ3G5XyZqNS3oq8wnp0RPwwIsZFxCTgP4GtWzyXpOJsVMqbjUp5s1EpY0XeJfj/i4hTgV8Ca4D3ppR+1vLJJBVio+KHPyx7AtVho1LebFQ/POigskdQHUV2Cd4LOBe4HugCzoiIUS2eS1JBNipGjaosypKNSnmzUY0aNoxRw4aVPYZqKLJL8M3Ap1NK5wB/BCwDftHSqSQ1wkY73ZVXVhblykalvNloh7vyqae48qmnyh5DNRT5HNbDUkovAaSUEvDliLiptWNJaoCNdrprr62czppV7hyqxUalvNloh7v22WcBmDV+fMmTqC9FNlj/JCL6unzZIM8iqTk2KuXNRqW82aiUsSIbrFN7fD0SOBq4F/h2SyaS1CgblfJmo1LebFTKWJF3Cf5Iz/MRsS3wnZZNJKkhNirlzUalvNmolLcib7rU2xpgr8EeRNKgsVEpbzYq5c1GpYz0+wprRNwMpOrZzYD9gWtbOZSk4mxULFhQ9gSqw0alvNmoFhxySNkjqI4ix7Be3uPrdcATKaUVLZpHUuNsVMqbjUp5s1EpY0WOYf2PZlYcEbtSOVh9F2A9MCel9I8RsQMwH5hI5cOZT0kpvdjMY0iyUQGXV/+tdcEF5c6hPtmolDcb1eVPPgnABbvtVvIk6kvNY1gj4uWIeKnH6Us9zxdY9zrg/JTSfsBbgA9FxP7AJ4A7U0p7AXdWz0tqkI2q2y23VBZlxUalvNmoNrhl5UpuWbmy7DFUQ71XWA9OKT3W7IpTSk8DT1e/fjkilgLjgROBadWbfQtYAHy82ceROpiNSnmzUSlvNioNAfXeJfh7ABFx50AfJCImAocAdwM7VwPfEPpONe5zdkQsiohFzz333EBHkDZFNirlzUalvNmoNATUe4V1s4i4GNg7Iv6q95UppX8o8gARMRq4HjgvpfRSRBQaLKU0B5gDMGXKlNTPzaVOZKNS3mxUypuNSkNAvVdYTwPWUtmo3bqPpV8RMZxKwFenlG6oXvxMRIyrXj8OeLa50aWOZ6Oq2HLLyqLc2KiUNxsVAFsOG8aWw4aVPYZqqPkKa0rpEeALEfFASulHja44Kn9e+gawtNdfqG4CzgQ+Xz39fqPrlmSj6uFHDf/nVxvYqJQ3G9UGPzrooLJHUB1FPtam2X8JvRU4A/hlRCypXvZJKvFeGxF/ATwJvKfJ9UvCRqXc2aiUNxuV8tbvBmuzUkoLgVo78R/dqseVVIyNbkIuvbRy+qlPlTuHBpWNSnmz0U3HpV1dAHxq4sRS51Df6n0O67h2DiKpMTaqbnfeWVmUFRuV8maj2uDOF1/kzhdfLHsM1VDvFdZvRsT2VD476lZgYUppXVumklSEjUp5s1EpbzYqDQH13nTpXRExksoHH/8pcHlEPEkl6FtTSk+2Z0RJfbFRKW82KuXNRqWhoe4xrCmltVSjBYiIPYB3AV+LiF1SSoe1fkRJtdiolDcblfJmo1L+GnrTpZTS48CVwJURMaI1I0lqlo12qB13LHsCFWSjUt5stDPtOHx42SOojqbfJTil9NpgDiJpcNloB7n++rInUBNsVMqbjXaO6ydNKnsE1VHzXYIlSZIkSSpTvxusEbFTH5ft05pxJDXKRsWFF1YWZclGpbzZqC587DEufOyxssdQDUVeYf1pRJyy4UxEnA/8v9aNJKlBNtrp7rqrsihXNirlzUY73F2//S13/fa3ZY+hGoocwzoNmBMR7wF2BpYCvmOalI9p2KiUs2nYqJSzadiolK1+X2FNKT1N5a2+DwcmAt9OKa1u8VySCrJRKW82KuXNRqW89fsKa0TcDjwNTAImAN+MiJ+klC5o9XCS+mejUt5sVMqbjUp5K3IM6xUppT9LKa1KKT0IHAG4k7eUDxvtdBMmVBblykalvNloh5uwxRZM2GKLssdQDf2+wppSurHX+XXApa0aSFJjbFTMnVv2BKrDRqW82ajm7r9/2SOojiK7BL8MpOrZEcBwYHVKadtWDiapGBuV8majUt5sVMpbkVdYt+55PiL+BN85TcqGjYrzzquczp5d5hSqwUalvNmozlu2DIDZe+1V8iTqS5FjWDdS3W3i7YM/iqTBYKMdaMmSyqIhwUalvNlo51myejVLVvvG0Lkqskvwu3uc3QyYwu93m5BUMhuV8majUt5sVMpbvxuswIweX68DuoATWzKNpGbYqJQ3G5XyZqNSxoocw/rn7RhEUnNsVMqbjUp5s1EpbzU3WCPiq9TZHSKl9NGWTCSpEBtVt733LnsC9cFGpbzZqDbYe9SoskdQHfVeYV3UtikkNcNGVTFnTtkTqG82KuXNRgXAnH32KXsE1VFvg/WYlNIZEXFuSukf2zaRpKJsVMqbjUp5s1FpCKj3sTZviojdgbMiYvuI2KHn0q4BJdVko6o4++zKotzYqJQ3GxUAZz/yCGc/8kjZY6iGeq+w/gtwK7AnsBiIHtel6uWSymOjqviv/yp7AvXNRqW82agA+K81a8oeQXXUfIU1pfRPKaX9gG+mlPZMKe3RYzFgqWQ2KuXNRqW82ag0NNTbJRiAlNJftmMQSc2xUSlvNirlzUalvPW7wSpJkiRJUhnqHcMqSRoKJk8uewJJkoasyaNHlz2C6nCDVZKGutmzy55AkqQha/Zee5U9guooZZfgiHhnRDwSEY9GxCfKmEFSbTYq5c1GpbzZqDR42v4Ka0QMA64A3gGsAH4RETellB5u9ywqyN0NO4qNSnmz0aHH3Q07i41Kg6uMXYIPAx5NKT0GEBHzgBMBI86Vuxt2GhuV8majQ4y7G3YcG5UGURm7BI8Hlvc4v6J62UYi4uyIWBQRi5577rm2DSfJRqXM2aiUNxuVBlEZG6zRx2XpDRekNCelNCWlNGXs2LFtGEtSlY1KebNRKW82Kg2iMjZYVwC79jg/Afh1CXNI6puNSnmzUSlvNioNojI2WH8B7BURe0TECOA04KYS5pDUNxuV8majUt5sVBpEbX/TpZTSuoj4MPBvwDDgmymlh9o9h6S+2aiUNxuV8maj0uCKlN6wS312IuI54Ik2PNQY4Pk2PE5ROc3jLLW1a57dU0pZHuRio1lwlr61cxYbzeu/PeQ1j7P0zUbp2EZzmgXymqdTZ+mz0SGxwdouEbEopTSl7Dk2yGkeZ6ktt3k2Zbk91znN4yx9y2mWTpDb853TPM7St5xm6QQ5Pd85zQJ5zeMsGyvjGFZJkiRJkvrlBqskSZIkKUtusG5sTtkD9JLTPM5SW27zbMpye65zmsdZ+pbTLJ0gt+c7p3mcpW85zdIJcnq+c5oF8prHWXrwGFZJkiRJUpZ8hVWSJEmSlKWO22CNiB0i4vaIWFY93b7G7d4ZEY9ExKMR8Yk+rr8gIlJEjClrloj4UkT8KiIeiIj/FxHbNTFDf99nRMQ/Va9/ICIOLXrfZjQ7T0TsGhE/joilEfFQRJxb1iw9rh8WEfdFxC0DnaWT2Gixdfe43kZttK1stNi6e1xvozbaVjZabN09rm9bozn1OZB5elzfnkZTSh21AF8EPlH9+hPAF/q4zTDgv4E9gRHA/cD+Pa7flcqHQT8BjClrFuBYYPPq11/o6/79PH7d77N6m+OBHwEBvAW4u+h9m3g+BjLPOODQ6tdbA/81kHkGMkuP6/8K+C5wS9k/90NpsdFBa8JGbbQli40OWhM2aqMtWWx00JoY1EZz6nOg8/S4vi2NdtwrrMCJwLeqX38L+JM+bnMY8GhK6bGU0mvAvOr9NvgK8DfAQA8AHtAsKaXbUkrrqrf7T2BCg4/f3/e5YcZvp4r/BLaLiHEF79uopudJKT2dUroXIKX0MrAUGF/GLAARMQH4Y+D/DGCGTmWjBdbda0YbtdF2stEC6+41o43aaDvZaIF195qxHY3m1OeA5oH2NtqJG6w7p5SeBqie7tTHbcYDy3ucX1G9jIg4AXgqpXR/2bP0chaVv4A0osi6a92m6FztmqdbREwEDgHuLnGW2VT+R79+ADN0KhttbN022twss7HRZtloY+u20eZmmY2NNstGG1t3uxrNqc/BmGc2bWp081Y/QBki4g5glz6u+tuiq+jjshQRo6rrOLbsWXo9xt8C64Cri85VdN11blPkvo0ayDyVKyNGA9cD56WUXipjloiYDjybUlocEdMGMMMmy0YLs9EWzGKj/bPRwmy0BbPYaP9stLCcGs2pzwHN0+5GN8kN1pTSMbWui4hnNry0Xn1J+9k+braCyr77G0wAfg38L2AP4P6I2HD5vRFxWErpN22eZcM6zgSmA0enlBqNqO66+7nNiAL3bdRA5iEihlOJ+OqU0g0lznIycEJEHA+MBLaJiLkppfcNcKZNho0WZqOtmcVG+2Gjhdloa2ax0X7YaGE5NZpTnwOdp72NphYeIJvjAnyJjQ/+/mIft9kceIxKsBsOQj6gj9t1MbAD0Qc0C/BO4GFgbJOP3+/3SWXf9J4HW9/TyHPUxnkC+DYwe5B+TpqepddtpuGbRTT63NvoIPwc2qiNtmqx0cH5ObRRG23VYqOD83M42I3m1OdA5+l1m5Y3WnpU7V6AHYE7gWXV0x2ql/8B8MMetzueyjtw/TfwtzXWNdCIBzQL8CiV/cqXVJd/aWKGN6wb+CDwwerXAVxRvf6XwJRGnqN2zQMcSWU3hgd6PB/HlzFLr3W0POJNbbHRwfs5tFEbbcVio4P3c2ijNtqKxUYH7+dwsBvNqc+BPjc91tHyRqP6QJIkSZIkZaUT3yVYkiRJkjQEuMEqSZIkScqSG6ySJEmSpCy5wSpJkiRJypIbrJIkSZKkLLnB2sEiYvVQWKfUqWxUypuNSnmz0U2DG6ySJEmSpCy5wSoAIuKvI+IXEfFARPxd9bIvRMSsHre5JCLOr3V7Sa1jo1LebFTKm40OXW6wiog4FtgLOAyYDLwpIo4C5gGn9rjpKcD36txeUgvYqJQ3G5XyZqND2+ZlD6AsHFtd7queHw3slVL6RkTsFBF/AIwFXkwpPRkRH+3r9sBP2jy31ClsVMqbjUp5s9EhzA1WAQTwuZTS1/u47jrgZGAXKn+F6u/2kgafjUp5s1EpbzY6hLlLsAD+DTgrIkYDRMT4iNipet084DQqIV9X4PaSBp+NSnmzUSlvNjqE+QqrSCndFhH7AXdFBMBq4H3AsymlhyJia+CplNLT/d2+lG9A2sTZqJQ3G5XyZqNDW6SUyp5BkiRJkqQ3cJdgSZIkSVKW3GCVJEmSJGXJDVZJkiRJUpbcYJUkSZIkZckNVkmSJElSltxglSRJkiRlyQ1WSZIkSVKW3GCVJEmSJGXp/wfFSa82q2yL2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x432 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig1 = plt.figure('compare_NN_results', figsize=(16,6))\n",
    "for ix1,model_name in enumerate(model_dict.keys()):\n",
    "    ax1 = fig1.add_subplot(2,4,(4*ix1) + 1,title=f'results for {model_name}')\n",
    "    ax1.plot(me_p[1:],yax, '-r', label='ME emu')\n",
    "    ax1.set_xlabel('level')\n",
    "    ax1.set_ylabel('flux / flux div. difference')\n",
    "    ax1.legend(loc='upper right')\n",
    "  \n",
    "    ax1 = fig1.add_subplot(2,4,(4*ix1) + 2, title=f'results for {model_name}')\n",
    "    ax1.plot(mae_p[1:],yax, '--r', label='MAE emu')\n",
    "    ax1.set_xlabel('level')\n",
    "    ax1.set_ylabel('flux / flux div. difference')\n",
    "    ax1.legend(loc='upper right')\n",
    "    \n",
    "    ax1 = fig1.add_subplot(2,4, (4*ix1) + 3, title=f'results for {model_name}')\n",
    "    plt.plot(me_ctl[1:],yax, '-c', label='ME ga7')\n",
    "    ax1.set_xlabel('level')\n",
    "    ax1.set_ylabel('flux / flux div. difference')\n",
    "    ax1.legend(loc='upper right')\n",
    "    \n",
    "    ax1 = fig1.add_subplot(2,4, (4*ix1) + 4, title=f'results for {model_name}')\n",
    "    ax1.plot(mae_ctl[1:],yax, '--c', label='MAE ga7')\n",
    "    ax1.set_xlabel('level')\n",
    "    ax1.set_ylabel('flux / flux div. difference')\n",
    "    ax1.legend(loc='upper right')\n",
    "    \n"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}