a short notebook to look at the distribution of heavy atom changes as a function of default core mapping in timemachine rbfes.

tm_rbfe_heavy_atom_change_distribution_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "608563cb-172d-4ef7-b43f-f9af1e4426be",
   "metadata": {},
   "source": [
    "this is a small internal consistency check to make sure I am actually parameterizing molecules with espaloma (there was a mild inconsistency previously.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e7a3f1c6-1995-46e5-b545-4de5dc7fc93b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "from importlib import resources\n",
    "from typing import no_type_check\n",
    "import pickle\n",
    "\n",
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "import pytest\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "\n",
    "from timemachine import potentials\n",
    "from timemachine.fe import topology\n",
    "from timemachine.fe.topology import BaseTopology, DualTopology, DualTopologyMinimization\n",
    "from timemachine.fe.utils import get_mol_name, get_romol_conf, read_sdf, set_romol_conf\n",
    "from timemachine.ff import Forcefield, make_mol_omm_sys\n",
    "from timemachine.ff.handlers.openmm_deserializer import deserialize_system\n",
    "from timemachine.ff.handlers.bonded import annotate_mol_sys_torsions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6a64a544-73fb-4ec1-8e53-39b1828a85cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import espaloma as esp\n",
    "esp_model = esp.get_model(\"latest\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "48b839e6-9c52-4dd9-ab59-6b85e58cbeae",
   "metadata": {},
   "outputs": [],
   "source": [
    "ligand_file = '/data1/choderaj/rufad/tm/jak2_esp1/ligands_from_FEP_cations_SUBSET_IC50_uM_largercore.sdf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e0b65ba4-58f3-418a-aa3d-eceb6583cc6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from openff.toolkit.topology import Molecule\n",
    "\n",
    "mols = read_sdf(ligand_file) # read mols from sdf\n",
    "mol_dict = {get_mol_name(mol): mol for mol in mols} # make a mol dict\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6a9dfdfe-8444-4749-b667-81679341a1f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rufad/micromamba/envs/timemachine_off/lib/python3.11/site-packages/dgl/heterograph.py:92: DGLWarning: Recommend creating graphs by `dgl.graph(data)` instead of `dgl.DGLGraph(data)`.\n",
      "  dgl_warning(\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9c4a4d208b894d79a2117d49446729fc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning: Cannot perform Hydrogen sampling with GPU-Omega: GPU-Omega disabled.\n"
     ]
    }
   ],
   "source": [
    "charged_mol, omm_sys, tm_ff, molecule_graph = make_mol_omm_sys(mols[0], charge_spec = 'nn', esp_model=esp_model)\n",
    "off_charged_mol, off_omm_sys, off_tm_ff, off_molecule_graph = make_mol_omm_sys(mols[0], charge_spec='am1bcc', esp_model=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e41873bb-eee3-4dfb-aad7-729a0ab63a0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "annotate_mol_sys_torsions(charged_mol, omm_sys, molecule_graph, tm_ff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6eb76335-3372-4cff-a9dd-0abee7bdfa5b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n"
     ]
    }
   ],
   "source": [
    "bt = BaseTopology(charged_mol, tm_ff)\n",
    "ann_sys = bt.setup_end_state()\n",
    "\n",
    "from timemachine.ff.handlers.openmm_deserializer import deserialize_system\n",
    "manual_sys, masses = deserialize_system(omm_sys, 1.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7214071c-c134-4ba8-aa3a-ad8d1bf26112",
   "metadata": {},
   "outputs": [],
   "source": [
    "conf = get_romol_conf(charged_mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "05881688-282e-48de-aa26-c149627b2716",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'HarmonicBond': Array(28.115944, dtype=float32),\n",
       " 'HarmonicAngle': Array(615.1383, dtype=float32),\n",
       " 'PeriodicTorsion': Array(224.76894, dtype=float32),\n",
       " 'Nonbonded': Array(93.094185, dtype=float32)}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{bp.potential.__class__.__name__: bp.potential(conf, bp.params, None) for bp in manual_sys}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7dff3b7-7b78-4457-a856-b00c9d0a713c",
   "metadata": {},
   "source": [
    "make a manual system from the 2.1.0 parameterized system...and check energies..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "69dc03a4-8b81-40eb-9133-6913dc36b7a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "off_manual_sys, off_masses = deserialize_system(off_omm_sys, 1.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c05e267b-7e48-438e-b0e3-b556941e42ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'HarmonicBond': Array(23.45145, dtype=float32),\n",
       " 'HarmonicAngle': Array(121.89885, dtype=float32),\n",
       " 'PeriodicTorsion': Array(4.4882736, dtype=float32),\n",
       " 'Nonbonded': Array(81.44742, dtype=float32)}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{bp.potential.__class__.__name__: bp.potential(conf, bp.params, None) for bp in off_manual_sys}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "3f6b5951-1c40-4b2e-aadf-e6d7f5bd60ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(28.115944, dtype=float32)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ann_sys.bond(conf, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "de2a5182-feb9-4bac-9b6c-f06f3effa351",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(615.1383, dtype=float32)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ann_sys.angle(conf, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "267e8a7b-737d-4d2e-acb3-ed80b26dadf8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(225.09843, dtype=float32)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ann_sys.torsion(conf, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "23b7c5f1-9fb7-4c23-9dac-60f1fe90b3c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(93.09424, dtype=float32)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ann_sys.nonbonded(conf, None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd5270ce-4b8a-4133-99a5-cc55d84daa35",
   "metadata": {},
   "source": [
    "this has appropriately shown consistencies/inconsistencies where they would be expected"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40b2d332-f5fe-455a-acc4-3fdfdd9641f9",
   "metadata": {},
   "source": [
    "now i want to pull out the actual (hybrid) system and make sure that that esp charges are actually being used..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8b94f5f4-bd46-40bb-be98-7228be25307e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets take this from the most recent calcs...\n",
    "# first edge is this:CHEMBL2062809 CHEMBL3645117\n",
    "pkl_path = '/data1/choderaj/rufad/tm/jak2_esp2/success_rbfe_result_CHEMBL2062809_CHEMBL3645117.pkl'\n",
    "\n",
    "# Load the pickle file\n",
    "with open(pkl_path, 'rb') as file:\n",
    "    data = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "31e70f33-db83-4343-9292-478416208a0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets just pull out solvent result\n",
    "\n",
    "solvent_res = data[4].final_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "91dd017e-86c8-405e-b5b1-2f9816e9dc97",
   "metadata": {},
   "outputs": [],
   "source": [
    "sis0 = solvent_res.initial_states[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "b9e11a70-bd1f-4078-8766-21aa0ca1b2c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6282, 6283, 6284, 6285, 6286, 6287, 6288, 6289, 6290, 6291, 6292,\n",
       "       6293, 6294, 6295, 6296, 6297, 6298, 6299, 6300, 6301, 6302, 6303,\n",
       "       6304, 6305, 6306, 6307, 6308, 6309, 6310, 6311, 6312, 6313, 6314,\n",
       "       6315, 6316, 6317, 6318, 6319, 6320, 6321, 6322, 6323, 6324, 6325,\n",
       "       6326, 6327, 6328, 6329, 6330, 6331, 6332, 6333, 6334, 6335, 6336,\n",
       "       6337, 6338, 6339, 6340, 6341, 6342, 6343, 6344, 6345, 6346, 6347,\n",
       "       6348, 6349, 6350, 6351, 6352, 6353, 6354, 6355, 6356, 6357, 6358,\n",
       "       6359], dtype=int32)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sis0.ligand_idxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "d0f36623-0725-42eb-8471-ee7e7edf744e",
   "metadata": {},
   "outputs": [],
   "source": [
    "potentials = {potential.potential.__class__.__name__: potential for potential in sis0.potentials}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "81b41685-d79a-4cd2-8985-49493cdb3ad3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1347, 4)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "potentials['PeriodicTorsion'].potential.idxs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "75b46762-4b11-4dde-acc6-367ffc773361",
   "metadata": {},
   "outputs": [],
   "source": [
    "# only esp makes so many torsion terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "c1513459-c874-4d4d-af2f-010bb356f5c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# now, i just want to look at charge parameterization..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "7e179dbc-c45c-4793-a843-6b541dd2789f",
   "metadata": {},
   "outputs": [],
   "source": [
    "mola = data[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "dfcb8646-e5d2-41a4-b99c-3724f30cc46e",
   "metadata": {},
   "outputs": [],
   "source": [
    "esp_mola_charges = []\n",
    "for atom in mola.GetAtoms():\n",
    "    esp_mola_charges.append(atom.GetProp(\"PartialCharge\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "9e4d985c-31aa-4ef6-b0b0-d8f55d528062",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rufad/micromamba/envs/timemachine_off/lib/python3.11/site-packages/dgl/heterograph.py:92: DGLWarning: Recommend creating graphs by `dgl.graph(data)` instead of `dgl.DGLGraph(data)`.\n",
      "  dgl_warning(\n"
     ]
    }
   ],
   "source": [
    "my_mola = \n",
    "charged_mola, omm_sys, tm_ff, molecule_graph = make_mol_omm_sys(mola, charge_spec = 'nn', esp_model=esp_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "529f1d1b-5b70-434b-8694-7b48ce40979a",
   "metadata": {},
   "outputs": [],
   "source": [
    "esp_mola_charges_2 = []\n",
    "for atom in charged_mola.GetAtoms():\n",
    "    esp_mola_charges_2.append(atom.GetProp(\"PartialCharge\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "b026a116-8026-46cc-8042-7812a203a9c3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.14006730914115906'"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "esp_mola_charges[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "37655aac-d111-4004-be5f-568fa9abd150",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.14006730914115906'"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "esp_mola_charges_2[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a88d569b-bb64-4d1b-8a69-db6d1bb6a9d0",
   "metadata": {},
   "source": [
    "so the partial charges came from espaloma...that's good to know...now we are done with this part."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96f514d3-25de-4ee3-b61a-4f5410de6587",
   "metadata": {},
   "source": [
    "now, I want to look at the distribution of core atom maps..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "ced8113d-1238-492a-bec4-79a4230acdb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from timemachine.fe.utils import plot_atom_mapping_grid, read_sdf, get_mol_name\n",
    "import typing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "db58f6da-7267-40fb-8062-17b9355e2bf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_edges_txt(edges_txt_file: str, **unused_kwargs) -> typing.Tuple[str]:\n",
    "    with open(edges_txt_file) as file: \n",
    "        lines = file.readlines() \n",
    "        lines = [line.rstrip() for line in lines] \n",
    "    as_tuple = [q.split() for q in lines]\n",
    "    return as_tuple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "670fd8c9-ef78-4772-b7db-f27b780636d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "edges_txt = '/data1/choderaj/rufad/tm/jak2_esp1/edges.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "1fab9805-e688-4b7e-a4d8-375358c99010",
   "metadata": {},
   "outputs": [],
   "source": [
    "edges = read_edges_txt(edges_txt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "eae479d4-40ef-4945-b465-e9790e0cebcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mol_dict = {get_mol_name(mol): mol for mol in mols} # make a mol dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "205f96c3-5597-4128-8fa3-432cf404c78b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mol_edge_dict = {(e[0], e[1]): (mol_dict[e[0]], mol_dict[e[1]]) for e in edges}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "43c13791-e096-44bf-b07f-79bfc40f3f3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from timemachine.fe import atom_mapping, model_utils\n",
    "from timemachine.constants import DEFAULT_ATOM_MAPPING_KWARGS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "a54d5e18-e117-4a27-af13-795225f5c601",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_num_uniques(mol_a, mol_b):\n",
    "    mol_a_natoms, mol_b_natoms = mol_a.GetNumAtoms(), mol_b.GetNumAtoms()\n",
    "    \n",
    "    mol_a_nums = [atom.GetAtomicNum() for atom in mol_a.GetAtoms()]\n",
    "    mol_b_nums = [atom.GetAtomicNum() for atom in mol_b.GetAtoms()]\n",
    "    \n",
    "    core = atom_mapping.get_cores(mol_a, mol_b, **DEFAULT_ATOM_MAPPING_KWARGS)[0]\n",
    "    \n",
    "    unique_old_indices = [i for i in range(mol_a_natoms) if i not in core[:,0]]\n",
    "    unique_new_indices = [i for i in range(mol_b_natoms) if i not in core[:,1]]\n",
    "    \n",
    "    num_heavy_olds = [i for i in unique_old_indices if mol_a_nums[i] > 1]\n",
    "    num_heavy_news = [i for i in unique_new_indices if mol_b_nums[i] > 1]\n",
    "\n",
    "    return len(num_heavy_olds) + len(num_heavy_news)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "01af148a-13d3-4b1b-894f-29f833bd49fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_olds_news_map = {key: get_num_uniques(*val) for key, val in mol_edge_dict.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "58d5ef66-7191-49b9-be29-4d5828b459b8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('CHEMBL2062809', 'CHEMBL3645117'): 11,\n",
       " ('CHEMBL2062809', 'CHEMBL3642402'): 0,\n",
       " ('CHEMBL3645054', 'CHEMBL3645117'): 12,\n",
       " ('CHEMBL3639411', 'CHEMBL3642389'): 10,\n",
       " ('CHEMBL3639411', 'CHEMBL3642281'): 12,\n",
       " ('CHEMBL3642281', 'CHEMBL3642389'): 11,\n",
       " ('CHEMBL3642297', 'CHEMBL3642389'): 7,\n",
       " ('CHEMBL3642389', 'CHEMBL3642402'): 8,\n",
       " ('CHEMBL3642389', 'CHEMBL3645060'): 4,\n",
       " ('CHEMBL3642285', 'CHEMBL3642297'): 3,\n",
       " ('CHEMBL3642285', 'CHEMBL3642291'): 5,\n",
       " ('CHEMBL3642285', 'CHEMBL3642290'): 5,\n",
       " ('CHEMBL3642291', 'CHEMBL3642297'): 4,\n",
       " ('CHEMBL3642290', 'CHEMBL3642297'): 2,\n",
       " ('CHEMBL3642291', 'CHEMBL3644979'): 10,\n",
       " ('CHEMBL3642294', 'CHEMBL3642402'): 6,\n",
       " ('CHEMBL3642294', 'CHEMBL3645060'): 8,\n",
       " ('CHEMBL3645002', 'CHEMBL3645060'): 6,\n",
       " ('CHEMBL3642350', 'CHEMBL3645060'): 8,\n",
       " ('CHEMBL3645060', 'CHEMBL3645065'): 5,\n",
       " ('CHEMBL3644981', 'CHEMBL3645026'): 3,\n",
       " ('CHEMBL3644981', 'CHEMBL3645084'): 8,\n",
       " ('CHEMBL3645002', 'CHEMBL3645026'): 10,\n",
       " ('CHEMBL3645002', 'CHEMBL3645084'): 6,\n",
       " ('CHEMBL3645084', 'CHEMBL3645091'): 4,\n",
       " ('CHEMBL3645024', 'CHEMBL3645084'): 7,\n",
       " ('CHEMBL3644983', 'CHEMBL3645032'): 5,\n",
       " ('CHEMBL3644979', 'CHEMBL3644983'): 5,\n",
       " ('CHEMBL3644979', 'CHEMBL3645032'): 0,\n",
       " ('CHEMBL3645032', 'CHEMBL3645041'): 5,\n",
       " ('CHEMBL3645032', 'CHEMBL3645079'): 8,\n",
       " ('CHEMBL3644986', 'CHEMBL3645054'): 9,\n",
       " ('CHEMBL3644986', 'CHEMBL3645051'): 9,\n",
       " ('CHEMBL3644986', 'CHEMBL3645035'): 14,\n",
       " ('CHEMBL3644987', 'CHEMBL3645054'): 11,\n",
       " ('CHEMBL3645049', 'CHEMBL3645054'): 17,\n",
       " ('CHEMBL3644991', 'CHEMBL3645054'): 12,\n",
       " ('CHEMBL3645035', 'CHEMBL3645054'): 13,\n",
       " ('CHEMBL3645036', 'CHEMBL3645054'): 3,\n",
       " ('CHEMBL3645054', 'CHEMBL3645082'): 5,\n",
       " ('CHEMBL3645054', 'CHEMBL3645055'): 12,\n",
       " ('CHEMBL3645054', 'CHEMBL3645092'): 11,\n",
       " ('CHEMBL3644987', 'CHEMBL3645051'): 13,\n",
       " ('CHEMBL3645002', 'CHEMBL3645024'): 5,\n",
       " ('CHEMBL3642350', 'CHEMBL3645002'): 14,\n",
       " ('CHEMBL3645002', 'CHEMBL3645085'): 7,\n",
       " ('CHEMBL3645030', 'CHEMBL3645079'): 0,\n",
       " ('CHEMBL3645030', 'CHEMBL3645036'): 3,\n",
       " ('CHEMBL3645079', 'CHEMBL3645085'): 6,\n",
       " ('CHEMBL3645079', 'CHEMBL3645091'): 2,\n",
       " ('CHEMBL3645049', 'CHEMBL3645079'): 15,\n",
       " ('CHEMBL3644979', 'CHEMBL3645041'): 5,\n",
       " ('CHEMBL3645036', 'CHEMBL3645045'): 0,\n",
       " ('CHEMBL3645045', 'CHEMBL3645091'): 0,\n",
       " ('CHEMBL3645041', 'CHEMBL3645080'): 2,\n",
       " ('CHEMBL3645041', 'CHEMBL3645046'): 5,\n",
       " ('CHEMBL3645051', 'CHEMBL3645056'): 5,\n",
       " ('CHEMBL3644991', 'CHEMBL3645051'): 14,\n",
       " ('CHEMBL3645051', 'CHEMBL3645055'): 14,\n",
       " ('CHEMBL3645056', 'CHEMBL3645071'): 3,\n",
       " ('CHEMBL3642287', 'CHEMBL3645056'): 1,\n",
       " ('CHEMBL3645082', 'CHEMBL3645092'): 5,\n",
       " ('CHEMBL3645065', 'CHEMBL3645071'): 2,\n",
       " ('CHEMBL3642287', 'CHEMBL3645071'): 2,\n",
       " ('CHEMBL3645071', 'CHEMBL3645074'): 0,\n",
       " ('CHEMBL3645065', 'CHEMBL3645074'): 2,\n",
       " ('CHEMBL3645085', 'CHEMBL3645087'): 3,\n",
       " ('CHEMBL3645029', 'CHEMBL3645085'): 0,\n",
       " ('CHEMBL3645046', 'CHEMBL3645080'): 3,\n",
       " ('CHEMBL3645029', 'CHEMBL3645087'): 0}"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_olds_news_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "e703641d-7a27-49cd-83ac-8c6e0e9d7426",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "919d4d8c-b453-4f62-b97c-bfe012600c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(list(n_olds_news_map.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "91ffb1d6-2423-4e4a-a847-88d8ba6f357a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11,  0, 12, 10, 12, 11,  7,  8,  4,  3,  5,  5,  4,  2, 10,  6,  8,\n",
       "        6,  8,  5,  3,  8, 10,  6,  4,  7,  5,  5,  0,  5,  8,  9,  9, 14,\n",
       "       11, 17, 12, 13,  3,  5, 12, 11, 13,  5, 14,  7,  0,  3,  6,  2, 15,\n",
       "        5,  0,  0,  2,  5,  5, 14, 14,  3,  1,  5,  2,  2,  0,  2,  3,  0,\n",
       "        3,  0])"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "04424780-339c-4830-96bf-9baf94ae3150",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'n_old + n_new heavy atom mapping distribution')"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAHFCAYAAAAHcXhbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9bklEQVR4nO3de5xNdf///+dmzmMMhjEzmBnJKeScUzKTnElUjiGHLqEifC6qb6ETcTWpFHVVg09ICrmoGDkWQqT4dCGNQw4paWaYZszh/fuj3+xscx57Zu81Hvfbbd9u9nu/13u91lrv2fO09lqzbcYYIwAAAIsq4+oCAAAArgdhBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBg4WLlwom82m48eP59s3KipKUVFRxV7TjSjrOOzdu9fVpbhUcnKypk+fri1btri6lFLHZrNp+vTpri7D7vjx47LZbFq4cKG9bfr06bLZbIUap6hzJqd1RUZGqmfPnoUaJz9Lly7V3Llzc3zN3Y6JlXi4ugAAyE1ycrJmzJghSQRnJ9u5c6eqV6/u6jLyNGrUKHXt2rVQyxR1zhRlXUWxdOlSHTx4UBMmTMj2mhWOibsizKDEbdmyRdHR0YqPj1dkZKSrywFuSK1bt3Z1CfmqXr16sf9yT05Olp+fX4msKz9WOCbuio+Z3FjWac9Dhw5p4MCBCgwMVNWqVTVixAglJCQUerw1a9aoTZs28vPzU0BAgDp16qSdO3fmu5wxRrNnz1ZERIR8fHzUrFkzffbZZ0XZpOuWddr3888/V7NmzeTr66t69erpvffeK/RYNptNjzzyiP73f/9X9evXl5+fnxo3bqy1a9dm63v06FENGjRIwcHB8vb2Vv369fXGG2/YXzfGqGrVqho3bpy9LSMjQxUrVlSZMmX0yy+/2NtjYmLk4eGhP/74I98ak5KSNGbMGFWuXFlBQUHq27evzpw5k63f8uXL1aZNG/n7+6tcuXLq0qWL9u/f79Bn7969GjBggCIjI+Xr66vIyEgNHDhQJ06csPc5cOCAbDab3n333Wzr+Oyzz2Sz2bRmzRpt375dNptNy5Yty9Zv8eLFstls2rNnT67b9euvv2rs2LG65ZZbVK5cOQUHB+vOO+/U9u3b7X2OHz+uKlWqSJJmzJghm80mm82mBx980N7nyy+/VMeOHRUQECA/Pz+1bdtW69atc1hX1kd2mzZt0kMPPaSgoCCVL19eQ4cO1eXLl3Xu3Dn169dPFSpUUGhoqCZPnqy0tLRca8+SNRfXrl2rpk2bytfXV/Xr17fPn4ULF6p+/fry9/fXbbfdlu0jw4Icj6vrj4uL0/Dhw1WpUiX5+/urV69e+umnnxz6RkVFqWHDhtq+fbtat24tX19fVatWTU8//bQyMjIc+l77kUbWejZv3pzvnEtNTdWkSZMUEhIiPz8/3XHHHfrmm28UGRnpcHxyc+bMGfXr108BAQEKDAxU//79de7cuWz9cvroZ9OmTYqKilJQUJB8fX0VHh6ue++9V8nJyfnOmazx9u3bp/vuu08VK1ZUrVq1cl1XllWrVunWW2+Vj4+PbrrpJr322msOr+f28fyWLVtks9nsH3lFRUVp3bp1OnHihL22q9eZ08dMBw8eVO/evVWxYkX5+PioSZMmWrRoUY7rWbZsmZ566imFhYWpfPnyuuuuu3T48OEct6nUMXBb06ZNM5JM3bp1zTPPPGPi4uJMTEyM8fb2NsOHDy/UWEuWLDGSTOfOnc3q1avN8uXLTfPmzY2Xl5fZvn27vV9sbKyRZOLj47PVMXLkSPPZZ5+Zt99+21SrVs2EhISYDh06FHq7Nm/enG0dBRUREWGqV69ubrnlFrN48WKzfv16c//99xtJZuvWrYUaS5KJjIw0t912m/nwww/Np59+aqKiooyHh4c5duyYvd+hQ4dMYGCgadSokVm8eLHZsGGDmTRpkilTpoyZPn26vd+AAQNMnTp17M937dplJBlfX1+zZMkSe3u3bt3MbbfdlmdtWcfhpptuMo8++qhZv369eeedd0zFihVNdHS0Q98XXnjB2Gw2M2LECLN27VqzcuVK06ZNG+Pv728OHTpk77dixQrzzDPPmFWrVpmtW7eaDz74wHTo0MFUqVLF/Prrr/Z+TZs2Ne3atctWU79+/UxwcLBJS0vLs1/Lli1Ny5Yt89y+//73v2bMmDHmgw8+MFu2bDFr1641I0eONGXKlDGbN282xhiTkpJiPv/8c/vc27lzp9m5c6f58ccfjTHGbNmyxXh6eprmzZub5cuXm9WrV5vOnTsbm81mPvjgg2z7smbNmmbSpElmw4YN5qWXXjJly5Y1AwcONM2aNTPPP/+8iYuLM1OmTDGSzMsvv5xn/cb8PRcbNmxoli1bZj799FPTqlUr4+npaZ555hnTrl07s3LlSrNq1SpTp04dU7VqVZOcnFzo45FVf40aNcyIESPsP4PBwcGmRo0a5uLFi/a+HTp0MEFBQSYsLMy89tprZv369eaxxx4zksy4ceMc6pdkpk2blm09BZlzAwcONGXKlDFTp041GzZsMHPnzjU1atQwgYGBZtiwYXnut+TkZFO/fn0TGBhoXn/9dXuN4eHhRpKJjY21981678kSHx9vfHx8TKdOnczq1avNli1bzJIlS8yQIUPMxYsX850zWeNFRESYKVOmmLi4OLN69eoc15V1jKtVq2bCw8PNe++9Zz799FMzePBgI8nMmTMn27679j0t670ua04fOnTItGvXzoSEhNhr27lzZ67H5L///a8JCAgwtWrVMosXLzbr1q0zAwcONJLMSy+9lG09kZGRZvDgwWbdunVm2bJlJjw83NSuXdukp6fneUxKA8KMG8v64Zo9e7ZD+9ixY42Pj4/JzMws0DgZGRkmLCzMNGrUyGRkZNjbk5KSTHBwsGnbtq297dofyosXLxofHx/Tp08fhzG/+uorI6lAYSY9Pd2kpaXZHxs3bjSSzI8//ujQfnVtuYmIiDA+Pj7mxIkT9rY///zTVKpUyYwePTrf5a8myVStWtUkJiba286dO2fKlCljZs6caW/r0qWLqV69uklISHBY/pFHHjE+Pj7m999/N8YY88477xhJ5uTJk8YYY55//nlTr149c/fdd9vD55UrV4y/v7958skn86wt6ziMHTvWoX327NlGkjl79qwxxpiTJ08aDw8P8+ijjzr0S0pKMiEhIaZfv365riM9Pd1cunTJ+Pv7m1dffdXe/tprrxlJ5vDhw/a233//3Xh7e5tJkyZlq3H//v32tt27dxtJZtGiRXluX061pKWlmY4dOzrMtV9//TXbG3yW1q1bm+DgYJOUlOQwTsOGDU316tXtPx9ZdV67j+655x4jycTExDi0N2nSxDRr1izfmiMiIoyvr6/5+eef7W3ffvutkWRCQ0PN5cuX7e2rV682ksyaNWvy3Ac5HY+s+nP7GXz++eftbR06dDCSzCeffOLQ96GHHjJlypRx+LnJLczkN+cOHTpkJJkpU6Y49Fu2bJmRlG+YmT9/fq415hdmPvroIyPJfPvtt7mOn9ecyRrvmWeeyfW1q0VERBibzZZtfZ06dTLly5e3H+OChhljjOnRo4eJiIjIsfZr6x4wYIDx9va2v6dk6datm/Hz8zN//PGHw3q6d+/u0O/DDz80khwCU2nFx0wWcPfddzs8v/XWW5WSkqLz588XaPnDhw/rzJkzGjJkiMqU+fuQlytXTvfee6927dql5OTkHJfduXOnUlJSNHjwYIf2tm3bKiIiokDr79ixozw9Pe2Pu+66S5J08803O7SPGDGiQOM1adJE4eHh9uc+Pj6qU6dOttPzBREdHa2AgAD786pVqyo4ONg+VkpKir744gv16dNHfn5+Sk9Ptz+6d++ulJQU7dq1S5Ls27Vx40ZJUlxcnDp16qS77rpLcXFxkv7an5cvX7b3zU9Ox16Svb7169crPT1dQ4cOdajNx8dHHTp0cLij49KlS5oyZYpuvvlmeXh4yMPDQ+XKldPly5f1ww8/2PsNHjxY3t7eDneVLFu2TKmpqRo+fLi9beDAgQoODnb4uO31119XlSpV1L9//3y3bcGCBWrWrJl8fHzk4eEhT09PffHFFw615Oby5cv6+uuvdd9996lcuXL29rJly2rIkCH6+eefs51ev/aulPr160uSevToka29oHOpSZMmqlatWrYxo6Ki5Ofnl6396nELejyy5PYzuHnzZof2gICAbPNm0KBByszM1LZt2/Ldpvzm3NatWyVJ/fr1c+h33333ycMj/8swN2/enGuN+WnSpIm8vLz0j3/8Q4sWLcr2MVtB3XvvvQXu26BBAzVu3NihbdCgQUpMTNS+ffuKtP6C2rRpkzp27KgaNWo4tD/44INKTk7OdplAfseuNCPMWEBQUJDDc29vb0nSn3/+WaDlL1y4IEkKDQ3N9lpYWJgyMzN18eLFPJcNCQnJ9lpObTl56623tGfPHvtjwYIFkv66hufq9oLeknjt/pD+2icF3R+FGevChQtKT0/X66+/7hC8PD091b17d0nSb7/9JkmKiIhQrVq1tHHjRvsbTVaYyfrlunHjRvn6+qpt27ZFqu/aY591LU7Lli2z1bd8+XJ7bdJfb8Dz5s3TqFGjtH79eu3evVt79uxRlSpVHPZdpUqVdPfdd2vx4sX26ywWLlyo2267TQ0aNHCoZfTo0Vq6dKn++OMP/frrr/rwww81atQoe525iYmJ0ZgxY9SqVSt9/PHH2rVrl/bs2aOuXbsW6DhevHhRxphc57T099y9eruu5uXllWt7SkpKvjUUdkxJDuMW9Hhkye1n8NrtrFq1aq7LXts3J/nNuawxrl2Ph4dHjj9P17pw4UKeNeYl6+crODhY48aNU61atVSrVi29+uqr+S57tZzmTW7yeu8ryP68HhcuXCjUHL/e3xVWxt1MN4CsCX727Nlsr505c0ZlypRRxYoV81w2p4vzzp07V6C7kerWrevw/NKlS5KkRo0auf3dTBUrVrT/b//qi3uvVrNmTfu/O3bsqE8++URbt25VZmamoqKiFBAQoLCwMMXFxWnjxo1q3759vr/sC6py5cqSpI8++ijPM2UJCQlau3atpk2bpqlTp9rbU1NT9fvvv2frP3z4cK1YsUJxcXEKDw/Xnj17NH/+/Gz9xowZo1mzZum9995TSkqK0tPT9fDDD+db9/vvv6+oqKhsYyYlJeW7rCT7hdW5zWnp733jjgp7PKTcfwZvvvlmh7arLza/dtmChI38ZI3xyy+/OJyVSk9PL3BY2r17d6415qd9+/Zq3769MjIytHfvXr3++uuaMGGCqlatqgEDBhRojML87Zrc9rv0977w8fGR9Nfxu9rV/5koiqCgIMvO8ZLGmZkbQN26dVWtWjUtXbpUxhh7++XLl/Xxxx/b73DKSevWreXj46MlS5Y4tO/YseOGOHXp5+en6Oho7d+/X7feeqtatGiR7XH1L4i77rpLv/zyi+bOnavWrVvbP8Lq2LGjVq1apT179hT4I6aC6NKlizw8PHTs2LEca2vRooWkv968jTHZQtQ777yT7S4XSercubOqVaum2NhYxcbGysfHRwMHDszWLzQ0VPfff7/efPNNLViwQL169XL4CDA3NpstWy3fffddttPmuf3P0t/fX61atdLKlSsdXsvMzNT777+v6tWrq06dOvnW4SqFPR6Scv0ZvPZvqSQlJWnNmjUObUuXLlWZMmV0xx13XHftWWMsX77cof2jjz5Senp6vstHR0fnWmNhlC1bVq1atbJ/zJn1kY+zz0YcOnRIBw4ccGhbunSpAgIC1KxZM0my/6fsu+++c+h37TZm1VfQ2jp27KhNmzZlu5ts8eLF8vPz41buq3Bm5gZQpkwZzZ49W4MHD1bPnj01evRopaamas6cOfrjjz80a9asXJetWLGiJk+erOeff16jRo3S/fffr1OnTmn69OkF/pjJ6l599VXdfvvtat++vcaMGaPIyEglJSXpxx9/1H/+8x9t2rTJ3vfOO++UzWbThg0b7H+4S/or5AwbNsz+b2eJjIzUs88+q6eeeko//fSTunbtqooVK+qXX37R7t275e/vrxkzZqh8+fK64447NGfOHFWuXFmRkZHaunWr3n33XVWoUCHbuGXLltXQoUMVExOj8uXLq2/fvgoMDMyxhvHjx6tVq1aSpNjY2ALV3bNnTz333HOaNm2aOnTooMOHD+vZZ59VzZo1HX4hBgQEKCIiQp988ok6duyoSpUq2eufOXOmOnXqpOjoaE2ePFleXl568803dfDgQS1btqzQfzm2JBX2eEh/3cp99c/gU089pWrVqmns2LEO/YKCgjRmzBidPHlSderU0aeffqp///vfGjNmTIGCZn4aNGiggQMH6uWXX1bZsmV155136tChQ3r55ZcVGBjocF1eToYOHapXXnlFQ4cO1QsvvKDatWvr008/1fr16/Nd94IFC7Rp0yb16NFD4eHhSklJsf9Zhqyfq7zmTFGEhYXp7rvv1vTp0xUaGqr3339fcXFxeumll+z/CWzZsqXq1q2ryZMnKz09XRUrVtSqVav05ZdfZhuvUaNGWrlypebPn6/mzZurTJky9v90XGvatGlau3atoqOj9cwzz6hSpUpasmSJ1q1bp9mzZ+f6M3lDcu31x8hL1tX1V9+maUzuV87nZ/Xq1aZVq1bGx8fH+Pv7m44dO5qvvvoq37EzMzPNzJkzTY0aNYyXl5e59dZbzX/+8x/ToUMHl9ya3aNHj2ztRalFOdyumrWOa+/IiI+PNyNGjDDVqlUznp6epkqVKqZt27YOd5Jkadq0qZHksG9Pnz5tJJmgoKAC3YWWdRz27Nnj0J7T3RHG/HVso6OjTfny5Y23t7eJiIgw9913n9m4caO9z88//2zuvfdeU7FiRRMQEGC6du1qDh48mOP2GmPMkSNHjCQjycTFxeVZb2RkpKlfv36+25UlNTXVTJ482VSrVs34+PiYZs2amdWrV5thw4Zlu9Nj48aNpmnTpsbb2zvb3TLbt283d955p/H39ze+vr6mdevW5j//+Y/D8rnty9x+voYNG2b8/f3z3Ybc5mJO8yo+Pj7b7bwFPR5Z9W/YsMEMGTLEVKhQwfj6+pru3bubo0ePOqynQ4cOpkGDBmbLli2mRYsWxtvb24SGhponn3zSfkv91XXmdDdTQeZcSkqKmThxogkODjY+Pj6mdevWZufOnSYwMNA8/vjj+e67rG0vV66cCQgIMPfee6/ZsWNHvncz7dy50/Tp08dEREQYb29vExQUZDp06JDtLrHc5kxuxzyndRnz9zH+6KOPTIMGDYyXl5eJjIzMdgecMX/9vHTu3NmUL1/eVKlSxTz66KNm3bp12fbd77//bu677z5ToUIFY7PZHNZ57TExxpjvv//e9OrVywQGBhovLy/TuHFjh31kzN/HaMWKFQ7tWfPu2v6lkc2Yqz53AIBC+u6779S4cWO98cYb2c4S4PotXLhQw4cP1549e3L9H3yWqKgo/fbbbzp48GAJVfe3HTt2qF27dlqyZEmB7kwCnImPmQAUybFjx3TixAk9+eSTCg0NLdBffkXpEBcXp507d6p58+by9fXVgQMHNGvWLNWuXVt9+/Z1dXm4ARFmLC4zM1OZmZl59inI334oTfK7CLFMmTL5fq6P/D333HP2r4JYsWJFrheRo/QpX768NmzYoLlz5yopKUmVK1dWt27dNHPmTPudPUBJ4mMmi3vwwQezfU/HtW60Q5zfhZ/Dhg1z+INwAABrI8xY3PHjx/P9Wwb5fc5e2lz7hX7Xup47GwAA7ocwAwAALI0LBwAAgKWV+itDMzMzdebMGQUEBLj1H9ECAAB/M8YoKSlJYWFh+d60UerDzJkzZ7J94ygAALCGU6dOqXr16nn2KfVhJuu7cU6dOqXy5cu7uBoAAFAQiYmJqlGjhv33eF5KfZjJ+mipfPnyhBkAACymIJeIcAEwAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNJeGmW3btqlXr14KCwuTzWbT6tWr7a+lpaVpypQpatSokfz9/RUWFqahQ4fqzJkzrisYAAC4HZeGmcuXL6tx48aaN29etteSk5O1b98+Pf3009q3b59WrlypI0eO6O6773ZBpQAAwF3ZjDHG1UVIf32R1KpVq3TPPffk2mfPnj267bbbdOLECYWHhxdo3MTERAUGBiohIYEvmgQAwCIK8/vbUtfMJCQkyGazqUKFCq4uBQAAuAkPVxdQUCkpKZo6daoGDRqUZ0JLTU1Vamqq/XliYmJJlAcAAFzEEmEmLS1NAwYMUGZmpt588808+86cOVMzZswoocqAnEVOXef0MY/P6uH0MQGgNHD7j5nS0tLUr18/xcfHKy4uLt/PzZ544gklJCTYH6dOnSqhSgEAgCu49ZmZrCBz9OhRbd68WUFBQfku4+3tLW9v7xKoDgAAuAOXhplLly7pxx9/tD+Pj4/Xt99+q0qVKiksLEz33Xef9u3bp7Vr1yojI0Pnzp2TJFWqVEleXl6uKhsAALgRl4aZvXv3Kjo62v584sSJkqRhw4Zp+vTpWrNmjSSpSZMmDstt3rxZUVFRJVUmAABwYy4NM1FRUcrrz9y4yZ/AAQAAbsztLwAGAADIC2EGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYmkvDzLZt29SrVy+FhYXJZrNp9erVDq8bYzR9+nSFhYXJ19dXUVFROnTokGuKBQAAbsmlYeby5ctq3Lix5s2bl+Prs2fPVkxMjObNm6c9e/YoJCREnTp1UlJSUglXCgAA3JWHK1ferVs3devWLcfXjDGaO3eunnrqKfXt21eStGjRIlWtWlVLly7V6NGjS7JUAADgptz2mpn4+HidO3dOnTt3trd5e3urQ4cO2rFjR67LpaamKjEx0eEBAABKL7cNM+fOnZMkVa1a1aG9atWq9tdyMnPmTAUGBtofNWrUKNY6AQCAa7ltmMlis9kcnhtjsrVd7YknnlBCQoL9cerUqeIuEQAAuJBLr5nJS0hIiKS/ztCEhoba28+fP5/tbM3VvL295e3tXez1AQAA9+C2Z2Zq1qypkJAQxcXF2duuXLmirVu3qm3bti6sDAAAuBOXnpm5dOmSfvzxR/vz+Ph4ffvtt6pUqZLCw8M1YcIEvfjii6pdu7Zq166tF198UX5+fho0aJALqwYAAO7EpWFm7969io6Otj+fOHGiJGnYsGFauHCh/vnPf+rPP//U2LFjdfHiRbVq1UobNmxQQECAq0oGAABuxmaMMa4uojglJiYqMDBQCQkJKl++vKvLwQ0icuo6p495fFYPp48JAO6qML+/3faaGQAAgIIgzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEtz6zCTnp6u//f//p9q1qwpX19f3XTTTXr22WeVmZnp6tIAAICb8HB1AXl56aWXtGDBAi1atEgNGjTQ3r17NXz4cAUGBmr8+PGuLg8AALgBtw4zO3fuVO/evdWjRw9JUmRkpJYtW6a9e/e6uDIAAOAu3Ppjpttvv11ffPGFjhw5Ikk6cOCAvvzyS3Xv3j3XZVJTU5WYmOjwAAAApZdbn5mZMmWKEhISVK9ePZUtW1YZGRl64YUXNHDgwFyXmTlzpmbMmFFiNUZOXef0MY/P6uH0MQEAKK3c+szM8uXL9f7772vp0qXat2+fFi1apH/9619atGhRrss88cQTSkhIsD9OnTpVghUDAICS5tZnZv7nf/5HU6dO1YABAyRJjRo10okTJzRz5kwNGzYsx2W8vb3l7e1dkmUCAAAXcuszM8nJySpTxrHEsmXLcms2AACwc+szM7169dILL7yg8PBwNWjQQPv371dMTIxGjBjh6tIAAICbcOsw8/rrr+vpp5/W2LFjdf78eYWFhWn06NF65plnXF0aAABwE24dZgICAjR37lzNnTvX1aUAAAA35dbXzAAAAOSHMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACytSGEmPj7e2XUAAAAUSZHCzM0336zo6Gi9//77SklJcXZNAAAABVakMHPgwAE1bdpUkyZNUkhIiEaPHq3du3c7uzYAAIB8FSnMNGzYUDExMTp9+rRiY2N17tw53X777WrQoIFiYmL066+/OrtOAACAHF3XBcAeHh7q06ePPvzwQ7300ks6duyYJk+erOrVq2vo0KE6e/ass+oEAADI0XWFmb1792rs2LEKDQ1VTEyMJk+erGPHjmnTpk06ffq0evfu7aw6AQAAcuRRlIViYmIUGxurw4cPq3v37lq8eLG6d++uMmX+ykY1a9bUW2+9pXr16jm1WAAAgGsVKczMnz9fI0aM0PDhwxUSEpJjn/DwcL377rvXVRwAAEB+ihRmjh49mm8fLy8vDRs2rCjDAwAAFFiRrpmJjY3VihUrsrWvWLFCixYtuu6iAAAACqpIYWbWrFmqXLlytvbg4GC9+OKL110UAABAQRUpzJw4cUI1a9bM1h4REaGTJ09ed1EAAAAFVaRrZoKDg/Xdd98pMjLSof3AgQMKCgpyRl24QUROXVcs4x6f1aNYxgUAuJ8inZkZMGCAHnvsMW3evFkZGRnKyMjQpk2bNH78eA0YMMDZNQIAAOSqSGdmnn/+eZ04cUIdO3aUh8dfQ2RmZmro0KFcMwMAAEpUkcKMl5eXli9frueee04HDhyQr6+vGjVqpIiICGfXBwAAkKcihZksderUUZ06dZxVCwAAQKEVKcxkZGRo4cKF+uKLL3T+/HllZmY6vL5p0yanFAcAAJCfIoWZ8ePHa+HCherRo4caNmwom83m7LoAAAAKpEhh5oMPPtCHH36o7t27O7seAACAQinSrdleXl66+eabnV0LAABAoRUpzEyaNEmvvvqqjDHOrgcAAKBQivQx05dffqnNmzfrs88+U4MGDeTp6enw+sqVK51SHAAAQH6KFGYqVKigPn36OLsWAACAQitSmImNjXV2HQAAAEVSpGtmJCk9PV0bN27UW2+9paSkJEnSmTNndOnSJacVBwAAkJ8inZk5ceKEunbtqpMnTyo1NVWdOnVSQECAZs+erZSUFC1YsMDZdQIAAOSoSGdmxo8frxYtWujixYvy9fW1t/fp00dffPGF04oDAADIT5HvZvrqq6/k5eXl0B4REaHTp087pTAAAICCKNKZmczMTGVkZGRr//nnnxUQEHDdRQEAABRUkcJMp06dNHfuXPtzm82mS5cuadq0aXzFAQAAKFFF+pjplVdeUXR0tG655RalpKRo0KBBOnr0qCpXrqxly5Y5u0YAAIBcFSnMhIWF6dtvv9WyZcu0b98+ZWZmauTIkRo8eLDDBcEAAADFrUhhRpJ8fX01YsQIjRgxwpn1AAAAFEqRwszixYvzfH3o0KFFKgYAAKCwihRmxo8f7/A8LS1NycnJ8vLykp+fH2EGAACUmCLdzXTx4kWHx6VLl3T48GHdfvvtXAAMAABKVJG/m+latWvX1qxZs7Kdtblep0+f1gMPPKCgoCD5+fmpSZMm+uabb5y6DgAAYF1FvgA4J2XLltWZM2ecNt7FixfVrl07RUdH67PPPlNwcLCOHTumChUqOG0dAADA2ooUZtasWePw3Bijs2fPat68eWrXrp1TCpOkl156STVq1FBsbKy9LTIy0mnjAwAA6ytSmLnnnnscnttsNlWpUkV33nmnXn75ZWfUJemv0NSlSxfdf//92rp1q6pVq6axY8fqoYcecto6AACAtRUpzGRmZjq7jhz99NNPmj9/viZOnKgnn3xSu3fv1mOPPSZvb+9c75hKTU1Vamqq/XliYmKJ1AoAAFzDqdfMOFtmZqZatGihF198UZLUtGlTHTp0SPPnz881zMycOVMzZswoyTJLncip64pl3OOzehTLuACAG1uRwszEiRML3DcmJqYoq5AkhYaG6pZbbnFoq1+/vj7++ONcl3niiScc6ktMTFSNGjWKXAMAAHBvRQoz+/fv1759+5Senq66detKko4cOaKyZcuqWbNm9n42m+26imvXrp0OHz7s0HbkyBFFRETkuoy3t7e8vb2va70AAMA6ihRmevXqpYCAAC1atEgVK1aU9Ndt1MOHD1f79u01adIkpxT3+OOPq23btnrxxRfVr18/7d69W2+//bbefvttp4wPAACsr0h/NO/ll1/WzJkz7UFGkipWrKjnn3/eqXcztWzZUqtWrdKyZcvUsGFDPffcc5o7d64GDx7stHUAAABrK9KZmcTERP3yyy9q0KCBQ/v58+eVlJTklMKy9OzZUz179nTqmAAAoPQo0pmZPn36aPjw4froo4/0888/6+eff9ZHH32kkSNHqm/fvs6uEQAAIFdFOjOzYMECTZ48WQ888IDS0tL+GsjDQyNHjtScOXOcWiAAAEBeihRm/Pz89Oabb2rOnDk6duyYjDG6+eab5e/v7+z6AAAA8nRd35p99uxZnT17VnXq1JG/v7+MMc6qCwAAoECKFGYuXLigjh07qk6dOurevbvOnj0rSRo1apTTbssGAAAoiCKFmccff1yenp46efKk/Pz87O39+/fX559/7rTiAAAA8lOka2Y2bNig9evXq3r16g7ttWvX1okTJ5xSGAAAQEEU6czM5cuXHc7IZPntt9/4KgEAAFCiihRm7rjjDi1evNj+3GazKTMzU3PmzFF0dLTTigMAAMhPkT5mmjNnjqKiorR3715duXJF//znP3Xo0CH9/vvv+uqrr5xdIwAAQK6KdGbmlltu0XfffafbbrtNnTp10uXLl9W3b1/t379ftWrVcnaNAAAAuSr0mZm0tDR17txZb731lmbMmFEcNQEAABRYoc/MeHp66uDBg7LZbMVRDwAAQKEU6WOmoUOH6t1333V2LQAAAIVWpAuAr1y5onfeeUdxcXFq0aJFtu9kiomJcUpxAAAA+SlUmPnpp58UGRmpgwcPqlmzZpKkI0eOOPTh4ycAAFCSChVmateurbNnz2rz5s2S/vr6gtdee01Vq1YtluIAAADyU6hrZq79VuzPPvtMly9fdmpBAAAAhVGkC4CzXBtuAAAASlqhwozNZst2TQzXyAAAAFcq1DUzxhg9+OCD9i+TTElJ0cMPP5ztbqaVK1c6r0IAAIA8FCrMDBs2zOH5Aw884NRiAAAACqtQYSY2Nra46gAAACiSIv3RPAAoTpFT1xXLuMdn9SiWcQG41nXdzQQAAOBqhBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBphBkAAGBplgozM2fOlM1m04QJE1xdCgAAcBOWCTN79uzR22+/rVtvvdXVpQAAADdiiTBz6dIlDR48WP/+979VsWJFV5cDAADciCXCzLhx49SjRw/ddddd+fZNTU1VYmKiwwMAAJReHq4uID8ffPCB9u3bpz179hSo/8yZMzVjxoxirgpwH5FT1zl9zOOzejh9TJSskpwXpXEOFsc2Sa7frtLKrc/MnDp1SuPHj9f7778vHx+fAi3zxBNPKCEhwf44depUMVcJAABcya3PzHzzzTc6f/68mjdvbm/LyMjQtm3bNG/ePKWmpqps2bIOy3h7e8vb27ukSwUAAC7i1mGmY8eO+v777x3ahg8frnr16mnKlCnZggwAALjxuHWYCQgIUMOGDR3a/P39FRQUlK0dAADcmNz6mhkAAID8uPWZmZxs2bLF1SUAAAA3wpkZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaR6uLgAoKZFT1zl9zOOzejh9THdVHPtPYh9erxtp/wG54cwMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNLcOMzNnzlTLli0VEBCg4OBg3XPPPTp8+LCrywIAAG7ErcPM1q1bNW7cOO3atUtxcXFKT09X586ddfnyZVeXBgAA3ISHqwvIy+eff+7wPDY2VsHBwfrmm290xx13uKgqAADgTtw6zFwrISFBklSpUqVc+6Smpio1NdX+PDExsdjrAgAArmOZMGOM0cSJE3X77berYcOGufabOXOmZsyYUYKVAbCyyKnrimXc47N6FMu4NzKOFXLj1tfMXO2RRx7Rd999p2XLluXZ74knnlBCQoL9cerUqRKqEAAAuIIlzsw8+uijWrNmjbZt26bq1avn2dfb21ve3t4lVBkAAHA1tw4zxhg9+uijWrVqlbZs2aKaNWu6uiQAAOBm3DrMjBs3TkuXLtUnn3yigIAAnTt3TpIUGBgoX19fF1cHAADcgVtfMzN//nwlJCQoKipKoaGh9sfy5ctdXRoAAHATbn1mxhjj6hIAAICbc+szMwAAAPkhzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEsjzAAAAEvzcHUBKJjIqeuKZdzjs3oUy7gAgIIpyff34liXO/we4cwMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNMIMAACwNEuEmTfffFM1a9aUj4+Pmjdvru3bt7u6JAAA4CbcPswsX75cEyZM0FNPPaX9+/erffv26tatm06ePOnq0gAAgBtw+zATExOjkSNHatSoUapfv77mzp2rGjVqaP78+a4uDQAAuAG3DjNXrlzRN998o86dOzu0d+7cWTt27HBRVQAAwJ14uLqAvPz222/KyMhQ1apVHdqrVq2qc+fO5bhMamqqUlNT7c8TEhIkSYmJicVSY2ZqstPHzKnW4lgP6yqe9ZTWdVn9WJXWdbl6XpTkuqx+rErruorr92vWuMaY/DsbN3b69GkjyezYscOh/fnnnzd169bNcZlp06YZSTx48ODBgwePUvA4depUvnnBrc/MVK5cWWXLls12Fub8+fPZztZkeeKJJzRx4kT788zMTP3+++8KCgqSzWYr1npzk5iYqBo1aujUqVMqX768S2pwF+wLR+yPv7Ev/sa++Bv74m832r4wxigpKUlhYWH59nXrMOPl5aXmzZsrLi5Offr0sbfHxcWpd+/eOS7j7e0tb29vh7YKFSoUZ5kFVr58+RtiAhYE+8IR++Nv7Iu/sS/+xr742420LwIDAwvUz63DjCRNnDhRQ4YMUYsWLdSmTRu9/fbbOnnypB5++GFXlwYAANyA24eZ/v3768KFC3r22Wd19uxZNWzYUJ9++qkiIiJcXRoAAHADbh9mJGns2LEaO3asq8soMm9vb02bNi3bx183IvaFI/bH39gXf2Nf/I198Tf2Re5sxhTknicAAAD35NZ/NA8AACA/hBkAAGBphBkAAGBphBkAAGBphBknefPNN1WzZk35+PioefPm2r59e579t27dqubNm8vHx0c33XSTFixYUEKVFp+ZM2eqZcuWCggIUHBwsO655x4dPnw4z2W2bNkim82W7fHf//63hKouPtOnT8+2XSEhIXkuUxrnhSRFRkbmeJzHjRuXY//SNC+2bdumXr16KSwsTDabTatXr3Z43Rij6dOnKywsTL6+voqKitKhQ4fyHffjjz/WLbfcIm9vb91yyy1atWpVMW2B8+S1L9LS0jRlyhQ1atRI/v7+CgsL09ChQ3XmzJk8x1y4cGGOcyUlJaWYt+b65DcvHnzwwWzb1Lp163zHteK8cAbCjBMsX75cEyZM0FNPPaX9+/erffv26tatm06ePJlj//j4eHXv3l3t27fX/v379eSTT+qxxx7Txx9/XMKVO9fWrVs1btw47dq1S3FxcUpPT1fnzp11+fLlfJc9fPiwzp49a3/Url27BCoufg0aNHDYru+//z7XvqV1XkjSnj17HPZDXFycJOn+++/Pc7nSMC8uX76sxo0ba968eTm+Pnv2bMXExGjevHnas2ePQkJC1KlTJyUlJeU65s6dO9W/f38NGTJEBw4c0JAhQ9SvXz99/fXXxbUZTpHXvkhOTta+ffv09NNPa9++fVq5cqWOHDmiu+++O99xy5cv7zBPzp49Kx8fn+LYBKfJb15IUteuXR226dNPP81zTKvOC6e47m+DhLntttvMww8/7NBWr149M3Xq1Bz7//Of/zT16tVzaBs9erRp3bp1sdXoCufPnzeSzNatW3Pts3nzZiPJXLx4seQKKyHTpk0zjRs3LnD/G2VeGGPM+PHjTa1atUxmZmaOr5fWeSHJrFq1yv48MzPThISEmFmzZtnbUlJSTGBgoFmwYEGu4/Tr18907drVoa1Lly5mwIABTq+5uFy7L3Kye/duI8mcOHEi1z6xsbEmMDDQucWVsJz2xbBhw0zv3r0LNU5pmBdFxZmZ63TlyhV988036ty5s0N7586dtWPHjhyX2blzZ7b+Xbp00d69e5WWllZstZa0hIQESVKlSpXy7du0aVOFhoaqY8eO2rx5c3GXVmKOHj2qsLAw1axZUwMGDNBPP/2Ua98bZV5cuXJF77//vkaMGJHvl7+W1nmRJT4+XufOnXM47t7e3urQoUOu7x9S7nMlr2WsKCEhQTabLd/v17t06ZIiIiJUvXp19ezZU/v37y+ZAovZli1bFBwcrDp16uihhx7S+fPn8+x/o8yLnBBmrtNvv/2mjIyMbN/iXbVq1Wzf9p3l3LlzOfZPT0/Xb7/9Vmy1liRjjCZOnKjbb79dDRs2zLVfaGio3n77bX388cdauXKl6tatq44dO2rbtm0lWG3xaNWqlRYvXqz169fr3//+t86dO6e2bdvqwoULOfa/EeaFJK1evVp//PGHHnzwwVz7lOZ5cbWs94jCvH9kLVfYZawmJSVFU6dO1aBBg/L8UsV69epp4cKFWrNmjZYtWyYfHx+1a9dOR48eLcFqna9bt25asmSJNm3apJdffll79uzRnXfeqdTU1FyXuRHmRW4s8XUGVnDt/zCNMXn+rzOn/jm1W9Ujjzyi7777Tl9++WWe/erWrau6devan7dp00anTp3Sv/71L91xxx3FXWax6tatm/3fjRo1Ups2bVSrVi0tWrRIEydOzHGZ0j4vJOndd99Vt27dFBYWlmuf0jwvclLY94+iLmMVaWlpGjBggDIzM/Xmm2/m2bd169YOF8a2a9dOzZo10+uvv67XXnutuEstNv3797f/u2HDhmrRooUiIiK0bt069e3bN9flSvO8yAtnZq5T5cqVVbZs2WzJ9/z589kScpaQkJAc+3t4eCgoKKjYai0pjz76qNasWaPNmzerevXqhV6+devWlv9fVU78/f3VqFGjXLettM8LSTpx4oQ2btyoUaNGFXrZ0jgvsu5uK8z7R9ZyhV3GKtLS0tSvXz/Fx8crLi4uz7MyOSlTpoxatmxZ6uZKaGioIiIi8tyu0jwv8kOYuU5eXl5q3ry5/e6MLHFxcWrbtm2Oy7Rp0yZb/w0bNqhFixby9PQstlqLmzFGjzzyiFauXKlNmzapZs2aRRpn//79Cg0NdXJ1rpeamqoffvgh120rrfPiarGxsQoODlaPHj0KvWxpnBc1a9ZUSEiIw3G/cuWKtm7dmuv7h5T7XMlrGSvICjJHjx7Vxo0bixTijTH69ttvS91cuXDhgk6dOpXndpXWeVEgLrv0uBT54IMPjKenp3n33XfN//3f/5kJEyYYf39/c/z4cWOMMVOnTjVDhgyx9//pp5+Mn5+fefzxx83//d//mXfffdd4enqajz76yFWb4BRjxowxgYGBZsuWLebs2bP2R3Jysr3PtfvilVdeMatWrTJHjhwxBw8eNFOnTjWSzMcff+yKTXCqSZMmmS1btpiffvrJ7Nq1y/Ts2dMEBATccPMiS0ZGhgkPDzdTpkzJ9lppnhdJSUlm//79Zv/+/UaSiYmJMfv377ffoTNr1iwTGBhoVq5cab7//nszcOBAExoaahITE+1jDBkyxOHuyK+++sqULVvWzJo1y/zwww9m1qxZxsPDw+zatavEt68w8toXaWlp5u677zbVq1c33377rcN7SGpqqn2Ma/fF9OnTzeeff26OHTtm9u/fb4YPH248PDzM119/7YpNLLC89kVSUpKZNGmS2bFjh4mPjzebN282bdq0MdWqVSuV88IZCDNO8sYbb5iIiAjj5eVlmjVr5nA78rBhw0yHDh0c+m/ZssU0bdrUeHl5mcjISDN//vwSrtj5JOX4iI2Ntfe5dl+89NJLplatWsbHx8dUrFjR3H777WbdunUlX3wx6N+/vwkNDTWenp4mLCzM9O3b1xw6dMj++o0yL7KsX7/eSDKHDx/O9lppnhdZt5lf+xg2bJgx5q/bs6dNm2ZCQkKMt7e3ueOOO8z333/vMEaHDh3s/bOsWLHC1K1b13h6epp69epZIujltS/i4+NzfQ/ZvHmzfYxr98WECRNMeHi48fLyMlWqVDGdO3c2O3bsKPmNK6S89kVycrLp3LmzqVKlivH09DTh4eFm2LBh5uTJkw5jlJZ54Qw2Y/7/KwwBAAAsiGtmAACApRFmAACApRFmAACApRFmAACApRFmAACApRFmAACApRFmAACApRFmAFhSVFSUJkyY4OoyALgBwgyAEterVy/dddddOb62c+dO2Ww27du3r4SrAmBVhBkAJW7kyJHatGmTTpw4ke219957T02aNFGzZs1cUBkAKyLMAChxPXv2VHBwsBYuXOjQnpycrOXLl+uee+7RwIEDVb16dfn5+alRo0ZatmxZnmPabDatXr3aoa1ChQoO6zh9+rT69++vihUrKigoSL1799bx48eds1EAXIYwA6DEeXh4aOjQoVq4cKGu/nq4FStW6MqVKxo1apSaN2+utWvX6uDBg/rHP/6hIUOG6Ouvvy7yOpOTkxUdHa1y5cpp27Zt+vLLL1WuXDl17dpVV65cccZmAXARwgwAlxgxYoSOHz+uLVu22Nvee+899e3bV9WqVdPkyZPVpEkT3XTTTXr00UfVpUsXrVixosjr++CDD1SmTBm98847atSokerXr6/Y2FidPHnSoQYA1uPh6gIA3Jjq1auntm3b6r333lN0dLSOHTum7du3a8OGDcrIyNCsWbO0fPlynT59WqmpqUpNTZW/v3+R1/fNN9/oxx9/VEBAgEN7SkqKjh07dr2bA8CFCDMAXGbkyJF65JFH9MYbbyg2NlYRERHq2LGj5syZo1deeUVz585Vo0aN5O/vrwkTJuT5cZDNZnP4yEqS0tLS7P/OzMxU8+bNtWTJkmzLVqlSxXkbBaDEEWYAuEy/fv00fvx4LV26VIsWLdJDDz0km82m7du3q3fv3nrggQck/RVEjh49qvr16+c6VpUqVXT27Fn786NHjyo5Odn+vFmzZlq+fLmCg4NVvnz54tsoACWOa2YAuEy5cuXUv39/Pfnkkzpz5owefPBBSdLNN9+suLg47dixQz/88INGjx6tc+fO5TnWnXfeqXnz5mnfvn3au3evHn74YXl6etpfHzx4sCpXrqzevXtr+/btio+P19atWzV+/Hj9/PPPxbmZAIoZYQaAS40cOVIXL17UXXfdpfDwcEnS008/rWbNmqlLly6KiopSSEiI7rnnnjzHefnll1WjRg3dcccdGjRokCZPniw/Pz/7635+ftq2bZvCw8PVt29f1a9fXyNGjNCff/7JmRrA4mzm2g+ZAQAALIQzMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNIIMwAAwNL+P+htWSeveCN8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.hist(data, bins=range(min(data), max(data) + 1), align='left', rwidth=0.8)\n",
    "\n",
    "# Add labels and title\n",
    "plt.xlabel('Value')\n",
    "plt.ylabel('Frequency')\n",
    "plt.title('n_old + n_new heavy atom mapping distribution')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "934c397c-61c8-45ab-9d5b-feac8a3446e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'n_old + n_new heavy atom mapping distribution')"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHFCAYAAAAOmtghAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/8UlEQVR4nO3deXxN1/7/8feJSA6RxJgJGUqRolTUWDXHVFW0lH7Neg2t1nS/pR4XdXsvHeTqhPYaUt+iqobblpaouWgN0UFdVVOUpGooqVSIrN8ffjkcmSPJSbbX8/E4j4ezztprf/Y+K/HOHs6xGWOMAAAALMLN1QUAAADkJ8INAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINAACwFMINshUdHS2bzabjx49n27dVq1Zq1apVgdd0N0p7H/bs2ePqUlwqKSlJU6dO1ebNm11diuXYbDZNnTrV1WU4HD9+XDabTdHR0Y62qVOnymaz5WqcvM6ZjNYVGhqqRx55JFfjZGfJkiWaNWtWhq8VtfekuHB3dQEAkBtJSUl66aWXJIkgnc927typKlWquLqMLA0dOlQdO3bM1TJ5nTN5WVdeLFmyRD/88INGjx6d7rXi8J4URYQbFAmbN29W69atdezYMYWGhrq6HOCu1KRJE1eXkK0qVaoU+H/2SUlJKl26dKGsKzvF4T0pijgtVcykHSY9cOCA+vTpI19fX/n7+2vw4MG6ePFirsf75JNP1LRpU5UuXVre3t5q3769du7cme1yxhi9+uqrCgkJkd1uV4MGDfT555/nZZPuWNph4i+++EINGjRQqVKlVKtWLS1YsCDXY9lsNj377LP6v//7P4WHh6t06dKqV6+ePvvss3R9Dx8+rL59+8rPz0+enp4KDw/XO++843jdGCN/f38988wzjrbr16+rXLlycnNz06+//upoj4qKkru7u37//fdsa0xMTNSIESNUsWJFVahQQT169NDp06fT9Vu2bJmaNm0qLy8vlSlTRh06dFBsbKxTnz179ujJJ59UaGioSpUqpdDQUPXp00cnTpxw9Pn2229ls9k0f/78dOv4/PPPZbPZ9Mknn2jbtm2y2WxaunRpun6LFi2SzWbT7t27M92u3377TSNHjtR9992nMmXKyM/PT23atNG2bdscfY4fP65KlSpJkl566SXZbDbZbDYNHDjQ0Wf79u1q27atvL29Vbp0aTVr1kxr1qxxWlfaKb6NGzfq6aefVoUKFeTj46P+/fvr8uXLSkhIUK9evVS2bFkFBgZq/PjxunbtWqa1p0mbi5999pkeeOABlSpVSuHh4Y75Ex0drfDwcHl5ealRo0bpTjHm5P24tf6YmBgNGjRI5cuXl5eXl7p27aqjR4869W3VqpXq1Kmjbdu2qUmTJipVqpQqV66sv/3tb7p+/bpT39tPgaStZ9OmTdnOueTkZI0bN04BAQEqXbq0Hn74Ye3du1ehoaFO709mTp8+rV69esnb21u+vr7q3bu3EhIS0vXL6FTRxo0b1apVK1WoUEGlSpVScHCwevbsqaSkpGznTNp4+/bt0+OPP65y5cqpWrVqma4rzapVq3T//ffLbrfrnnvu0Ztvvun0eman8zdv3iybzeY4RdaqVSutWbNGJ06ccNR26zozOi31ww8/qFu3bipXrpzsdrvq16+v999/P8P1LF26VJMmTVJQUJB8fHzUrl07HTp0KMNtshSDYmXKlClGkqlZs6aZPHmyiYmJMVFRUcbT09MMGjQoV2MtXrzYSDKRkZFm9erVZtmyZSYiIsJ4eHiYbdu2OfotXLjQSDLHjh1LV8eQIUPM559/bt577z1TuXJlExAQYFq2bJnr7dq0aVO6deRUSEiIqVKlirnvvvvMokWLzLp168wTTzxhJJktW7bkaixJJjQ01DRq1Mh89NFHZu3ataZVq1bG3d3dHDlyxNHvwIEDxtfX19StW9csWrTIrF+/3owbN864ubmZqVOnOvo9+eSTpkaNGo7nu3btMpJMqVKlzOLFix3tnTp1Mo0aNcqytrT34Z577jGjRo0y69atM/PmzTPlypUzrVu3dur7j3/8w9hsNjN48GDz2WefmZUrV5qmTZsaLy8vc+DAAUe/5cuXm8mTJ5tVq1aZLVu2mA8//NC0bNnSVKpUyfz222+Ofg888IBp3rx5upp69epl/Pz8zLVr17Ls9+CDD5oHH3wwy+3773//a0aMGGE+/PBDs3nzZvPZZ5+ZIUOGGDc3N7Np0yZjjDFXrlwxX3zxhWPu7dy50+zcudP8/PPPxhhjNm/ebEqWLGkiIiLMsmXLzOrVq01kZKSx2Wzmww8/TLcvw8LCzLhx48z69evNK6+8YkqUKGH69OljGjRoYF5++WUTExNjXnjhBSPJzJw5M8v6jbk5F+vUqWOWLl1q1q5daxo3bmxKlixpJk+ebJo3b25WrlxpVq1aZWrUqGH8/f1NUlJSrt+PtPqrVq1qBg8e7PgZ9PPzM1WrVjUXLlxw9G3ZsqWpUKGCCQoKMm+++aZZt26dee6554wk88wzzzjVL8lMmTIl3XpyMuf69Olj3NzczIQJE8z69evNrFmzTNWqVY2vr68ZMGBAlvstKSnJhIeHG19fX/PWW285agwODjaSzMKFCx190373pDl27Jix2+2mffv2ZvXq1Wbz5s1m8eLFpl+/fubChQvZzpm08UJCQswLL7xgYmJizOrVqzNcV9p7XLlyZRMcHGwWLFhg1q5da5566ikjybz22mvp9t3tv9PSftelzekDBw6Y5s2bm4CAAEdtO3fuzPQ9+e9//2u8vb1NtWrVzKJFi8yaNWtMnz59jCTzyiuvpFtPaGioeeqpp8yaNWvM0qVLTXBwsLn33ntNSkpKlu9JcUe4KWbSftheffVVp/aRI0cau91uUlNTczTO9evXTVBQkKlbt665fv26oz0xMdH4+fmZZs2aOdpu/yG9cOGCsdvtpnv37k5jfvXVV0ZSjsJNSkqKuXbtmuOxYcMGI8n8/PPPTu231paZkJAQY7fbzYkTJxxtf/75pylfvrwZNmxYtsvfSpLx9/c3ly5dcrQlJCQYNzc3M336dEdbhw4dTJUqVczFixedln/22WeN3W4358+fN8YYM2/ePCPJxMXFGWOMefnll02tWrXMo48+6gijV69eNV5eXubFF1/Msra092HkyJFO7a+++qqRZOLj440xxsTFxRl3d3czatQop36JiYkmICDA9OrVK9N1pKSkmD/++MN4eXmZN954w9H+5ptvGknm0KFDjrbz588bT09PM27cuHQ1xsbGOtq++eYbI8m8//77WW5fRrVcu3bNtG3b1mmu/fbbb+l+4adp0qSJ8fPzM4mJiU7j1KlTx1SpUsXx85FW5+376LHHHjOSTFRUlFN7/fr1TYMGDbKtOSQkxJQqVcr88ssvjrb9+/cbSSYwMNBcvnzZ0b569WojyXzyySdZ7oOM3o+0+jP7GXz55ZcdbS1btjSSzH/+8x+nvk8//bRxc3Nz+rnJLNxkN+cOHDhgJJkXXnjBqd/SpUuNpGzDzZw5czKtMbtw8/HHHxtJZv/+/ZmOn9WcSRtv8uTJmb52q5CQEGOz2dKtr3379sbHx8fxHuc03BhjTJcuXUxISEiGtd9e95NPPmk8PT0dv1PSdOrUyZQuXdr8/vvvTuvp3LmzU7+PPvrISHIKUFbEaali6tFHH3V6fv/99+vKlSs6c+ZMjpY/dOiQTp8+rX79+snN7eY0KFOmjHr27Kldu3YpKSkpw2V37typK1eu6KmnnnJqb9asmUJCQnK0/rZt26pkyZKOR7t27SRJ1atXd2ofPHhwjsarX7++goODHc/tdrtq1KiR7nB+TrRu3Vre3t6O5/7+/vLz83OMdeXKFX355Zfq3r27SpcurZSUFMejc+fOunLlinbt2iVJju3asGGDJCkmJkbt27dXu3btFBMTI+nG/rx8+bKjb3Yyeu8lOepbt26dUlJS1L9/f6fa7Ha7WrZs6XTHyB9//KEXXnhB1atXl7u7u9zd3VWmTBldvnxZBw8edPR76qmn5Onp6XTXytKlS5WcnKxBgwY52vr06SM/Pz+n03NvvfWWKlWqpN69e2e7bXPnzlWDBg1kt9vl7u6ukiVL6ssvv3SqJTOXL1/W119/rccff1xlypRxtJcoUUL9+vXTL7/8ku5w/O13vYSHh0uSunTpkq49p3Opfv36qly5croxW7VqpdKlS6drv3XcnL4faTL7Gdy0aZNTu7e3d7p507dvX6Wmpmrr1q3ZblN2c27Lli2SpF69ejn1e/zxx+Xunv2lnZs2bcq0xuzUr19fHh4e+stf/qL3338/3Wm5nOrZs2eO+9auXVv16tVzauvbt68uXbqkffv25Wn9ObVx40a1bdtWVatWdWofOHCgkpKS0l1WkN17Z1WEm2KqQoUKTs89PT0lSX/++WeOlj937pwkKTAwMN1rQUFBSk1N1YULF7JcNiAgIN1rGbVl5N1339Xu3bsdj7lz50q6cQ3Qre05vQXy9v0h3dgnOd0fuRnr3LlzSklJ0VtvveUUxEqWLKnOnTtLks6ePStJCgkJUbVq1bRhwwbHL560cJP2n+2GDRtUqlQpNWvWLE/13f7ep13L8+CDD6arb9myZY7apBu/kN9++20NHTpU69at0zfffKPdu3erUqVKTvuufPnyevTRR7Vo0SLHdRrR0dFq1KiRateu7VTLsGHDtGTJEv3+++/67bff9NFHH2no0KGOOjMTFRWlESNGqHHjxlqxYoV27dql3bt3q2PHjjl6Hy9cuCBjTKZzWro5d2/drlt5eHhk2n7lypVsa8jtmJKcxs3p+5Ems5/B27fT398/02Vv75uR7OZc2hi3r8fd3T3Dn6fbnTt3Lssas5L28+Xn56dnnnlG1apVU7Vq1fTGG29ku+ytMpo3mcnqd19O9uedOHfuXK7m+J3+X1FccbfUXSptwsfHx6d77fTp03Jzc1O5cuWyXDaji/0SEhJydLdTzZo1nZ7/8ccfkqS6desW+bulypUr5zgacOvFwrcKCwtz/Ltt27b6z3/+oy1btig1NVWtWrWSt7e3goKCFBMTow0bNqhFixbZ/uefUxUrVpQkffzxx1keSbt48aI+++wzTZkyRRMmTHC0Jycn6/z58+n6Dxo0SMuXL1dMTIyCg4O1e/duzZkzJ12/ESNGaMaMGVqwYIGuXLmilJQUDR8+PNu6P/jgA7Vq1SrdmImJidkuK8lxoXZmc1q6uW+Koty+H1LmP4PVq1d3arv14vXbl81J+MhO2hi//vqr01GrlJSUHIenb775JtMas9OiRQu1aNFC169f1549e/TWW29p9OjR8vf315NPPpmjMXLz2TmZ7Xfp5r6w2+2Sbrx/t7r1j4u8qFChQrGd44WJIzd3qZo1a6py5cpasmSJjDGO9suXL2vFihWOO6gy0qRJE9ntdi1evNipfceOHZY/1ClJpUuXVuvWrRUbG6v7779fDRs2TPe49T+Mdu3a6ddff9WsWbPUpEkTxymvtm3batWqVdq9e3eOT0nlRIcOHeTu7q4jR45kWFvDhg0l3fhlboxJF6rmzZuX7i4aSYqMjFTlypW1cOFCLVy4UHa7XX369EnXLzAwUE888YRmz56tuXPnqmvXrk6nDDNjs9nS1fLdd9+lO8ye2V+eXl5eaty4sVauXOn0Wmpqqj744ANVqVJFNWrUyLYOV8nt+yEp05/B2z/LJTExUZ988olT25IlS+Tm5qaHH374jmtPG2PZsmVO7R9//LFSUlKyXb5169aZ1pgbJUqUUOPGjR2nRdNOEeX30YoDBw7o22+/dWpbsmSJvL291aBBA0ly/JH23XffOfW7fRvT6stpbW3bttXGjRvT3a22aNEilS5dmlvH/z+O3Nyl3Nzc9Oqrr+qpp57SI488omHDhik5OVmvvfaafv/9d82YMSPTZcuVK6fx48fr5Zdf1tChQ/XEE0/o5MmTmjp1ao5PSxV3b7zxhh566CG1aNFCI0aMUGhoqBITE/Xzzz/r008/1caNGx1927RpI5vNpvXr1zs+SEy6EXoGDBjg+Hd+CQ0N1bRp0zRp0iQdPXpUHTt2VLly5fTrr7/qm2++kZeXl1566SX5+Pjo4Ycf1muvvaaKFSsqNDRUW7Zs0fz581W2bNl045YoUUL9+/dXVFSUfHx81KNHD/n6+mZYw/PPP6/GjRtLkhYuXJijuh955BH9/e9/15QpU9SyZUsdOnRI06ZNU1hYmNN/kN7e3goJCdF//vMftW3bVuXLl3fUP336dLVv316tW7fW+PHj5eHhodmzZ+uHH37Q0qVLc/3JtoUpt++HdOPW8Vt/BidNmqTKlStr5MiRTv0qVKigESNGKC4uTjVq1NDatWv173//WyNGjMhR8MxO7dq11adPH82cOVMlSpRQmzZtdODAAc2cOVO+vr5O1/VlpH///vrXv/6l/v376x//+IfuvfderV27VuvWrct23XPnztXGjRvVpUsXBQcH68qVK46PgUj7ucpqzuRFUFCQHn30UU2dOlWBgYH64IMPFBMTo1deecXxR+GDDz6omjVravz48UpJSVG5cuW0atUqbd++Pd14devW1cqVKzVnzhxFRETIzc3N8UfI7aZMmaLPPvtMrVu31uTJk1W+fHktXrxYa9as0auvvprpz+Rdx7XXMyO30q7ev/W2UGMyvzI/O6tXrzaNGzc2drvdeHl5mbZt25qvvvoq27FTU1PN9OnTTdWqVY2Hh4e5//77zaeffmpatmzpklvBu3Tpkq49L7Uog9tj09Zx+x0fx44dM4MHDzaVK1c2JUuWNJUqVTLNmjVzulMlzQMPPGAkOe3bU6dOGUmmQoUKObrLLe192L17t1N7RndfGHPjvW3durXx8fExnp6eJiQkxDz++ONmw4YNjj6//PKL6dmzpylXrpzx9vY2HTt2ND/88EOG22uMMT/99JORZCSZmJiYLOsNDQ014eHh2W5XmuTkZDN+/HhTuXJlY7fbTYMGDczq1avNgAED0t1JsmHDBvPAAw8YT0/PdHfjbNu2zbRp08Z4eXmZUqVKmSZNmphPP/3UafnM9mVmP18DBgwwXl5e2W5DZnMxo3l17NixdLcP5/T9SKt//fr1pl+/fqZs2bKmVKlSpnPnzubw4cNO62nZsqWpXbu22bx5s2nYsKHx9PQ0gYGB5sUXX3Tcwn9rnRndLZWTOXflyhUzduxY4+fnZ+x2u2nSpInZuXOn8fX1NWPGjMl236Vte5kyZYy3t7fp2bOn2bFjR7Z3S+3cudN0797dhISEGE9PT1OhQgXTsmXLdHehZTZnMnvPM1qXMTff448//tjUrl3beHh4mNDQ0HR32Blz4+clMjLS+Pj4mEqVKplRo0aZNWvWpNt358+fN48//rgpW7assdlsTuu8/T0xxpjvv//edO3a1fj6+hoPDw9Tr149p31kzM33aPny5U7tafPu9v5WYzPmlnMSAJAPvvvuO9WrV0/vvPNOuqMIuHPR0dEaNGiQdu/enelf+GlatWqls2fP6ocffiik6m7asWOHmjdvrsWLF+fozicgv3BaCkC+OXLkiE6cOKEXX3xRgYGBOfpkWlhDTEyMdu7cqYiICJUqVUrffvutZsyYoXvvvVc9evRwdXm4yxBuLCg1NVWpqalZ9snJZ09YSXYXNbq5uWV7XQCy9/e//93x1RXLly/P9KJ0WI+Pj4/Wr1+vWbNmKTExURUrVlSnTp00ffp0x51DQGHhtJQFDRw4MN33jNzubnvbs7uQdMCAAU4fUAcAKL4INxZ0/PjxbD9LIbvz9FZz+xcU3u5O7pwAABQthBsAAGApXGQAAAAs5e66qlQ3LrY9ffq0vL29i/QHegEAgJuMMUpMTFRQUFC2N4DcdeHm9OnT6b5NFQAAFA8nT55UlSpVsuxz14WbtO/1OXnypHx8fFxcDQAAyIlLly6patWqjv/Hs3LXhZu0U1E+Pj6EGwAAipmcXFLCBcUAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSCDcAAMBSXBputm7dqq5duyooKEg2m02rV6/OdpktW7YoIiJCdrtd99xzj+bOnVvwhQIAgGLDpeHm8uXLqlevnt5+++0c9T927Jg6d+6sFi1aKDY2Vi+++KKee+45rVixooArBQAAxYVLvzizU6dO6tSpU477z507V8HBwZo1a5YkKTw8XHv27NHrr7+unj17FlCVAACgOClW19zs3LlTkZGRTm0dOnTQnj17dO3aNRdVBQAAihKXHrnJrYSEBPn7+zu1+fv7KyUlRWfPnlVgYGC6ZZKTk5WcnOx4funSpQKvEwAAuE6xCjeSZLPZnJ4bYzJsTzN9+nS99NJLBV4XAADFTeiENQUy7vEZXQpk3JwqVqelAgIClJCQ4NR25swZubu7q0KFChkuM3HiRF28eNHxOHnyZGGUCgAAXKRYHblp2rSpPv30U6e29evXq2HDhipZsmSGy3h6esrT07MwygMAAEWAS4/c/PHHH9q/f7/2798v6cat3vv371dcXJykG0dd+vfv7+g/fPhwnThxQmPHjtXBgwe1YMECzZ8/X+PHj3dF+QAAoAhy6ZGbPXv2qHXr1o7nY8eOlSQNGDBA0dHRio+PdwQdSQoLC9PatWs1ZswYvfPOOwoKCtKbb77JbeAAAMDBZtKuyL1LXLp0Sb6+vrp48aJ8fHxcXQ4AAOkUxIW+GV3kW5wuKM7N/9/F6oJiAACA7BBuAACApRBuAACApRBuAACApRSrz7kBAMBVCusiX9w5jtwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLIdwAAABLcXd1AQAA5FXohDX5PubxGV3yfUwULo7cAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAAS+ETigEA+YpPDYarceQGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYirurCwAAFLzQCWsKZNzjM7oUyLjAneDIDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBSXh5vZs2crLCxMdrtdERER2rZtW5b9Fy9erHr16ql06dIKDAzUoEGDdO7cuUKqFgAAFHUuDTfLli3T6NGjNWnSJMXGxqpFixbq1KmT4uLiMuy/fft29e/fX0OGDNGBAwe0fPly7d69W0OHDi3kygEAQFHl0nATFRWlIUOGaOjQoQoPD9esWbNUtWpVzZkzJ8P+u3btUmhoqJ577jmFhYXpoYce0rBhw7Rnz55CrhwAABRVLgs3V69e1d69exUZGenUHhkZqR07dmS4TLNmzfTLL79o7dq1Msbo119/1ccff6wuXTL/EKnk5GRdunTJ6QEAAKzLZeHm7Nmzun79uvz9/Z3a/f39lZCQkOEyzZo10+LFi9W7d295eHgoICBAZcuW1VtvvZXpeqZPny5fX1/Ho2rVqvm6HQAAoGhx+QXFNpvN6bkxJl1bmh9//FHPPfecJk+erL179+qLL77QsWPHNHz48EzHnzhxoi5evOh4nDx5Ml/rBwAARYvLvluqYsWKKlGiRLqjNGfOnEl3NCfN9OnT1bx5c/31r3+VJN1///3y8vJSixYt9PLLLyswMDDdMp6envL09Mz/DQAAAEWSy47ceHh4KCIiQjExMU7tMTExatasWYbLJCUlyc3NueQSJUpIunHEBwAAwKWnpcaOHat58+ZpwYIFOnjwoMaMGaO4uDjHaaaJEyeqf//+jv5du3bVypUrNWfOHB09elRfffWVnnvuOTVq1EhBQUGu2gwAAFCEuOy0lCT17t1b586d07Rp0xQfH686depo7dq1CgkJkSTFx8c7febNwIEDlZiYqLffflvjxo1T2bJl1aZNG73yyiuu2gQAAFDEuDTcSNLIkSM1cuTIDF+Ljo5O1zZq1CiNGjWqgKsCAADFlcvvlgIAAMhPhBsAAGAphBsAAGApLr/mBgDuVqET1hTIuMdnZP6VNMDdgCM3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUgg3AADAUtxdXQAAFDWhE9bk+5jHZ3TJ9zEBZIwjNwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFL4hGIAxQKfGgwgpzhyAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALMXd1QUAKL5CJ6zJ9zGPz+iS72MCuLtw5AYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFiKy8PN7NmzFRYWJrvdroiICG3bti3L/snJyZo0aZJCQkLk6empatWqacGCBYVULQAAKOpceiv4smXLNHr0aM2ePVvNmzfXu+++q06dOunHH39UcHBwhsv06tVLv/76q+bPn6/q1avrzJkzSklJKeTKAQBAUeXScBMVFaUhQ4Zo6NChkqRZs2Zp3bp1mjNnjqZPn56u/xdffKEtW7bo6NGjKl++vCQpNDS0MEsGAABFnMtOS129elV79+5VZGSkU3tkZKR27NiR4TKffPKJGjZsqFdffVWVK1dWjRo1NH78eP3555+FUTIAACgGXHbk5uzZs7p+/br8/f2d2v39/ZWQkJDhMkePHtX27dtlt9u1atUqnT17ViNHjtT58+czve4mOTlZycnJjueXLl3Kv40AAABFjssvKLbZbE7PjTHp2tKkpqbKZrNp8eLFatSokTp37qyoqChFR0dnevRm+vTp8vX1dTyqVq2a79sAAACKDpeFm4oVK6pEiRLpjtKcOXMm3dGcNIGBgapcubJ8fX0dbeHh4TLG6JdffslwmYkTJ+rixYuOx8mTJ/NvIwAAQJHjsnDj4eGhiIgIxcTEOLXHxMSoWbNmGS7TvHlznT59Wn/88Yej7aeffpKbm5uqVKmS4TKenp7y8fFxegAAAOvKU7g5duxYvqx87NixmjdvnhYsWKCDBw9qzJgxiouL0/DhwyXdOOrSv39/R/++ffuqQoUKGjRokH788Udt3bpVf/3rXzV48GCVKlUqX2oCAADFW54uKK5evboefvhhDRkyRI8//rjsdnueVt67d2+dO3dO06ZNU3x8vOrUqaO1a9cqJCREkhQfH6+4uDhH/zJlyigmJkajRo1Sw4YNVaFCBfXq1Usvv/xyntYPAACsJ0/h5ttvv9WCBQs0btw4Pfvss+rdu7eGDBmiRo0a5XqskSNHauTIkRm+Fh0dna6tVq1a6U5lAbgpdMKafB/z+Iwu+T4mABSUPJ2WqlOnjqKionTq1CktXLhQCQkJeuihh1S7dm1FRUXpt99+y+86AQAAcuSOLih2d3dX9+7d9dFHH+mVV17RkSNHNH78eFWpUkX9+/dXfHx8ftUJAACQI3cUbvbs2aORI0cqMDBQUVFRGj9+vI4cOaKNGzfq1KlT6tatW37VCQAAkCN5uuYmKipKCxcu1KFDh9S5c2ctWrRInTt3lpvbjawUFhamd999V7Vq1crXYgEAALKTp3AzZ84cDR48WIMGDVJAQECGfYKDgzV//vw7Kg4AACC38hRuYmJiFBwc7DhSk8YYo5MnTyo4OFgeHh4aMGBAvhQJAACQU3m65qZatWo6e/Zsuvbz588rLCzsjosCAADIqzyFG2NMhu1//PFHnj/QDwAAID/k6rTU2LFjJd34Ju/JkyerdOnSjteuX7+ur7/+WvXr18/XAgEAAHIjV+EmNjZW0o0jN99//708PDwcr3l4eKhevXoaP358/lYIAACQC7kKN5s2bZIkDRo0SG+88QbfsA0AAIqcPN0ttXDhwvyuAwAAIF/kONz06NFD0dHR8vHxUY8ePbLsu3LlyjsuDAAAIC9yHG58fX1ls9kc/wYAACiKchxubj0VxWkpAABQVOXpc27+/PNPJSUlOZ6fOHFCs2bN0vr16/OtMAAAgLzIU7jp1q2bFi1aJEn6/fff1ahRI82cOVPdunXTnDlz8rVAAACA3MhTuNm3b59atGghSfr4448VEBCgEydOaNGiRXrzzTfztUAAAIDcyFO4SUpKkre3tyRp/fr16tGjh9zc3NSkSROdOHEiXwsEAADIjTx9zk316tW1evVqde/eXevWrdOYMWMkSWfOnOGD/YAMhE5YUyDjHp/RpUDGBYDiLE9HbiZPnqzx48crNDRUjRs3VtOmTSXdOIrzwAMP5GuBAAAAuZGnIzePP/64HnroIcXHx6tevXqO9rZt26p79+75VhwAAEBu5SncSFJAQIACAgKc2ho1anTHBQEAANyJPIWby5cva8aMGfryyy915swZpaamOr1+9OjRfCkOAAAgt/IUboYOHaotW7aoX79+CgwMdHwtAwAAgKvlKdx8/vnnWrNmjZo3b57f9QAAANyRPN0tVa5cOZUvXz6/awEAALhjeQo3f//73zV58mSn75cCAAAoCvJ0WmrmzJk6cuSI/P39FRoaqpIlSzq9vm/fvnwpDgAAILfyFG4ee+yxfC4DAAAgf+Qp3EyZMiW/6wAAAMgXebrmRpJ+//13zZs3TxMnTtT58+cl3TgdderUqXwrDgAAILfydOTmu+++U7t27eTr66vjx4/r6aefVvny5bVq1SqdOHFCixYtyu86AQAAciRPR27Gjh2rgQMH6vDhw7Lb7Y72Tp06aevWrflWHAAAQG7lKdzs3r1bw4YNS9deuXJlJSQk3HFRAAAAeZWncGO323Xp0qV07YcOHVKlSpXuuCgAAIC8ylO46datm6ZNm6Zr165Jkmw2m+Li4jRhwgT17NkzXwsEAADIjTyFm9dff12//fab/Pz89Oeff6ply5aqXr26vL299Y9//CO/awQAAMixPN0t5ePjo+3bt2vTpk3au3evUlNT1aBBA7Vr1y6/6wMAAMiVXIeb1NRURUdHa+XKlTp+/LhsNpvCwsIUEBAgY4xsNltB1AkAAJAjuTotZYzRo48+qqFDh+rUqVOqW7euateurRMnTmjgwIHq3r17QdUJAACQI7k6chMdHa2tW7fqyy+/VOvWrZ1e27hxox577DEtWrRI/fv3z9ciAQAAcipXR26WLl2qF198MV2wkaQ2bdpowoQJWrx4cb4VBwAAkFu5CjffffedOnbsmOnrnTp10rfffnvHRQEAAORVrsLN+fPn5e/vn+nr/v7+unDhwh0XBQAAkFe5CjfXr1+Xu3vml+mUKFFCKSkpd1wUAABAXuXqgmJjjAYOHChPT88MX09OTs6XogAAAPIqV+FmwIAB2fbhTikAAOBKuQo3CxcuLKg6AAAA8kWevlsKAACgqCLcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAAS3F5uJk9e7bCwsJkt9sVERGhbdu25Wi5r776Su7u7qpfv37BFggAAIoVl4abZcuWafTo0Zo0aZJiY2PVokULderUSXFxcVkud/HiRfXv319t27YtpEoBAEBx4dJwExUVpSFDhmjo0KEKDw/XrFmzVLVqVc2ZMyfL5YYNG6a+ffuqadOmhVQpAAAoLlwWbq5evaq9e/cqMjLSqT0yMlI7duzIdLmFCxfqyJEjmjJlSo7Wk5ycrEuXLjk9AACAdbks3Jw9e1bXr1+Xv7+/U7u/v78SEhIyXObw4cOaMGGCFi9eLHd39xytZ/r06fL19XU8qlatese1AwCAosvlFxTbbDan58aYdG2SdP36dfXt21cvvfSSatSokePxJ06cqIsXLzoeJ0+evOOaAQBA0ZWzwx8FoGLFiipRokS6ozRnzpxJdzRHkhITE7Vnzx7Fxsbq2WeflSSlpqbKGCN3d3etX79ebdq0Sbecp6enPD09C2YjAABAkeOyIzceHh6KiIhQTEyMU3tMTIyaNWuWrr+Pj4++//577d+/3/EYPny4atasqf3796tx48aFVToAACjCXHbkRpLGjh2rfv36qWHDhmratKnee+89xcXFafjw4ZJunFI6deqUFi1aJDc3N9WpU8dpeT8/P9nt9nTtAADg7uXScNO7d2+dO3dO06ZNU3x8vOrUqaO1a9cqJCREkhQfH5/tZ94AAADcyqXhRpJGjhypkSNHZvhadHR0lstOnTpVU6dOzf+iAABAseXyu6UAAADyE+EGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYiss/xA9wldAJawpk3OMzuhTIuACAnOHIDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBTCDQAAsBR3VxcA3Cp0wpoCGff4jC4FMi4AoOjhyA0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUwg0AALAUd1cXgOIhdMKafB/z+Iwu+T4mAAAuP3Ize/ZshYWFyW63KyIiQtu2bcu078qVK9W+fXtVqlRJPj4+atq0qdatW1eI1QIAgKLOpeFm2bJlGj16tCZNmqTY2Fi1aNFCnTp1UlxcXIb9t27dqvbt22vt2rXau3evWrdura5duyo2NraQKwcAAEWVS8NNVFSUhgwZoqFDhyo8PFyzZs1S1apVNWfOnAz7z5o1S//7v/+rBx98UPfee6/++c9/6t5779Wnn35ayJUDAICiymXh5urVq9q7d68iIyOd2iMjI7Vjx44cjZGamqrExESVL18+0z7Jycm6dOmS0wMAAFiXy8LN2bNndf36dfn7+zu1+/v7KyEhIUdjzJw5U5cvX1avXr0y7TN9+nT5+vo6HlWrVr2jugEAQNHm8guKbTab03NjTLq2jCxdulRTp07VsmXL5Ofnl2m/iRMn6uLFi47HyZMn77hmAABQdLnsVvCKFSuqRIkS6Y7SnDlzJt3RnNstW7ZMQ4YM0fLly9WuXbss+3p6esrT0/OO6wUAAMWDy47ceHh4KCIiQjExMU7tMTExatasWabLLV26VAMHDtSSJUvUpQufkwIAAJy59EP8xo4dq379+qlhw4Zq2rSp3nvvPcXFxWn48OGSbpxSOnXqlBYtWiTpRrDp37+/3njjDTVp0sRx1KdUqVLy9fV12XYAAICiw6Xhpnfv3jp37pymTZum+Ph41alTR2vXrlVISIgkKT4+3ukzb959912lpKTomWee0TPPPONoHzBggKKjowu7fAAAUAS5/OsXRo4cqZEjR2b42u2BZfPmzQVfEAAAKNZcfrcUAABAfiLcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAAS3F3dQHIu9AJa/J9zOMzuuT7mAAAFCaO3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEvhguJ8xkW+AAC4FkduAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApbg83MyePVthYWGy2+2KiIjQtm3bsuy/ZcsWRUREyG6365577tHcuXMLqVIAAFAcuDTcLFu2TKNHj9akSZMUGxurFi1aqFOnToqLi8uw/7Fjx9S5c2e1aNFCsbGxevHFF/Xcc89pxYoVhVw5AAAoqlwabqKiojRkyBANHTpU4eHhmjVrlqpWrao5c+Zk2H/u3LkKDg7WrFmzFB4erqFDh2rw4MF6/fXXC7lyAABQVLks3Fy9elV79+5VZGSkU3tkZKR27NiR4TI7d+5M179Dhw7as2ePrl27VmC1AgCA4sPdVSs+e/asrl+/Ln9/f6d2f39/JSQkZLhMQkJChv1TUlJ09uxZBQYGplsmOTlZycnJjucXL16UJF26dOlONyFDqclJ+T5mZrVacV0FsR7WVTDrYV13vq7iPi8Kc12ufq+suq7CnBf5NaYxJvvOxkVOnTplJJkdO3Y4tb/88sumZs2aGS5z7733mn/+859Obdu3bzeSTHx8fIbLTJkyxUjiwYMHDx48eFjgcfLkyWwzhsuO3FSsWFElSpRId5TmzJkz6Y7OpAkICMiwv7u7uypUqJDhMhMnTtTYsWMdz1NTU3X+/HlVqFBBNpvtDrci7y5duqSqVavq5MmT8vHxcVkdRQH74ib2xU3si5vYFzexL2662/aFMUaJiYkKCgrKtq/Lwo2Hh4ciIiIUExOj7t27O9pjYmLUrVu3DJdp2rSpPv30U6e29evXq2HDhipZsmSGy3h6esrT09OprWzZsndWfD7y8fG5KyZlTrAvbmJf3MS+uIl9cRP74qa7aV/4+vrmqJ9L75YaO3as5s2bpwULFujgwYMaM2aM4uLiNHz4cEk3jrr079/f0X/48OE6ceKExo4dq4MHD2rBggWaP3++xo8f76pNAAAARYzLjtxIUu/evXXu3DlNmzZN8fHxqlOnjtauXauQkBBJUnx8vNNn3oSFhWnt2rUaM2aM3nnnHQUFBenNN99Uz549XbUJAACgiHFpuJGkkSNHauTIkRm+Fh0dna6tZcuW2rdvXwFXVfA8PT01ZcqUdKfM7kbsi5vYFzexL25iX9zEvriJfZE5mzE5uacKAACgeHD5d0sBAADkJ8INAACwFMINAACwFMINAACwFMJNAZk9e7bCwsJkt9sVERGhbdu2Zdl/y5YtioiIkN1u1z333KO5c+cWUqUFa/r06XrwwQfl7e0tPz8/PfbYYzp06FCWy2zevFk2my3d47///W8hVV0wpk6dmm6bAgICslzGqvMiNDQ0w/f4mWeeybC/lebE1q1b1bVrVwUFBclms2n16tVOrxtjNHXqVAUFBalUqVJq1aqVDhw4kO24K1as0H333SdPT0/dd999WrVqVQFtQf7Jal9cu3ZNL7zwgurWrSsvLy8FBQWpf//+On36dJZjRkdHZzhXrly5UsBbc2eymxcDBw5Mt01NmjTJdtziOC/yA+GmACxbtkyjR4/WpEmTFBsbqxYtWqhTp05On9lzq2PHjqlz585q0aKFYmNj9eKLL+q5557TihUrCrny/LdlyxY988wz2rVrl2JiYpSSkqLIyEhdvnw522UPHTqk+Ph4x+Pee+8thIoLVu3atZ226fvvv8+0r5Xnxe7du532Q0xMjCTpiSeeyHI5K8yJy5cvq169enr77bczfP3VV19VVFSU3n77be3evVsBAQFq3769EhMTMx1z586d6t27t/r166dvv/1W/fr1U69evfT1118X1Gbki6z2RVJSkvbt26e//e1v2rdvn1auXKmffvpJjz76aLbj+vj4OM2T+Ph42e32gtiEfJPdvJCkjh07Om3T2rVrsxyzuM6LfJGD77hELjVq1MgMHz7cqa1WrVpmwoQJGfb/3//9X1OrVi2ntmHDhpkmTZoUWI2ucubMGSPJbNmyJdM+mzZtMpLMhQsXCq+wQjBlyhRTr169HPe/m+bF888/b6pVq2ZSU1MzfN2qc0KSWbVqleN5amqqCQgIMDNmzHC0Xblyxfj6+pq5c+dmOk6vXr1Mx44dndo6dOhgnnzyyXyvuaDcvi8y8s033xhJ5sSJE5n2WbhwofH19c3f4gpZRvtiwIABplu3brkaxwrzIq84cpPPrl69qr179yoyMtKpPTIyUjt27MhwmZ07d6br36FDB+3Zs0fXrl0rsFpd4eLFi5Kk8uXLZ9v3gQceUGBgoNq2batNmzYVdGmF4vDhwwoKClJYWJiefPJJHT16NNO+d8u8uHr1qj744AMNHjw42y+zteKcuNWxY8eUkJDg9L57enqqZcuWmf7+kDKfK1ktUxxdvHhRNpst2+8H/OOPPxQSEqIqVarokUceUWxsbOEUWMA2b94sPz8/1ahRQ08//bTOnDmTZf+7ZV5khHCTz86ePavr16+n+2Zzf3//dN9oniYhISHD/ikpKTp79myB1VrYjDEaO3asHnroIdWpUyfTfoGBgXrvvfe0YsUKrVy5UjVr1lTbtm21devWQqw2/zVu3FiLFi3SunXr9O9//1sJCQlq1qyZzp07l2H/u2VerF69Wr///rsGDhyYaR+rzonbpf2OyM3vj7TlcrtMcXPlyhVNmDBBffv2zfJLImvVqqXo6Gh98sknWrp0qex2u5o3b67Dhw8XYrX5r1OnTlq8eLE2btyomTNnavfu3WrTpo2Sk5MzXeZumBeZcfnXL1jV7X+BGmOy/Ks0o/4ZtRdnzz77rL777jtt3749y341a9ZUzZo1Hc+bNm2qkydP6vXXX9fDDz9c0GUWmE6dOjn+XbduXTVt2lTVqlXT+++/r7Fjx2a4zN0wL+bPn69OnTopKCgo0z5WnROZye3vj7wuU1xcu3ZNTz75pFJTUzV79uws+zZp0sTpQtvmzZurQYMGeuutt/Tmm28WdKkFpnfv3o5/16lTRw0bNlRISIjWrFmjHj16ZLqcledFVjhyk88qVqyoEiVKpEvGZ86cSZeg0wQEBGTY393dXRUqVCiwWgvTqFGj9Mknn2jTpk2qUqVKrpdv0qRJsf/L63ZeXl6qW7duptt1N8yLEydOaMOGDRo6dGiul7XinEi7ey43vz/SlsvtMsXFtWvX1KtXLx07dkwxMTFZHrXJiJubmx588EHLzZXAwECFhIRkuV1WnhfZIdzkMw8PD0VERDju/kgTExOjZs2aZbhM06ZN0/Vfv369GjZsqJIlSxZYrYXBGKNnn31WK1eu1MaNGxUWFpancWJjYxUYGJjP1blWcnKyDh48mOl2WXlepFm4cKH8/PzUpUuXXC9rxTkRFhamgIAAp/f96tWr2rJlS6a/P6TM50pWyxQHacHm8OHD2rBhQ55CvTFG+/fvt9xcOXfunE6ePJnldll1XuSIyy5ltrAPP/zQlCxZ0syfP9/8+OOPZvTo0cbLy8scP37cGGPMhAkTTL9+/Rz9jx49akqXLm3GjBljfvzxRzN//nxTsmRJ8/HHH7tqE/LNiBEjjK+vr9m8ebOJj493PJKSkhx9bt8f//rXv8yqVavMTz/9ZH744QczYcIEI8msWLHCFZuQb8aNG2c2b95sjh49anbt2mUeeeQR4+3tfVfOC2OMuX79ugkODjYvvPBCutesPCcSExNNbGysiY2NNZJMVFSUiY2NddwBNGPGDOPr62tWrlxpvv/+e9OnTx8TGBhoLl265BijX79+TndffvXVV6ZEiRJmxowZ5uDBg2bGjBnG3d3d7Nq1q9C3Lzey2hfXrl0zjz76qKlSpYrZv3+/0++P5ORkxxi374upU6eaL774whw5csTExsaaQYMGGXd3d/P111+7YhNzLKt9kZiYaMaNG2d27Nhhjh07ZjZt2mSaNm1qKleubMl5kR8INwXknXfeMSEhIcbDw8M0aNDA6dbnAQMGmJYtWzr137x5s3nggQeMh4eHCQ0NNXPmzCnkiguGpAwfCxcudPS5fX+88sorplq1asZut5ty5cqZhx56yKxZs6bwi89nvXv3NoGBgaZkyZImKCjI9OjRwxw4cMDx+t00L4wxZt26dUaSOXToULrXrDwn0m5rv/0xYMAAY8yN28GnTJliAgICjKenp3n44YfN999/7zRGy5YtHf3TLF++3NSsWdOULFnS1KpVq1gEv6z2xbFjxzL9/bFp0ybHGLfvi9GjR5vg4GDj4eFhKlWqZCIjI82OHTsKf+NyKat9kZSUZCIjI02lSpVMyZIlTXBwsBkwYICJi4tzGsMq8yI/2Iz5/1coAgAAWADX3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3ACwhFatWmn06NGuLgNAEUC4AeByXbt2Vbt27TJ8befOnbLZbNq3b18hVwWguCLcAHC5IUOGaOPGjTpx4kS61xYsWKD69eurQYMGLqgMQHFEuAHgco888oj8/PwUHR3t1J6UlKRly5bpscceU58+fVSlShWVLl1adevW1dKlS7Mc02azafXq1U5tZcuWdVrHqVOn1Lt3b5UrV04VKlRQt27ddPz48fzZKAAuQ7gB4HLu7u7q37+/oqOjdevX3S1fvlxXr17V0KFDFRERoc8++0w//PCD/vKXv6hfv376+uuv87zOpKQktW7dWmXKlNHWrVu1fft2lSlTRh07dtTVq1fzY7MAuAjhBkCRMHjwYB0/flybN292tC1YsEA9evRQ5cqVNX78eNWvX1/33HOPRo0apQ4dOmj58uV5Xt+HH34oNzc3zZs3T3Xr1lV4eLgWLlyouLg4pxoAFD/uri4AACSpVq1aatasmRYsWKDWrVvryJEj2rZtm9avX6/r169rxowZWrZsmU6dOqXk5GQlJyfLy8srz+vbu3evfv75Z3l7ezu1X7lyRUeOHLnTzQHgQoQbAEXGkCFD9Oyzz+qdd97RwoULFRISorZt2+q1117Tv/71L82aNUt169aVl5eXRo8eneXpI5vN5nSKS5KuXbvm+HdqaqoiIiK0ePHidMtWqlQp/zYKQKEj3AAoMnr16qXnn39eS5Ys0fvvv6+nn35aNptN27ZtU7du3fQ///M/km4Ek8OHDys8PDzTsSpVqqT4+HjH88OHDyspKcnxvEGDBlq2bJn8/Pzk4+NTcBsFoNBxzQ2AIqNMmTLq3bu3XnzxRZ0+fVoDBw6UJFWvXl0xMTHasWOHDh48qGHDhikhISHLsdq0aaO3335b+/bt0549ezR8+HCVLFnS8fpTTz2lihUrqlu3btq2bZuOHTumLVu26Pnnn9cvv/xSkJsJoIARbgAUKUOGDNGFCxfUrl07BQcHS5L+9re/qUGDBurQoYNatWqlgIAAPfbYY1mOM3PmTFWtWlUPP/yw+vbtq/Hjx6t06dKO10uXLq2tW7cqODhYPXr0UHh4uAYPHqw///yTIzlAMWczt5+UBgAAKMY4cgMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACyFcAMAACzl/wGqgpHIlUKwZgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(data, bins=range(min(data), max(data) + 1), align='left', rwidth=0.8, cumulative=True, density=True)\n",
    "\n",
    "# Add labels and title\n",
    "plt.xlabel('Value')\n",
    "plt.ylabel('Density')\n",
    "plt.title('n_old + n_new heavy atom mapping distribution')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5949094f-add9-4bd7-baa0-e6a77306be54",
   "metadata": {},
   "source": [
    "this is the distribution of the sum of old/new heavy atom changes for each transform;\n",
    "lets randomly sample 10 w/o replacement below and above the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "20cf025a-1a96-4b95-bfb0-2734f03e42cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = np.mean(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "761b84e5-2fec-468c-b648-68c7bf00c5cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "belows = [key for key, val in n_olds_news_map.items() if val < mean]\n",
    "aboves = [key for key, val in n_olds_news_map.items() if val > mean]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "b8cef1f1-a00e-41d7-90aa-e4d9f8f4c9cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "41"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(belows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "a30027e3-08bc-4cab-beaa-9b07757dcdba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "29"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(aboves)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "1771add3-011f-4136-95e8-2c9a53333287",
   "metadata": {},
   "outputs": [],
   "source": [
    "n=10\n",
    "sample_below_idxs = np.random.choice(np.arange(len(belows)), size=n, replace=False)\n",
    "sample_above_idxs =  np.random.choice(np.arange(len(aboves)), size=n, replace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "70f9253e-15c1-4a5d-ae4a-ab5dd3fcf0a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_belows = [belows[i] for i in sample_below_idxs]\n",
    "sample_aboves = [aboves[i] for i in sample_above_idxs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "99770d18-f75a-453a-bdb1-473aea48a1c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "my_samples = sample_belows + sample_aboves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "1d29a22c-154d-4eda-afbf-8110fdca85a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('CHEMBL3645036', 'CHEMBL3645054'),\n",
       " ('CHEMBL3642285', 'CHEMBL3642290'),\n",
       " ('CHEMBL3645002', 'CHEMBL3645060'),\n",
       " ('CHEMBL3644981', 'CHEMBL3645026'),\n",
       " ('CHEMBL3645082', 'CHEMBL3645092'),\n",
       " ('CHEMBL3645032', 'CHEMBL3645041'),\n",
       " ('CHEMBL3645079', 'CHEMBL3645091'),\n",
       " ('CHEMBL3642290', 'CHEMBL3642297'),\n",
       " ('CHEMBL3645071', 'CHEMBL3645074'),\n",
       " ('CHEMBL3645079', 'CHEMBL3645085'),\n",
       " ('CHEMBL3639411', 'CHEMBL3642281'),\n",
       " ('CHEMBL3645054', 'CHEMBL3645092'),\n",
       " ('CHEMBL3642294', 'CHEMBL3645060'),\n",
       " ('CHEMBL3642350', 'CHEMBL3645002'),\n",
       " ('CHEMBL3644991', 'CHEMBL3645054'),\n",
       " ('CHEMBL3644987', 'CHEMBL3645054'),\n",
       " ('CHEMBL3645035', 'CHEMBL3645054'),\n",
       " ('CHEMBL3642350', 'CHEMBL3645060'),\n",
       " ('CHEMBL3642297', 'CHEMBL3642389'),\n",
       " ('CHEMBL3644987', 'CHEMBL3645051')]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "cdc690ca-497f-430d-885a-bcea5a468b5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the file path\n",
    "file_path = \"sample_edges.txt\"\n",
    "\n",
    "# Open the file in write mode\n",
    "with open(file_path, 'w') as file:\n",
    "    # Iterate over the list of tuples\n",
    "    for item in my_samples:\n",
    "        # Convert the tuple to a string and write it to the file\n",
    "        file.write(f\"{item[0]} {item[1]}\" + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53ce802c-a587-477a-89eb-9bdb0c455b60",
   "metadata": {},
   "source": [
    "after running, many of these jobs failed. I am going to tail the list of `.err` files to query which ones are bad.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "260fbc84-7178-42d1-91a4-d9a4c3dd05d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "dda23601-0abe-4507-827f-3dc1db31c5fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/data1/choderaj/rufad/tm'"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "4de8fd52-733a-4331-b09e-39f3e5ab622c",
   "metadata": {},
   "outputs": [],
   "source": [
    "run_dir = os.path.join(os.getcwd(), 'jak2_esp1', 'temp_fails')\n",
    "my_file_list = [os.path.join(run_dir, f) for f in os.listdir(run_dir) if f[-3:] == 'err']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "4f416ee8-26e2-43c3-9231-d1732636cdf8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_19_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_14_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_8_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_3_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_4_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_12_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_17_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_10_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_2_tm-jak2_esp1_isce002.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_15_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_13_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_9_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_20_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_21_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_18_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_5_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_7_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_16_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_11_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_6_tm-jak2_esp1_iscb025.err',\n",
       " '/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_1_tm-jak2_esp1_isce002.err']"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_file_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "e8c40b22-8316-42d2-8b89-428e446c0580",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_19_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_14_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_8_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_3_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_4_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_12_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_17_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_10_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_2_tm-jak2_esp1_isce002.err:\n",
      "slurmstepd: error: *** JOB 332390 ON isce002 CANCELLED AT 2024-04-30T14:43:08 ***\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_15_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_13_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_9_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_20_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_21_tm-jak2_esp1_iscb025.err:\n",
      "IndexError: list index out of range\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_18_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_5_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_7_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_16_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_11_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_6_tm-jak2_esp1_iscb025.err:\n",
      "GPUassert: no kernel image is available for execution on the device /home/rufad/github2/timemachine/timemachine/cpp/src/neighborlist.cu 289\n",
      "\n",
      "/data1/choderaj/rufad/tm/jak2_esp1/temp_fails/332387_1_tm-jak2_esp1_isce002.err:\n",
      "slurmstepd: error: *** JOB 332389 ON isce002 CANCELLED AT 2024-04-30T14:43:08 ***\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def tail_files(file_list, n):\n",
    "    for file_name in file_list:\n",
    "        try:\n",
    "            with open(file_name, 'r') as file:\n",
    "                lines = file.readlines()\n",
    "                last_n_lines = lines[-n:]\n",
    "                print(f\"{file_name}:\")\n",
    "                for line in last_n_lines:\n",
    "                    print(line, end='')\n",
    "        except FileNotFoundError:\n",
    "            print(f\"File '{file_name}' not found.\")\n",
    "        except IOError:\n",
    "            print(f\"Error reading file '{file_name}'.\")\n",
    "        print()\n",
    "\n",
    "# Example usage:\n",
    "n = 1  # Number of lines to tail\n",
    "tail_files(my_file_list, n)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eef57263-6902-450e-947e-38ff54f16599",
   "metadata": {},
   "source": [
    "at this point, I am ready to save _just_ the ligand coords...let's see if we can't modify an existing array for just ligand coordinates..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "895387f6-8b8f-4b05-b255-b1723d91b5d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from timemachine.fe.stored_arrays import StoredArrays\n",
    "from timemachine.fe.free_energy import Trajectory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77e4820d-8a8a-4d6c-a6d4-4f14ece652e8",
   "metadata": {},
   "source": [
    "it might be valuable to just write .npz array objects for ligand positions and blank out trajectories since stored arrays need access to memory every time we want the data. then it should be safe to empty the tmpdir."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34a8796f-5742-4e41-9f64-26495f1848c4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "timemachine_off",
   "language": "python",
   "name": "timemachine_off"
  },
  "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}