notebooks/Benchmark_Ursu_Example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8c743d91",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import pickle\n",
    "from scipy import sparse\n",
    "import scanpy as sc\n",
    "import anndata as ad\n",
    "from anndata import AnnData\n",
    "import scvi\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import torch\n",
    "import torch.nn.functional as F\n",
    "import torch.nn as nn\n",
    "\n",
    "from perturbnet.util import * \n",
    "from perturbnet.cinn.flow import * \n",
    "from perturbnet.chemicalvae.chemicalVAE import *\n",
    "from perturbnet.data_vae.vae import *\n",
    "from perturbnet.cinn.flow_generate import TFVAEZ_CheckNet2Net, SCVIZ_CheckNet2Net\n",
    "import pickle5\n",
    "from tqdm import tqdm \n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "from perturbnet.baselines.knn import *\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bed0670",
   "metadata": {},
   "source": [
    "# Load Data and PerturbNet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b1407c30",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(156632, 17)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(133231, 17)\n",
      "\u001b[34mINFO    \u001b[0m No batch_key inputted, assuming all cells are same batch                                                  \n",
      "\u001b[34mINFO    \u001b[0m No label_key inputted, assuming all cells have same label                                                 \n",
      "\u001b[34mINFO    \u001b[0m Using data from adata.layers\u001b[1m[\u001b[0m\u001b[32m\"counts\"\u001b[0m\u001b[1m]\u001b[0m                                                                    \n",
      "\u001b[34mINFO    \u001b[0m Computing library size prior per batch                                                                    \n",
      "\u001b[34mINFO    \u001b[0m Successfully registered anndata object containing \u001b[1;36m123653\u001b[0m cells, \u001b[1;36m1629\u001b[0m vars, \u001b[1;36m1\u001b[0m batches, \u001b[1;36m1\u001b[0m labels, and \u001b[1;36m0\u001b[0m     \n",
      "         proteins. Also registered \u001b[1;36m0\u001b[0m extra categorical covariates and \u001b[1;36m0\u001b[0m extra continuous covariates.               \n",
      "\u001b[34mINFO    \u001b[0m Please do not further modify adata until model is trained.                                                \n",
      "\u001b[34mINFO    \u001b[0m Using data from adata.layers\u001b[1m[\u001b[0m\u001b[32m\"counts\"\u001b[0m\u001b[1m]\u001b[0m                                                                    \n",
      "\u001b[34mINFO    \u001b[0m Computing library size prior per batch                                                                    \n",
      "\u001b[34mINFO    \u001b[0m Registered keys:\u001b[1m[\u001b[0m\u001b[32m'X'\u001b[0m, \u001b[32m'batch_indices'\u001b[0m, \u001b[32m'local_l_mean'\u001b[0m, \u001b[32m'local_l_var'\u001b[0m, \u001b[32m'labels'\u001b[0m\u001b[1m]\u001b[0m                           \n",
      "\u001b[34mINFO    \u001b[0m Successfully registered anndata object containing \u001b[1;36m123653\u001b[0m cells, \u001b[1;36m1629\u001b[0m vars, \u001b[1;36m1\u001b[0m batches, \u001b[1;36m1\u001b[0m labels, and \u001b[1;36m0\u001b[0m     \n",
      "         proteins. Also registered \u001b[1;36m0\u001b[0m extra categorical covariates and \u001b[1;36m0\u001b[0m extra continuous covariates.               \n",
      "Note: Conditioning network uses batch-normalization. Make sure to train with a sufficiently large batch size\n"
     ]
    }
   ],
   "source": [
    "some_path = \"\"\n",
    "path_cinn_model =  some_path + '/cinn/50ep/'\n",
    "path_data = \"\"\n",
    "path_scvi_model_train = some_path + '/cellvae/700ep/model/' \n",
    "\n",
    "    \n",
    "adata = ad.read_h5ad(os.path.join(path_data, \"ursu_400filter_standard.h5ad\"))\n",
    "\n",
    "adata_save_index = np.array(adata.obs[adata.obs[\"variant_seq_use\"] == \"Yes\"].index).astype(int)\n",
    "adata = adata[adata_save_index, :].copy()    \n",
    "print(adata.obs.shape)\n",
    "adata.obs = adata.obs.reset_index(drop = True)\n",
    "\n",
    "adata_save_index = np.array(adata.obs[adata.obs[\"if_control\"] != \"Wild Type Like\"].index).astype(int)\n",
    "adata = adata[adata_save_index, :].copy()    \n",
    "print(adata.obs.shape)\n",
    "adata.obs = adata.obs.reset_index(drop = True)\n",
    "    \n",
    "#scvi.data.setup_anndata(adata, layer = \"counts\")\n",
    "adata_test_index = np.array(adata.obs[adata.obs[\"1st_hengshi_holdout\"] == \"test\"].index).astype(int)\n",
    "adata_test = adata[adata_test_index, :].copy()\n",
    "    \n",
    "adata_train_index = np.array(adata.obs[adata.obs[\"1st_hengshi_holdout\"] == \"train\"].index).astype(int)\n",
    "adata_train = adata[adata_train_index, :].copy()\n",
    "\n",
    "scvi.data.setup_anndata(adata_train, layer = \"counts\")\n",
    "scvi_model_cinn = scvi.model.SCVI.load(path_scvi_model_train, adata_train, use_cuda = False)\n",
    "    \n",
    "perturb_with_onehot_kept, perturb_with_onehot_removed = np.array(adata_train.obs[\"variant_seq\"]), np.array(adata_test.obs[\"variant_seq\"])\n",
    "trt_obs_list, trt_unseen_list = np.unique(perturb_with_onehot_kept), np.unique(perturb_with_onehot_removed)    \n",
    "\n",
    "onehot_unseen_idx = []\n",
    "for i in range(len(trt_unseen_list)):\n",
    "    trt = trt_unseen_list[i]\n",
    "    onehot_idx = np.where(adata.uns[\"ordered_all_trt\"] == trt)[0][0]\n",
    "    \n",
    "    onehot_unseen_idx.append(onehot_idx)\n",
    "    unseen_onehot = adata.uns[\"ordered_all_embedding\"][onehot_unseen_idx]\n",
    "\n",
    "onehot_seen_idx = []\n",
    "for i in range(len(trt_obs_list)):\n",
    "    trt = trt_obs_list[i]\n",
    "    onehot_idx = np.where(adata.uns[\"ordered_all_trt\"] == trt)[0][0]\n",
    "    onehot_seen_idx.append(onehot_idx)\n",
    "    seen_onehot = adata.uns[\"ordered_all_embedding\"][onehot_seen_idx]\n",
    "\n",
    "device = 'cuda' if torch.cuda.is_available() else 'cpu'\n",
    "\n",
    "perturbToEmbed = {}\n",
    "for i in range(trt_obs_list.shape[0]):\n",
    "    perturbToEmbed[trt_obs_list[i]] = i\n",
    "        \n",
    "torch.manual_seed(42)\n",
    "flow_model = ConditionalFlatCouplingFlow(conditioning_dim=1280,embedding_dim=10,conditioning_depth=2,n_flows=20,in_channels=10,hidden_dim=1024,\n",
    "hidden_depth=2,activation=\"none\",conditioner_use_bn=True)\n",
    "\n",
    "model_c = Net2NetFlow_scVIFixFlow(configured_flow = flow_model,\n",
    "                                 cond_stage_data = perturb_with_onehot_kept,\n",
    "                                perturbToEmbedLib = perturbToEmbed,\n",
    "                                  embedData = seen_onehot,\n",
    "                                     scvi_model = scvi_model_cinn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d6b677fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_c.to(device = device)\n",
    "model_c.load(path_cinn_model)\n",
    "model_c.eval()\n",
    "scvi_model_de = scvi_predictive_z(scvi_model_cinn)\n",
    "perturbnet_model = SCVIZ_CheckNet2Net(model_c, device, scvi_model_de)\n",
    "if sparse.issparse(adata_test.X):\n",
    "    usedata = adata_test.X.A\n",
    "else:\n",
    "    usedata = adata_test.X\n",
    "\n",
    "if sparse.issparse(adata_test.layers[\"counts\"]):\n",
    "    usedata_count = adata_test.layers[\"counts\"].A\n",
    "else:\n",
    "    usedata_count = adata_test.layers[\"counts\"]\n",
    "\n",
    "fidscore_cal = fidscore()\n",
    "\n",
    "normModel = NormalizedRevisionRSquare(largeCountData = usedata_count)\n",
    "Lsample_obs = scvi_model_cinn.get_latent_library_size(adata = adata_train, give_mean = False)\n",
    "adata.var[\"gene_idx\"] = np.arange(0,adata.var.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f1291979",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 12/12 [00:03<00:00,  3.51it/s]\n"
     ]
    }
   ],
   "source": [
    "perturb = []\n",
    "ncell = []\n",
    "n_large = []\n",
    "r2 = []\n",
    "r2_deg = []\n",
    "r2_large = []\n",
    "fid_mu = []\n",
    "fid_sample = []\n",
    "pear = []\n",
    "pear_deg = []\n",
    "pear_large = []\n",
    "\n",
    "np.random.seed(2023)\n",
    "for indice_trt in tqdm(range(len(trt_unseen_list))):\n",
    "    trt_type = trt_unseen_list[indice_trt]\n",
    "    pert = trt_type\n",
    "    idx_trt_type = np.where(perturb_with_onehot_removed == trt_type)[0]\n",
    "    if idx_trt_type.shape[0] > 1000:\n",
    "        idx_trt_type  =  np.random.choice(idx_trt_type, 1000, replace = False)\n",
    "    Lsample_idx = np.random.choice(range(Lsample_obs.shape[0]), 500, replace=True)\n",
    "    DEG_gene = adata.uns[\"rank_genes_groups\"][\"names\"][trt_type]\n",
    "    DEG_idx = np.array(adata.var.loc[DEG_gene][\"gene_idx\"])\n",
    "\n",
    "    onehot_indice_trt = np.tile(unseen_onehot[[indice_trt]], (500, 1))\n",
    "\n",
    "\n",
    "\n",
    "    trt_onehot = onehot_indice_trt + np.random.normal(scale = 0.001, size = onehot_indice_trt.shape)\n",
    "    library_trt_latent =Lsample_obs[Lsample_idx]\n",
    "\n",
    "    _, fake_data = perturbnet_model.sample_data(trt_onehot, library_trt_latent)\n",
    "    real_data = usedata_count[idx_trt_type]\n",
    "    r2_value, real_norm, rfake_norm = normModel.calculate_r_square(real_data, fake_data)\n",
    "    pear_value = normModel.calculate_pearson(real_data, fake_data)\n",
    "\n",
    "        \n",
    "\n",
    "    real_data_deg = real_data[:,DEG_idx]\n",
    "    fake_data_deg = fake_data[:,DEG_idx]\n",
    "\n",
    "\n",
    "    r2_deg_value,_,_ = normModel.calculate_r_square(real_data_deg, fake_data_deg)\n",
    "    pear_deg_value = normModel.calculate_pearson(real_data_deg, fake_data_deg)\n",
    "\n",
    "   \n",
    "    large_effect_idx =   DEG_idx[abs(adata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][pert])>=1]\n",
    "    num_large =  sum(abs(adata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][pert])>=1)\n",
    "    \n",
    "    if num_large <=1:\n",
    "        hd_large_value = 1.5\n",
    "        r2_large_value = 1.5\n",
    "        pear_large_value = 1.5\n",
    "    else:\n",
    "        real_data_large = real_data[:,large_effect_idx ]\n",
    "        fake_data_large = fake_data[:,large_effect_idx ]\n",
    "\n",
    "        r2_large_value,_ ,_  = normModel.calculate_r_square(real_data_large, fake_data_large)\n",
    "        pear_large_value = normModel.calculate_pearson(real_data_large, fake_data_large)\n",
    "        \n",
    "\n",
    "    perturb.append(pert)\n",
    "    ncell.append(len(idx_trt_type))\n",
    "    n_large.append(num_large)\n",
    "    r2.append(r2_value)\n",
    "    r2_deg.append(r2_deg_value)\n",
    "    r2_large.append(r2_large_value)\n",
    "    pear.append(pear_value)\n",
    "    pear_deg.append(pear_deg_value)\n",
    "    pear_large.append(pear_large_value)\n",
    "\n",
    "results_PerturbNet = pd.DataFrame({\"perturbation\":perturb, \"number_real_cells_used\":ncell,\"n_large\":n_large,\n",
    "                        \"r2\":r2,\"pear\":pear,\"r2_deg\":r2_deg, \"pear_deg\":pear_deg,\n",
    "                        \"r2_large\":r2_large, \"pear_large\": pear_large\n",
    "                       })"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c07e17",
   "metadata": {},
   "source": [
    "# Training Sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "18c320f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/12 [00:00<?, ?it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      "  8%|▊         | 1/12 [00:00<00:06,  1.62it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 17%|█▋        | 2/12 [00:01<00:06,  1.56it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 25%|██▌       | 3/12 [00:01<00:05,  1.55it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 33%|███▎      | 4/12 [00:02<00:05,  1.55it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 42%|████▏     | 5/12 [00:03<00:04,  1.54it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 50%|█████     | 6/12 [00:03<00:03,  1.58it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 58%|█████▊    | 7/12 [00:04<00:03,  1.59it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 67%|██████▋   | 8/12 [00:05<00:02,  1.59it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 75%|███████▌  | 9/12 [00:05<00:01,  1.58it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 83%|████████▎ | 10/12 [00:06<00:01,  1.59it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 92%|█████████▏| 11/12 [00:06<00:00,  1.57it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      "100%|██████████| 12/12 [00:07<00:00,  1.58it/s]\n"
     ]
    }
   ],
   "source": [
    "adata.var[\"gene_idx\"] = np.arange(0,adata.var.shape[0],1)\n",
    "\n",
    "perturb = []\n",
    "ncell = []\n",
    "n_large = []\n",
    "r2 = []\n",
    "r2_deg = []\n",
    "r2_large = []\n",
    "fid_mu = []\n",
    "fid_sample = []\n",
    "pear = []\n",
    "pear_deg = []\n",
    "pear_large = []\n",
    "\n",
    "np.random.seed(2023)\n",
    "for indice_trt in tqdm(range(len(trt_unseen_list))):\n",
    "    pert = trt_unseen_list[indice_trt]\n",
    "    seen_data = adata_train.layers[\"counts\"].A\n",
    "    seen_data_idx = list(range(seen_data.shape[0]))\n",
    "\n",
    "    real_data_idx = adata.obs[adata.obs[\"variant_seq\"] == pert].index\n",
    "    if real_data_idx.shape[0] > 1000:\n",
    "        real_data_idx = np.random.choice(real_data_idx, 1000, replace = False)\n",
    "    \n",
    "    real_data = adata[real_data_idx, :].copy().layers[\"counts\"].A\n",
    "    norm_vec_real = adata[real_data_idx, :].obs[\"n_counts_total\"].to_numpy()\n",
    "    idx_rsample = np.random.choice(seen_data_idx, len(real_data_idx), replace=True)\n",
    "    fake_data = seen_data[idx_rsample]\n",
    "    norm_vec_fake = adata_train.obs[\"n_counts_total\"].to_numpy()[idx_rsample]\n",
    "    \n",
    "    r2_value = normModel.calculate_r_square_norm_by_all(real_data, fake_data,norm_vec_real,norm_vec_fake)\n",
    "    pear_value = normModel.calculate_pearson_norm_by_all(real_data, fake_data,norm_vec_real,norm_vec_fake)\n",
    "\n",
    "\n",
    "    \n",
    "    DEG_gene = adata.uns[\"rank_genes_groups\"][\"names\"][pert]\n",
    "    DEG_idx = np.array(adata.var.loc[DEG_gene][\"gene_idx\"])\n",
    "    \n",
    "    \n",
    "    real_data_deg = real_data[:,DEG_idx]\n",
    "    fake_data_deg = fake_data[:,DEG_idx]\n",
    "    r2_deg_value= normModel.calculate_r_square_norm_by_all(real_data_deg, fake_data_deg,norm_vec_real,norm_vec_fake)\n",
    "    pear_deg_value = normModel.calculate_pearson_norm_by_all(real_data_deg, fake_data_deg,norm_vec_real,norm_vec_fake)\n",
    "\n",
    "    \n",
    "    large_effect_idx =  DEG_idx[abs(adata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][pert])>=1]\n",
    "    num_large =  sum(abs(adata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][pert])>=1)\n",
    "    if num_large <=1:\n",
    "        hd_large_value = 1.5\n",
    "        r2_large_value = 1.5\n",
    "        pear_large_value = 1.5\n",
    "    else:\n",
    "        real_data_large = real_data[:,large_effect_idx ]\n",
    "        fake_data_large = fake_data[:,large_effect_idx ]\n",
    "\n",
    "        r2_large_value = normModel.calculate_r_square_norm_by_all(real_data_large, fake_data_large,norm_vec_real,norm_vec_fake)\n",
    "        pear_large_value = normModel.calculate_pearson_norm_by_all(real_data_large, fake_data_large,norm_vec_real,norm_vec_fake)\n",
    "\n",
    "    \n",
    "    perturb.append(pert)\n",
    "    ncell.append(real_data_idx.shape[0])\n",
    "    n_large.append(num_large)\n",
    "    r2.append(r2_value)\n",
    "    r2_deg.append(r2_deg_value)\n",
    "    r2_large.append(r2_large_value)\n",
    "    pear.append(pear_value)\n",
    "    pear_deg.append(pear_deg_value)\n",
    "    pear_large.append(pear_large_value)\n",
    "\n",
    "    \n",
    "results_sample = pd.DataFrame({\"perturbation\":perturb, \"number_real_cells_used\":ncell,\"n_large\":n_large,\n",
    "                        \"r2\":r2,\"pear\":pear,\"r2_deg\":r2_deg, \"pear_deg\":pear_deg,\n",
    "                        \"r2_large\":r2_large, \"pear_large\": pear_large\n",
    "                       })\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65997cd9",
   "metadata": {},
   "source": [
    "# KNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5853f5b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/12 [00:00<?, ?it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      "  8%|▊         | 1/12 [00:00<00:01,  6.61it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 17%|█▋        | 2/12 [00:00<00:01,  6.23it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 25%|██▌       | 3/12 [00:00<00:01,  6.26it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 33%|███▎      | 4/12 [00:00<00:01,  6.12it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 42%|████▏     | 5/12 [00:00<00:01,  6.19it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 50%|█████     | 6/12 [00:00<00:00,  6.87it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 58%|█████▊    | 7/12 [00:01<00:00,  7.39it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 67%|██████▋   | 8/12 [00:01<00:00,  7.18it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 75%|███████▌  | 9/12 [00:01<00:00,  7.14it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 83%|████████▎ | 10/12 [00:01<00:00,  7.21it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      " 92%|█████████▏| 11/12 [00:01<00:00,  6.82it/s]/home/wzqian/.local/lib/python3.7/site-packages/anndata/_core/anndata.py:121: ImplicitModificationWarning: Transforming to str index.\n",
      "100%|██████████| 12/12 [00:01<00:00,  6.83it/s]\n"
     ]
    }
   ],
   "source": [
    "neigh = NearestNeighbors(n_neighbors=5)\n",
    "neigh_fit = neigh.fit(seen_onehot)\n",
    "adata.var[\"gene_idx\"] = np.arange(0,adata.var.shape[0],1)\n",
    "\n",
    "perturb = []\n",
    "ncell = []\n",
    "n_large = []\n",
    "r2 = []\n",
    "r2_deg = []\n",
    "r2_large = []\n",
    "\n",
    "pear = []\n",
    "pear_deg = []\n",
    "pear_large = []\n",
    "np.random.seed(2023)\n",
    "seen_data = adata_train.layers[\"counts\"].A\n",
    "seen_data_idx = list(range(seen_data.shape[0]))\n",
    "\n",
    "for indice_trt in tqdm(range(len(trt_unseen_list))):\n",
    "    pert = trt_unseen_list[indice_trt]\n",
    "\n",
    "    #fake_data = pred\n",
    "    real_data_idx = adata.obs[adata.obs[\"variant_seq\"] == pert].index\n",
    "    if real_data_idx.shape[0] > 1000:\n",
    "        real_data_idx = np.random.choice(real_data_idx, 1000, replace = False)\n",
    "    \n",
    "    real_data = adata[real_data_idx, :].copy().layers[\"counts\"].A\n",
    "    norm_vec_real = adata[real_data_idx, :].obs[\"n_counts_total\"].to_numpy()\n",
    "    distances, other_trts = neigh_fit.kneighbors(unseen_onehot[[indice_trt]], 5, return_distance=True)\n",
    "    samplerNN = samplefromNeighbors(distances, other_trts)\n",
    "    idx_sample = samplerNN.samplingTrt(trt_obs_list, perturb_with_onehot_kept, len(real_data_idx))\n",
    "    \n",
    "    fake_data = seen_data[idx_sample]\n",
    "    norm_vec_fake = adata_train.obs[\"n_counts_total\"].to_numpy()[idx_sample]\n",
    "    \n",
    "    r2_value = normModel.calculate_r_square_norm_by_all(real_data, fake_data,norm_vec_real,norm_vec_fake)\n",
    "    pear_value = normModel.calculate_pearson_norm_by_all(real_data, fake_data,norm_vec_real,norm_vec_fake)\n",
    "\n",
    "\n",
    "    \n",
    "    DEG_gene = adata.uns[\"rank_genes_groups\"][\"names\"][pert]\n",
    "    DEG_idx = np.array(adata.var.loc[DEG_gene][\"gene_idx\"])\n",
    "    \n",
    "    \n",
    "    real_data_deg = real_data[:,DEG_idx]\n",
    "    fake_data_deg = fake_data[:,DEG_idx]\n",
    "    r2_deg_value= normModel.calculate_r_square_norm_by_all(real_data_deg, fake_data_deg,norm_vec_real,norm_vec_fake)\n",
    "    pear_deg_value = normModel.calculate_pearson_norm_by_all(real_data_deg, fake_data_deg,norm_vec_real,norm_vec_fake)\n",
    "\n",
    "    \n",
    "    large_effect_idx =  DEG_idx[abs(adata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][pert])>=1]\n",
    "    num_large =  sum(abs(adata.uns[\"rank_genes_groups\"][\"logfoldchanges\"][pert])>=1)\n",
    "    if num_large <=1:\n",
    "        r2_large_value = 1.5\n",
    "        pear_large_value = 1.5\n",
    "    else:\n",
    "        real_data_large = real_data[:,large_effect_idx ]\n",
    "        fake_data_large = fake_data[:,large_effect_idx ]\n",
    "        r2_large_value = normModel.calculate_r_square_norm_by_all(real_data_large, fake_data_large,norm_vec_real,norm_vec_fake)\n",
    "        pear_large_value = normModel.calculate_pearson_norm_by_all(real_data_large, fake_data_large,norm_vec_real,norm_vec_fake)\n",
    "    \n",
    "    \n",
    "    perturb.append(pert)\n",
    "    ncell.append(real_data_idx.shape[0])\n",
    "    n_large.append(num_large)\n",
    "    r2.append(r2_value)\n",
    "    r2_deg.append(r2_deg_value)\n",
    "    r2_large.append(r2_large_value)\n",
    "    pear.append(pear_value)\n",
    "    pear_deg.append(pear_deg_value)\n",
    "    pear_large.append(pear_large_value)\n",
    "\n",
    "    \n",
    "results_KNN = pd.DataFrame({\"perturbation\":perturb, \"number_real_cells_used\":ncell,\"n_large\":n_large,\n",
    "                        \"r2\":r2,\"pear\":pear,\"r2_deg\":r2_deg, \"pear_deg\":pear_deg,\n",
    "                        \"r2_large\":r2_large, \"pear_large\": pear_large\n",
    "                       })\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "30641538",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "../perturbnet/util.py:963: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuQAAAGuCAYAAADYjxPyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEKUlEQVR4nO3deVyVZf7/8ffhAAc4yCZCKIqWqaO2qVk53yZcytRGLTdySZwWrcmtUtNcwNRSK7WvZtSYK6ZmOqWjzlhhiy1mjWkzWuaGGpooirIIHM7vj36eb6So4DnnQng9Hw8fcN/3dd/354ZLeHOd676Pxel0OgUAAADACB/TBQAAAABVGYEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBBHIAAADAIF/TBVQUmZmZpkuoUiwWiwIDA5WXlyfemwqeRF+Dt9DXzIiMjDRdAnDFGCGHET4+PgoKCpKPD10QnkVfg7fQ1wCUFz81AAAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQb6mCzhn7dq1+uijj7R//37dcccdGjFiRKltO3fuLJvNJovFIklq3LixkpKSXNs3b96shQsX6sSJE2rUqJGGDBmiqKgoT18CAAAAUGYVJpBHRESoZ8+e2rZtm06fPn3J9jNmzFBsbOx56w8ePKhZs2Zp9OjRaty4sRYvXqxp06bppZde8kTZAAAAwBWpMFNWWrVqpdtvv10hISFXdJxNmzapWbNmuuWWW2Sz2dS7d2/t27dP6enpbqoUV2rTpk3q2rWrbrvtNnXv3l07duwwXRIqKfoavIW+BuBKVJhAXlZjx45Vv379NHHixBJh+8CBA6pXr55rOSgoSNdcc40OHDhgokz8ztatW/Xiiy8qJydHDodDp06d0siRI/mDCW5HX4O30NcAXKkKM2WlLKZMmaKGDRuqsLBQq1at0vjx4/Xaa68pKChI+fn5stvtJdrb7Xbl5eWVWJeRkaGMjAzXss1mU82aNb1Sf1W2cOHC89Y5nU4tXrxY48ePN1ARKiv6GryFvgbgSl2Vgbxp06aSJD8/P/Xt21dpaWnauXOnmjdvroCAAOXm5pZon5ubq8DAwBLrUlJSlJyc7FoeM2aMJk+e7Pniq7iCgoJS14eHh3u5GlRm9DV4C30NwJW6KgP575172ookxcXFae/eva7lvLw8HTlyRHFxcSX2GThwoDp37uxattlsysrK8nyxVVzTpk0vOH2oWbNmfP3hVvQ1eAt9zSz+6EFlUGECucPhkMPhUHFxsYqLi1VQUCAfHx/5+pYsMT09XYWFhapbt66Kior07rvvqqCgQA0bNpQkxcfH6+mnn9a2bdvUuHFjLV26VHXr1lWdOnVKHCcmJkYxMTGu5czMTDkcDs9faBX317/+Vf/973+1b98+17qWLVuqS5cufP3hVvQ1eAt9DcCVsjidTqfpIiRp6dKlWrZsWYl1bdq00bBhw9SzZ09NmDBBTZo00fbt2zV37lxlZmbK399f9evXV2JiYokbOT/77DMtXLhQWVlZatiwoYYOHXrJ55BnZmZ65LpwYf/+97+VkZGha6+9Vo0aNTJdDiox+hq8hb5mRmRkpOkSgCtWYQK5aQRy77JarQoPD1dWVhYjSPAo+hq8hb5mBoEclcFV+9hDAAAAoDIgkAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRwAAAAwiEAOAAAAGEQgBwAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBvqYLqCj8/f1ls9lMl1FlWCwWSZLdbpfT6TRcDSoz+hq8hb4GoLwI5P9fQUGBCgoKTJdRZVitVvn7+ysnJ0cOh8N0OajE6GvwFvqaGQymoTJgygoAAABgEIEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRwAAAAwiEAOAAAAGEQgBwAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADDI13QB56xdu1YfffSR9u/frzvuuEMjRoy4YLtdu3bp7bff1k8//SRJatiwoR555BHVrFlTkrRjxw6NHTtWNpvNtU/37t3Vs2dPz18EAAAAUEYVJpBHRESoZ8+e2rZtm06fPl1qu5ycHLVr104jR46Uv7+/UlNTNWnSJL322muuNqGhoVq0aJE3ygYAAACuSIUJ5K1atZIk7d2796KBvHnz5iWWu3btqlWrVik7O1shISEerREAAABwtwoTyMvr+++/V3h4eIkwfvr0aT300EPy8/NTs2bN9NBDD6latWoGqwQAAAAu7KoO5EeOHFFKSooee+wx17rY2FjNmjVLsbGxOn78uObOnauZM2dq3LhxJfbNyMhQRkaGa9lms7nmocPzrFZriY+Ap9DX4C30NQDlddUG8mPHjmncuHHq1q2b7rzzTtf68PBwhYeHS5Jq1Kihxx57TIMGDdLZs2dL3OiZkpKi5ORk1/KYMWM0efJk710AJIlpRvAa+hq8hb4GoKyuykCemZmpsWPHqn379uratetF2/r4+MjpdMrpdJZYP3DgQHXu3Nm1bLPZlJWV5YlycQFWq1UhISHKzs6Ww+EwXQ4qMfoavIW+Zsa5QTjgalZhArnD4ZDD4VBxcbGKi4tVUFAgHx8f+fqWLPH48eN67rnnFB8fr+7du593nO3btys6OlpRUVE6efKk3njjDd18880KCAgo0S4mJkYxMTGu5czMTH6AGnDu+w54Gn0N3kJfA1BWFSaQL1++XMuWLXMtb968WW3atNGwYcPUs2dPTZgwQU2aNNG//vUvZWRkaPXq1Vq9erWr/Zw5c1SjRg3t3btXM2fOVHZ2toKDg9WsWTP179/fxCUBAAAAl2Rx/n4uRxWVmZlpuoQqxWq1Kjw8XFlZWYwkwaPoa/AW+poZkZGRpksArpiP6QIAAACAqoxADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRwAAAAwiEAOAAAAGEQgBwAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAb5mi6govD395fNZjNdRpVhsVgkSXa7XU6n03A1qMzoa/AW+hqA8iKQ/38FBQUqKCgwXUaVYbVa5e/vr5ycHDkcDtPloBKjr8Fb6GtmMJiGyoApKwAAABXczTffrMTExDLts3//flksFq1cudIzRcFtCOQAAACAQQRyAAAAwCACOQAAQBkkJiaqadOm+uCDD3TjjTcqMDBQd911l/bv368TJ06oZ8+eCgkJ0XXXXafly5eX2DclJUUNGzaUzWZT3bp1NWnSJBUXF5do8/nnn6t58+YKCAhQ06ZNtX79+gvW8cUXX6hNmzay2+0KDQ1V79699csvv3jsuuE5BHIAAIAyOnLkiJ5++mk999xzSk1N1Z49e9SnTx/16tVLN9xwg9599101b95cffv21YEDByRJ//u//6tBgwapffv2WrNmjRITE5WUlKSRI0eWOG779u1ls9m0YsUKjRgxQo8//rgOHz5c4vxffPGF4uPjFRoaquXLl+uNN97Q119/rS5dunj16wD34CkrAAAAZXTixAl9/PHHatKkiSTp559/1uDBgzVq1CiNGzdOknTrrbdq1apV+vvf/64nn3xSEydOVEJCgl599VVJ0j333KOCggK9/PLLGj16tKpXr66ZM2fKYrFo/fr1Cg0NlSTVrl1bbdu2LXH+Z599Vi1atNCqVatcj9y84YYb1LRpU61bt04dO3b01pcCbsAIOQAAQBnVrFnTFcYlqUGDBpKkdu3audaFhYUpKipKBw8e1K5du5SZmakePXqUOE6vXr1UUFCgLVu2SJK++uortW7d2hXGJalNmzaKiIhwLefm5mrz5s3q0aOHHA6HioqKVFRUpAYNGqh27dr6+uuvPXLN8BwCOQAAQBmFhYWVWPb39y91fX5+vrKysiRJ0dHRJbafWz5x4oQkKSMjQ1FRUeed77frsrKy5HA4NHz4cPn5+ZX4l56eroMHD17RtcH7mLICAADgYedGuH9/0+XRo0dLbI+JibngjZm/XRcWFiaLxaIxY8aoa9eu57WNjIx0V9nwEgI5AACAhzVs2FA1atTQO++8o/vvv9+1fsWKFfL391fLli0lSS1bttTcuXN16tQp17SVjz76yDWCLkl2u1133HGHdu7cqUmTJnn3QuARBHIAAAAPs1qtGjdunIYMGaKoqCh17NhRX375paZOnaphw4apevXqkqRhw4Zpzpw56tChg5599lllZWVpwoQJru3nTJ8+XW3atFGvXr2UkJCg8PBwHTp0SBs3btSAAQMUHx9v4CpRXswhBwAA8ILBgwdr7ty5Wrdune677z7NmzdPSUlJmjZtmqtNTEyM1q9fr7y8PPXo0UNTp07VnDlzFBsbW+JYrVq10meffaYzZ85owIAB6tixoyZOnKigoCDVr1/f25eGK2RxOp1O00VUBJmZmaZLqFKsVqvCw8NdN6YAnkJfg7fQ18xgvjQqA0bIAQAAAIOuKJAXFhbq9ddf18MPP6x77rlHu3fvliQtX75cO3fudEuBAAAAQGVW7ps69+7dq3bt2ikzM1O33HKLPvvsM50+fVqS9Mknn2jDhg2aP3++2woFAAAAKqNyB/IhQ4aoRo0a2rJli8LCwlwPxJeku+66S6NHj3ZLgQAAABXBuYFHd6tWrZpHjourR7kD+aZNm/T2228rMjLyvJtXrrnmGmVkZFxxcQAAAEBlV+455L6+virtAS1Hjx5VcHBwuYsCAAC4Wp05c0YLFizQmTNnTJeCq0S5A/ldd92ll19+WYWFha51FotFTqdTb7zxhtq2beuWAgEAAK4mOTk5WrhwoXJyckyXgqtEuaesTJ06Va1atVLjxo3VuXNnWSwWzZkzR99//712796tLVu2uLNOAAAAoFIqdyBv1KiRvvnmGyUlJentt9+W1WrV2rVr1a5dO6Wmpuq6665zZ50AAAAVRuvWrS/ZJiEh4ZJt0tLS3FEOrnLlDuSSVK9ePS1cuNBdtQAAAFw1ljW9sreoT/j+p4tuj4+PV0JCggYNGiRJ2rZtm9q3b69Ro0bp/fff12effab//ve/atCggSRp165d+sMf/uC6xy8pKUnJycn629/+pocffth13GuuuUbLli1TfHz8FdUP9+GdOgEAACq4zz//XG3bttWkSZP01FNPSZJCQ0M1bty4i+4XERGhiRMn6uzZs94oE+VU7hHyNm3alLrNx8dHoaGhuuWWWzRgwADVqlWrvKcBAACo0j744AN1795dc+bMUZ8+fVzrBw8erJdffln//ve/dcstt1xw33bt2mnv3r2aO3euhg0b5qWKUVblHiEPDQ3VTz/9pE8//VTZ2dkKCAhQdna2Pv30U/3444/KysrSyy+/rMaNG+vbb791Z80AAABVwvr169WtWzctWLCgRBiXfp16MmTIEI0ZM6bU/S0Wi6ZMmaIpU6bwGMYKrNyBvEePHgoLC9NPP/2krVu3at26ddq6dat2796t0NBQ9e/fX3v37lX9+vV5104AAFDpZBYUXtG/y5GWlqY6deqUOjNhxIgR2rJliz755JNSj3H33XeradOmeuWVV8p1nfC8ck9ZSU5O1gsvvKC4uLgS6+vWrasJEybo2WefVf/+/fXMM8+4bkYAAACoLJ788YDHz5GcnKx//OMfuvfee/XPf/5T1apVK7E9LCxMo0aN0ujRozVv3rxSj/PCCy/onnvu0V//+ldPl4xyKHcgT09Pl8ViueA2i8Wiw4cPS5Jq1qypoqKi8p4GAACgQprdIO7SjS7icgJ9YGCg1qxZo06dOqlDhw7asGHDee+GPnjwYM2aNUtr164t9Ti33Xab2rRpoxdeeOGKaoZnlDuQ33rrrRo/frxatGih2rVru9YfOHBAEyZMUMuWLSVJ+/fv56ZOAABQ6UT6+3nlPIGBgVq7dm2JUP777ePHj9dzzz130eNMnjxZLVu2dD0WERVHueeQz507VydOnNB1112n5s2b695771Xz5s1Vv359nTx5Uq+//rok6ejRo3rsscfcVjAAAEBVExQUpLVr18pqtapjx47Kyckpsf3hhx9WeHj4RY/RuHFjde/eXbm5uZ4sFeVgcV7Bn0n5+fl66623tHXrVmVkZCgmJka33nqrBgwYoICAAHfW6XGZmZmmS6hSrFarwsPDlZWVJYfDYbocVGL0NXgLfc2MyMhIr53r9OnTrs9bt27tljcGSktLO29eOKqeck1Zyc/P18iRI9WvXz898cQT7q4JAACgwrvUO20Cl6tcgTwgIEBvvfWWunXr5u56AAAAKry0tLRStx09elQJCQlatmyZoqOjvVgVrlblnkPeqlUrffnll+6sxW3Wrl2rp556Sg888ICmT59uuhwAAACgVOV+ysrEiRPVp08f180F0dHR5z0GMSIi4ooLLI+IiAj17NlT27ZtKzHfCwAAwNPsdrv69+8vu91uuhRcJcp9U6ePz/8Nrpf2PHLTN7UsXbpUhw8f1ogRIy7Zlps6vYubn+At9DV4C33NDFM3dboTN3Wi3CPkb731VqlBHAAAAMDlKXcgT0xMdGMZ3peRkaGMjAzXss1mU82aNQ1WVLVYrdYSHwFPoa/BW+hrlR8j2fCUcgfyq11KSoqSk5Ndy2PGjNHkyZMNVvR/Tp48qfz8fNNleFxVemOCgIAAhYWFmS6jSgsJCTFdAqoI+hqAsrqiQP7JJ5/ojTfe0I8//njBALl9+/YrObxHDRw4UJ07d3Yt22w2ZWVlGazoV6dOnVLPHj1UzNvaVio+FotWvPOOQkNDTZdS5VitVoWEhCg7O5t5vfAo+poZl3p3SndiDjk8pdyB/J///Kc6deqkdu3aaevWrerQoYPy8vK0efNmxcbG6q677nJnnWXicDjkcDhUXFys4uJiFRQUyMfHR76+/3e5MTExiomJcS1nZmZWiB+gubm5KnY6NfHaWorwrbIvYFQqJ4qKNH7vYeXm5io4ONh0OVXWuZ8LgKfR1wCUVbkT34QJEzRs2DBNnTpVfn5+ev7559WsWTMdOHBA7du3V5s2bdxZZ5ksX75cy5Ytcy1v3rxZbdq00bBhw4zVVFZZhUVikLxyOFlUZLqEKi0nJ0crVqxQp06dFBAQYLocAADOU+5AvnPnTk2ePFk+Pj6yWCzKycmRJMXFxSkpKUlJSUnq27ev2woti969e6t3795Gzu0uMw4eNV0CUCnk5OTozTffVHx8PIEcgFecOXNGK1euVPfu3XllFJel3IE8ICBAxcXFslgsiomJ0Z49e3TnnXdK+nUu1MGDB91WJADPOHXqlM6ePWu6DI86d29IRZmW5mk2m417FQDDcnJytHDhQnXo0IFAjstS7kB+00036YcfftDdd9+ttm3bavLkyYqMjJSfn5/Gjh2rG264wZ11VjnDa0crjDnklcLJoqIK+YrHqVOnlPBggpzFVWNu1NChQ02X4BUWH4uWvb2MUA7A7ZKSkrRr164S04LhHuVOfMOGDdO+ffskSVOmTNGf//xn11NLYmNjtXr1avdUWEVdFxigSH8/02XADTILCk2XcEFnz56Vs9ipWhOvlW8Ef/xVBkUninR4/N5K/6oHUBG0bt36km0SEhIu2SYtLa3UbfHx8fryyy/l6+srm82mW2+9VbNmzVLDhg3LVOumTZuUkJCgI0eOlGm/skpMTNTChQv1wQcfqG3btpKk/Px8BQYGat++fapbt+4lj2GxWLRz5041atTIo7VWNOX+LdyxY0fX57Vq1dI333yjn376SXl5eWrUqJH8/f3dUiAAz/KN8JVfJP9fAaCs6i9rekX7/5Tw/SXbzJw5U4MGDdKZM2f06KOPKjExUV988cVln6PITQ8WuNzjRERE6LnnnnMFclweH3cdyGKx6Prrr9eNN95IGAcAAHCj4OBg9e3bVzt27NCRI0fUq1cvRUdHq3bt2kpKSlJxcbEkacGCBbr99ts1atQoRUVF6dFHH1WHDh30yy+/KDg4WMHBwdq5c6eSkpJKjODn5+fLYrFo//79kn4d7R40aJDuv/9+BQcH65133pEkFRQUqG/fvqpWrZqaNGmijz/+uESd/fr106FDh/T3v//9gtfhdDr1yiuvqEGDBoqIiFDHjh116NAhSVKrVq0kSc2bN1dwcLDeeOMNd34JK7QyjZC/8sorl93WYrFo+PDhZS4IvzrBo/IqDb6XAIArlZ2drcWLF6tZs2bq3Lmz7r77bi1cuFAnTpxQp06dVKtWLT366KOSpK1bt6pbt246fPiwioqK1L9//3JNWVmyZInWrFmjVatWKT8/X1OnTtX777+v+fPna8GCBUpNTVWXLl20b98+1xs0BQQEaMKECRo7dmyJN2A8Z/bs2UpNTdW//vUv1apVS8nJyUpISNBnn32mzz//XBaLRd988w1TVi7mmWeeuey2BPLysdls8rFYNH7vYdOlwI18LBbZbDbTZQAA3Kgws8Dj53jqqac0evRoBQYG6rbbbtP06dPVtWtXTZo0SRaLRTVr1tRTTz2l+fPnuwJ5VFSUnnnmGVksFvn5lf9+tPvuu881Vz4wMFDSrw/16NevnySpf//+mjlzpv7xj3+UeNT1gAED9NJLL2nJkiXq2bNniWPOnTtXL7/8sms+eVJSkux2u9LT01WnTp1y13q1K1MgP/dyCDwnNDRUby9bVulvysrKytLQoUM1a9Ysr77tsSk8ig4AKp8DT/7o8XO88sorGjRokGt5xYoVOnbsWInfncXFxapdu7ZruXbt2rJYLFd87ri4uPPW/T40x8XF6fDhkoOIvr6+ev755zVq1Ch17dq1xLb9+/erV69e8vH5v1nTPj4+OnToEIHc04qLi9WuXTulpKTo+uuv98Ypr2pVIbhZrVZJUmRkpKpXr264GgAAyi5udoMr2r88gb5OnTqKjY11zfW+kN+H8QuF8+DgYOXm5rqWLzSd5UL7paenn7fcvXv389r16NFDU6dOPW8eeJ06dfT6668rPj6+1PqrIrfd1HkxTqdTmzZt0unTp71xOgAAAI/zi/S/on/lceutt6pGjRp6/vnnlZOTo+LiYu3evfu8myt/Kzo6WllZWa43SpOkW265RZ9++qn27t2rnJwcJScnX9b5v/vuOy1dulRFRUVavHix9uzZU+LJe+dYLBZNmTJFL774Yon1jz/+uJ577jnt2bNH0q+vmK9YsaJEree2VSVeCeQAAAC4clarVWvWrNHu3bt1/fXXKzw8XD179lRGRkap+zRq1Eh9+/ZV/fr1FRYWpp07d6pt27Z66KGH1Lx5czVt2lTt2rW7rPN37txZ69atU3h4uF544QWtXr1aERERF2zbvn17NW1a8tGQgwcPVkJCgu677z6FhITopptu0j//+U/X9uTkZD3yyCMKCwvTm2++eVk1VQYWp9Pp8bfpczgc8vPz09atW9WsWTNPn65cMjMzTZdQpRw/flx9+vRRamoqU1YM+eWXX/TQQw8pbnYDnkNeSRRmFujAkz9q0aJFioqKMl1OlWO1WhUeHq6srCw5HA7T5VQZkZGRXjvXb1/pb926tVueQ56WlqZq1apdaWm4yjFCDgAAABjE+2UDAACUw+W80yZwOQjkAAAAZZSWllbqtqNHjyohIUHLli1TdHS0F6vC1YopKwAAAIBBXhkht1qt2rdvn2rWrOmN0wEAABhjt9vVv39/2e1206XgKlHmEfLCwkKtXr1ac+bM0UcffXTBNocPH9bEiRNLrIuLi7uit29F5WK32/Xoo4/ywwoAUOkEBwcrMTFRwcHBpkvBVaJMI+THjx/Xn/70J+3cuVPSrw99b9asmRYtWqQ//OEPrnaHDh1ScnKyxo8f795qUWnY7XYNHDiQx4MBAK4aPJ4QnlKmEfJx48YpJydHmzdvVm5urjZs2KCCggLdcccd+uCDDzxVIwAAAFBplSmQb9iwQc8//7zuuOMOBQQE6O6779aWLVvUtWtXderUSampqZ6qEwAAAKiUyjRl5ejRo4qLiyuxzmazacGCBYqNjVX//v2VkZGhO++8061FAgAAmPbbd+p0J6bCoEyBvG7dutq+fbv+9Kc/nbdt0qRJio6O1rBhw9S2bVu3FQgAAABUZmWastK6dWvNmzev1O2DBw/WkiVL9PHHH19xYQAAAFejM2fOaMGCBTpz5ozpUnCVKNMI+ZAhQ/SHP/xBx44d06effqqbb75Z1157bYk2Dz74oGrWrKlNmza5s04AAICrQk5OjhYuXKgOHTrw6ENcljIF8gYNGqhBgwaSpN69e2vDhg3nBXJJuuuuu3TXXXe5p0IAHpW/J09FxwtNlwE3KDpZZLoEAEA5lPudOhs1aqT09HR31gLAgKMzDpouAQDgQVOmTNGPP/6oBQsWXLJtamqq5s2bV+qbP1YW8fHxSkhI0KBBg0yXIukKAvkLL7ygoUOHqnHjxmrRooU7awIAAKjQWrdufck2CQkJl2yTlpZ2wfW/neqSl5cnPz8/+fr+GtvGjBmjMWPGXGalKlPbPn36qE+fPpfdviwOHz6sp556SmlpacrPz1dMTIz+8pe/aNSoUR4539Wk3IF85MiROn78uG677TZVr15d0dHRslgsru0Wi0XfffedW4oE4DnRw2vLN6zcPwpQgRSdLOIVD8CLmtZfdkX7f/9T6YH9tzeE3n777Ro0aJASExPPa1dUVOQK6hVdv3799Ic//EF79uxRUFCQfvjhB/3nP/8xXVaFUO7vYPPmzRkZByqBgOsC5Rfpb7oMuEFhZoHpEgB42P79+1WvXj3Nnz9fycnJCgwM1H//+1899dRTWrlypbKysnT99ddrxowZrvv5kpKStGvXLi1btsy1/6JFizR+/HidPHlSiYmJmjFjhiRpwYIFev311/Xll19K+nWANSUlRa+88ooyMjL05z//WW+99Zb8/X/9vTFjxgy99NJLKioq0qhRozRt2jQtW7ZM8fHx59W+ZcsWTZ061fXc9caNG6tx48au7Ze6hh07dig0NFQrV65UdHS0UlNTtWvXLtc7yT/33HMaPny4q/327dsVFBSk9957T3Xq1NFrr71W6j2OS5Ys0YsvvqhDhw7phhtu0Ouvv64mTZq44Tt2ecodyC9nHhIAAADcb8OGDfruu+/k5+cn6deB0jFjxigsLEyzZ89Wjx49tH//fgUFBV1w/w8++EDff/+9jh49qmbNmunPf/6z2rRpc8G2K1eu1ObNm+VwOHTHHXdoyZIl+stf/qKNGzdqypQp2rhxoxo1aqQRI0YoMzOz1JpbtWqlJ598UoMHD1arVq3OezDIpa5h7dq1euedd/Tmm2/q2WefVY8ePdShQwft2rVL27dv15133qnu3burdu3akqT3339f8+fP14IFC5SamqouXbpo3759Cg8PL3HeNWvWaOzYsVqzZo0aN26sv/3tb/rzn/+sXbt2uf7w8LQyPYcccJe8vDwdOnRIBQWM6AEArk4FhZlX9O9KJCUlKSQkRIGBgZJ+nfsdGRkpX19fDRs2TIWFhdq5c2ep+ycnJ8tut+vaa6/Vn/70J3377belth0zZoyqV6+uqKgoderUydX27bffVv/+/XXzzTcrICBAkydPVnFxcanHWbFihe6++25NmzbN9eS+DRs2uLZf6hruuOMOde7cWVarVQ8++KDS09OVlJSkwMBA3Xbbbapbt26J6dI33XST+vXrJ19fX/Xv31/16tXTP/7xj/Pqmjt3rkaNGqUbbrhBVqtVAwcOlMVicb1K4A1Xx6QjVCp///vf9cYbb6i4uFj+/v56+umneUwmAOCq8+OBJ42dOy4ursTySy+9pHnz5unnn3+WxWJRdnb2RUerr7nmGtfnQUFBF30To9+3PXLkiCTp559/1k033eTaFhISotDQ0FKPExYWpkmTJmnSpEk6deqUpkyZou7duys9PV0RERGXvIbf13Gp66hTp06J88fFxenw4cPn1bV//36NGDFCo0ePdq0rKCi4YFtPIZDDq7766iulpKTI6XRK+rXDT506VbVq1VL9+vUNVwcAwOVrEDf7iva/kkD/2wdpfPrpp3rhhReUlpampk2bysfHR+Hh4a7ftZ5Ss2ZNHTz4fzeSZ2dn69SpU5e1b2hoqMaNG6dp06Zp3759+s9//uP2a/j947nT09PVvXv389rVqVNHI0eOvOBNs97ClBV41ZYtW0r8EJEkX1/fi75UBgBAReTvF3lF/9zl9OnT8vX1VWRkpIqKijR58mRlZ2e77fil6dWrlxYtWqTt27fr7NmzGjdunHx8So+WI0eO1Pbt21VUVKS8vDzNmDFD4eHhatSokUeu4bvvvtPSpUtVVFSkxYsXa8+ePerYseN57R5//HG9+OKL+u677+R0OnXmzBmtWbNGp0+fvqLzlwWBHF51oUczOZ1OWa1WA9UAAHD1a9++vTp16qRGjRopLi5Ofn5+rhsbPX3ekSNH6t5771Xt2rVVp04dhYaGymazXbB9Xl6eunfvrvDwcNWuXVubNm3SunXrZLfbPXINnTt31rp16xQeHq4XXnhBq1evVkRExHntunTpogkTJqh///4KCwvT9ddfr8WLF1/RucvK4vT06xlXiezs7FI7ENxn586deuyxx1w3ffj4+Mjf319LlixRTEyM4eqqliNHjqhbt26Km92Axx5WEoWZBTrw5I969913S8yrhHdYLBb5+/uroKDA41MF8H+8+bv7tyOmrVu3dstzyNPS0lyPAbzaZWVlKSIiQgcPHlRsbKzRWn77qMerAXPI/7+CggKe+OEFsbGxmjhxol599VWdOHFCNWvW1NNPP63g4GCvvjQEKScnx3QJ8JCcnBz+PxlgtVrl7++vnJwcORwO0+VUGQymmbVq1Sp16tRJBQUFGj58uFq2bGk8jF+NCOTwuhYtWig1NVXh4eHKysriFxcA4Kp0sXfarCreeustDRgwQBaLRbfddpuWLl1quqSrEoEcAACgjNLS0krddvToUSUkJGjZsmWKjo72YlXet3btWtMlXFBSUpLpEsqEmzoBAAAAgwjkAAAAbmS329W/f3/Z7XbTpeAqwZQVAAAANwoODjb6JjO4+hDIAQAALkNleTwhKh6mrAAAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgbuoEqriiE0WmS4Cb8L0EgKsTgRyoomw2myw+Fh0ev9d0KXAji49FNpvNdBkAgDIgkANVVGhoqJa9vUxnz541XYpHZWVlaejQoZo1a5bCw8NNl+NxNptNoaGhpssAAJQBgRyowqpCcLNarZKkyMhIVa9e3XA1AACcj5s6AQAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRwAAAAwiEAOAAAAGEQgBwAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgkK/pAs45c+aM5syZo2+//VaBgYG6//771aVLl/Pabdq0Sa+99ppr2el06uzZs3r22WfVqlUr7dixQ2PHjpXNZnO16d69u3r27OmV6wAAAADKosIE8pSUFBUWFmr+/Pn65ZdfNG7cOMXGxqp58+Yl2sXHxys+Pt61/M0332j69Okl2oWGhmrRokXeKh0AAAAotwoxZSU/P1+bN29Wv379FBQUpLp16+qee+7Rxo0bL7nvxo0b9T//8z8lRsQBAACAq0WFCOSHDx+W0+lUXFyca129evWUnp5+0f2ys7O1ZcsWtWvXrsT606dP66GHHtLDDz+sOXPm6PTp0x6pGwAAALhSFWLKSn5+voKCgkqss9vtysvLu+h+H3/8sWJiYtSoUSPXutjYWM2aNUuxsbE6fvy45s6dq5kzZ2rcuHEl9s3IyFBGRoZr2WazqWbNmm64GlwOq9Va4iPgKT4+Pq6P9Dd4Ej/XAJRXhQjkAQEB54Xv3NxcBQYGXnS/Dz74QG3bti2xLjw8XOHh4ZKkGjVq6LHHHtOgQYN09uzZEtNaUlJSlJyc7FoeM2aMJk+efKWXgjIKCQkxXQIqubNnz0qSqlWr5vrZAHjCmTNnlJKSoj59+vCzDUCZVIhAXqtWLUlSenq66tSpI0nat2+f6/ML2bNnj9LT09W6deuLHtvHx0dOp1NOp7PE+oEDB6pz586uZZvNpqysrPJeAsrIarUqJCRE2dnZcjgcpstBJXZuytrp06e51wQelZmZqTfffFNt2rRRZGSk6XKqDP7QRmVQIQJ5QECA/vjHP2rx4sUaPny4jh07pn/9618aOnRoqft8+OGHat68+Xn/Ebdv367o6GhFRUXp5MmTeuONN3TzzTcrICCgRLuYmBjFxMS4ljMzMwmGBjgcDr7u8Kji4mLXR/oaPIm+BqC8KkQgl34dsZ49e7YSExMVGBiobt26uR5l2LNnT02YMEFNmjSRJBUWFurjjz/W4MGDzzvO3r17NXPmTGVnZys4OFjNmjVT//79vXotAAAAwOWqMIE8ODhYzz777AW3rVixosSyn5+fUlNTL9i2a9eu6tq1q7vLAwAAADyiQjz2EAAAAKiqKswIOQCg8jp16pTriTeV1bkHA1SVe5JsNptCQ0NNlwFUCgRyAIBHnTp1SgkJD8rpLDZdildc7IEElYnF4qNly94mlANuQCAHAHjU2bNn5XQW69paE+XrG2G6HLhBUdEJ7T08vtK/6gF4C4EcAOAVvr4R8vfj+dwA8Hvc1AkAAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCAeewgA8Iq8/D0qLDpuugy4QVHRSdMlAJUKgRwA4BUHj84wXQIAVEhMWQEAAAAMYoQcAOAVtaOHy9c3zHQZcIOiopO84gG4EYEcAOAVgQHXyd8v0nQZcIOCwkzTJQCVClNWAAAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRxApWa32/Xoo4/KbrebLgUAgAsikAOo1Ox2uwYOHEggBwBUWARyAAAAwCACOQAAAGAQgRwAAAAwiEAOAAAAGEQgBwAAAAwikAMAAAAG+ZouoKLw9/eXzWYzXUaVYbFYJP36SDqn02m4GlRm9DXzcnJyTJcAD7Hb7apWrZrpMoCrHoH8/ysoKFBBQYHpMqoMq9Uqf39/5eTkyOFwmC4HlRh9zTwCeeWVk5Oj06dPG62BwTRUBkxZAQAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEI89hNc5nU5t27ZNOTk5ioiIUIMGDUyXBAAAYAyBHF7lcDg0efJkffHFF/Lz81NhYaG6dOmiQYMGmS4NAADACKaswKvWr1+vr776Sk6nUwUFBXI6nXr//fe1ZcsW06UBAAAYQSCHV+3evVvFxcUl1vn6+uqnn34yVBEAAIBZBHJ4VWhoqKxWa4l1xcXFCg0NNVQRAACAWQRyeFXnzp0VEBDgCuW+vr6KjIxU69atDVcGAABgBjd1wqsiIyM1Z84cLVmyRMeOHVPt2rXVr18/BQUFmS4NAADACAI5vC46OlojR45UeHi4srKy5HA4TJcEAABgDFNWAAAAAIMI5AAAAIBBBHIAAADAIAI5AAAAYBCBHAAAADCIQA4AAAAYRCAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRwAAAAwyNd0AQCAqqGo6ITpEuAmfC8B9yKQAwA8ymazyWLx0d7D402XAjeyWHxks9lMlwFUCgRyAIBHhYaGatmyt3X27FnTpXhUVlaWhg4dqlmzZik8PNx0OR5ns9kUGhpqugygUiCQAwA8rioEN6vVKkmKjIxU9erVDVcD4GrCTZ0AAACAQQRyAAAAwCACOQAAAGBQhZlDfubMGc2ZM0fffvutAgMDdf/996tLly4XbNu5c+f/f9e+RZLUuHFjJSUlubZv3rxZCxcu1IkTJ9SoUSMNGTJEUVFR3rgMAAAAoEwqTCBPSUlRYWGh5s+fr19++UXjxo1TbGysmjdvfsH2M2bMUGxs7HnrDx48qFmzZmn06NFq3LixFi9erGnTpumll17y9CUAAAAAZVYhpqzk5+dr8+bN6tevn4KCglS3bl3dc8892rhxY5mPtWnTJjVr1ky33HKLbDabevfurX379ik9Pd0DlQMAAABXpkKMkB8+fFhOp1NxcXGudfXq1dMXX3xR6j5jx46Vw+HQ9ddfr8TERNWpU0eSdODAAV1//fWudkFBQbrmmmt04MABVxsAAACgoqgQgTw/P19BQUEl1tntduXl5V2w/ZQpU9SwYUMVFhZq1apVGj9+vF577TUFBQUpPz9fdrv9ksfKyMhQRkaGa9lms6lmzZpuuiJcyrnn9Z77CHgKfQ3e4uPj4/pIfwNQFhUikAcEBJwXmHNzcxUYGHjB9k2bNpUk+fn5qW/fvkpLS9POnTvVvHlzBQQEKDc395LHSklJUXJysmt5zJgxmjx5sjsuB2UQEhJiugRUEfQ1eNq5dyKtVq1alXinTgDuUyECea1atSRJ6enprmkl+/btu+wpJueetiJJcXFx2rt3r2s5Ly9PR44cKTEdRpIGDhyozp07u5ZtNpuysrLKfQ0oG6vVqpCQEGVnZ8vhcJguB5UYfQ3ecvr0addHm81muJqqgz9+UBlUiEAeEBCgP/7xj1q8eLGGDx+uY8eO6V//+peGDh16Xtv09HQVFhaqbt26Kioq0rvvvquCggI1bNhQkhQfH6+nn35a27ZtU+PGjbV06VLVrVv3vHAfExOjmJgY13JmZia/rA1wOBx83eEV9DV4WnFxsesjfQ1AWVSIQC79OmI9e/ZsJSYmKjAwUN26dXM98rBnz56aMGGCmjRpopMnT2ru3LnKzMyUv7+/6tevr+TkZAUHB0uSateurSFDhmjOnDnKyspSw4YNNXLkSJOXBgAAAJTK4nQ6naaLqAgyMzNNl1ClWK1WhYeHKysri5EkeBR9Dd5y/Phx9enTR6mpqapevbrpcqqMyMhI0yUAV6xCPIccAAAAqKoI5AAAAIBBBHIAAADAIAI5AABuYLfb9eijj5735nQAcCkEcgAA3MBut2vgwIEEcgBlRiAHAAAADCKQAwAAAAYRyAEAAACDCOQAAACAQQRyAAAAwCACOQAAAGAQgRwAAAAwiEAOAAAAGEQgBwAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcAAAAMMjidDqdpotA1ZORkaGUlBQNHDhQMTExpstBJUZfg7fQ1wCUFyPkMCIjI0PJycnKyMgwXQoqOfoavIW+BqC8COQAAACAQQRyAAAAwCACOYyIiYnRhAkTmGcJj6OvwVvoawDKi5s6AQAAAIMYIQcAAAAMIpADqHKWLl2q6dOnmy4DHrRixQrNnDnzstpu2rRJzz33nGcLqgDGjBmj9evXmy4DwAX4mi4AlceYMWN05513qkOHDpKkvXv3asKECerWrZu2bNmi//73v5ozZ45q1aolSTp06JCeeOIJvf/++5J+DUnLli3Tk08+qXvuucd13IceekgjRozQDTfc4P2LgluMGTNGP/zwg6xWq/z8/FS/fn09+uijio2NLdNxduzYoenTp2vRokUeqvRXM2fO1EcffaTnn39eN910kySpoKBA3bt315tvvqno6OhLHqNz58567bXXynyNVVnPnj1dnxcUFMhqtcpqtUqSunfvXmJ7WY51KfHx8YqPj7/s9mVx/PhxzZs3Tzt27FBBQYEiIiLUrl07devWzSPnA3B1IpDDI3bu3KlJkybpoYceUvv27bVlyxbZ7XalpqZq5MiRpe5XrVo1LV++XK1bt5afn58XK4anPfLII+rQoYPy8vI0e/ZszZo1q0yj1A6Hwy11XO5xqlWrpsWLF7sCOTxvxYoVrs+feeYZdejQQW3btj2vncPhcAX1iu6VV15R7dq1lZKSIpvNpsOHDys9Pd10WQAqGAI53G7btm2aOnWqBg4cWGLUqVOnTnrvvfe0Z88eXXfddRfc96abbtKRI0e0bt06denSxUsVw5sCAwMVHx+v6dOnKysrS2+++aZ27NghPz8/tWvXTgkJCfLx8dGHH36o9evXq2nTpvrwww/VokULffrppyosLHSNfr788sv69NNPdfjwYY0YMULS+SPZM2fOlL+/v06ePKlt27bpr3/9qySpsLBQL7/8srZs2aLIyEg9/vjjatq0qavO+Ph4ff755/ryyy91++23n3cdTqdT7733njZs2KDs7Gw1bNhQf/3rXxUZGen6o3P48OGyWCz6y1/+onvvvdfTX9pK6+jRo3r00Uc1dOhQvf3227LZbJozZ47mzZunzZs3KycnRzExMXrkkUdc38OlS5e6+sW5/YcPH67U1FTl5OSobdu2euSRRyTJ1ddeeuklSb++uvHEE0/ovffeU1ZWlm699VYNGTLENUjw3nvvafXq1XI4HOrWrZtWrVpV6qt4u3fvVmJiooKCgiRJderUUZ06dVzbL3UNBw4cUFBQkD7//HOFhYXp6aef1qFDh5Samqr8/Hz17NnT9bNy6dKl2r9/v2w2W6n9+rfS0tL07rvv6vjx44qLi9MTTzxRojYA3sMccrjVN998oxdffFFDhw497yXg8PBw3XfffVq8eHGp+1ssFvXr108rV65UXl6eh6uFCbm5udq0aZOuvfZaTZo0STExMZo3b56mT5+ur776Shs3bnS1/emnnxQSEqL58+dr0KBBmjBhgkJDQ7VixQqtWLFCtWvXvqxzbtq0Sffdd5+WL1/uCtdbtmxRs2bNtHTpUj3wwAOaPHmyzpw549rH399fCQkJWrJkiYqLi8875j/+8Q99/PHHSk5O1qJFi3Tddde5RvynTZsmSZoxY4ZWrFhBGHeTb7/9Vq+++qpeeeUVSdJ1112nmTNnaunSpWrdurWmTp2qs2fPlrr/tm3bNHv2bM2YMUMffvihvvvuu1Lbfv7555o6dapef/11/fDDD9q0aZMk6d///rfeeecdjR8/XvPmzdPRo0eVnZ1d6nEaNWqklJQUpaWl6ciRI+dtv9Q1fP3117r99tu1dOlS3X777Zo6dap27dql1157TePHj9fChQt17NgxV/tL9evftktNTdUzzzyjJUuWKD4+Xs8//7wKCwtLvRYAnkMgh1vt2LFDNWrU0I033njB7ffff792796t77//vtRj3HLLLapTp47ee+89T5UJA+bNm6cHH3xQTzzxhAoLCzVgwABlZmaqb9++8vf3V/Xq1dW1a1d98sknrn3CwsJ0//33y9fXVzabrdznvvXWW3XjjTfKYrG4jlOvXj21bt1aVqtVbdu2VVRUlL7++usS+7Vr104Oh8MVxn5r/fr16tu3r6Kjo+Xr66sHH3xQu3fvLhGO4F4PPviggoKCXN/D+Ph4hYSEyGq1qkuXLioqKtLBgwdL3b93794KCAjQNddcoyZNmmjv3r2ltu3Ro4dCQkIUFhamFi1aaM+ePZKkTz75RG3atNG1114rf39/9evXTxd7evDIkSN18803a9WqVRo0aJAGDRqkb775xrX9UtfQqFEj3XbbbbJarfrTn/6kY8eO6cEHH5TNZlPDhg0VFRWl/fv3u9pfTr+Wfu2/DzzwgOrWrSur1ap7771XFotFP/zwQ6nXAsBzmLICt+rdu7e2bt2qpKQkJSUluV6mPSc4OFgPPPCAFi1apCFDhpR6nIceekgTJkxQx44dPV0yvOThhx923fArSZ999plOnTql3r17u9YVFxcrMjLStRwZGSmLxXLF565Ro8Yl10VFRenEiRMl1lmtVvXp00cLFiw4b9rK0aNHNW3aNPn4/N+4ho+PjzIzMy94Ply5339dV69erY0bN+rEiROyWCzKzc296Gh1eHi463ObzXbRV+F+3zYrK0uSdOLECdWrV8+1LSgo6Lyfc78VHBysvn37qm/fvsrJydE777yjqVOnat68eapWrdolr+H3dVzqOi6nX0vSL7/8ogULFpR4xbKwsFDHjx8v9VoAeA6BHG7l7++vsWPHauLEiUpOTlZSUpICAwNLtLnvvvu0Zs2aC47anNOwYUPdeOONWrlypadLhiGRkZGKjIzU3/72t8ve50LhPDAwsMRL/OeC06X2+/1I9rFjx9SqVavz2v3xj3/Uu+++q3/+858l1teoUUNPPPEET//xot9+H//zn/9o5cqVmjRpkuLi4uTj46MHH3zwoqPV7hAREaHMzEzXcm5urnJzcy9rX7vdrl69emnVqlU6evSo0tPT3X4Nl9uvIyMj9cADD1zwplkA3seUFbidzWbT+PHj5evrq+TkZOXn55+3PSEh4ZJhu2/fvtqwYQNzySup66+/XiEhIVq2bJny8/NVXFysn3/++aLTmcLCwnTmzJkSc2KvvfZa/ec//9GRI0eUn5+vZcuWXdb59+3bp48//lgOh8M1v7dFixbntfvtfQ2/1aFDBy1evFgZGRmSpDNnzuizzz4rUeu5bXC/vLw8+fj4KCQkRA6HQytWrPDKz4o777xTaWlp2rdvnwoLC5WamnrRV3Hmz5+vffv2yeFw6OzZs3rvvfcUHBys2NhYj1zD5fbrDh06aOXKldq3b5+cTqfy8vK0ZcuWy/7jAoB7MUIOj7DZbBo3bpxrpPz3N1rdfffdWr16tU6fPl3qMerUqaNWrVrpo48+8nS5MMBqtWrcuHFasGCBBg0apPz8fEVHR1/0+cyxsbGKj4/XwIEDVVxcrGnTpummm25SmzZtNHz4cNntdvXt21cffvjhJc/fsmVLbd26Va+99poiIyM1evRoVatW7YJtmzVrpri4uBJ/LNx3333y8fHR888/r+PHj8tut+vmm2/W//zP/0j6dfrW7NmzVVBQoMTERLVv376MXyFczC233KIWLVroiSeeUEBAgDp37lxiupOnNGvWTA888ICSk5PlcDj0wAMPyG63l/qY1oKCAk2dOlVZWVny9fVVvXr1NH78eAUEBHjkGi63X99+++06e/asZs6cqaNHj8pms6lx48alPpEFgGdZnJ5+fQ8AgErqzJkz6t27t9566y2v/EFwMb991COAqwtTVgAAKIPPP/9chYWFys3N1d/+9jc1aNDAeBgHcHVjygoAAGXwwQcf6NVXX5XFYlGDBg309NNPmy4JwFWOKSsAAACAQUxZAQAAAAwikAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcQKV28803KzExsUz77N+/XxaLRStXrvRMUQAA/AaBHAAAADCIQA4AAAAYRCAH4DWJiYlq2rSpPvjgA914440KDAzUXXfdpf379+vEiRPq2bOnQkJCdN1112n58uUl9k1JSVHDhg1ls9lUt25dTZo0ScXFxSXafP7552revLkCAgLUtGlTrV+//oJ1fPHFF2rTpo3sdrtCQ0PVu3dv/fLLLx67bgAALoZADsCrjhw5oqefflrPPfecUlNTtWfPHvXp00e9evXSDTfcoHfffVfNmzdX3759deDAAUnS//7v/2rQoEFq37691qxZo8TERCUlJWnkyJEljtu+fXvZbDatWLFCI0aM0OOPP67Dhw+XOP8XX3yh+Ph4hYaGavny5XrjjTf09ddfq0uXLl79OgAAcI6v6QIAVC0nTpzQxx9/rCZNmkiSfv75Zw0ePFijRo3SuHHjJEm33nqrVq1apb///e968sknNXHiRCUkJOjVV1+VJN1zzz0qKCjQyy+/rNGjR6t69eqaOXOmLBaL1q9fr9DQUElS7dq11bZt2xLnf/bZZ9WiRQutWrVKFotFknTDDTeoadOmWrdunTp27OitLwUAAJIYIQfgZTVr1nSFcUlq0KCBJKldu3audWFhYYqKitLBgwe1a9cuZWZmqkePHiWO06tXLxUUFGjLli2SpK+++kqtW7d2hXFJatOmjSIiIlzLubm52rx5s3r06CGHw6GioiIVFRWpQYMGql27tr7++muPXDMAABdDIAfgVWFhYSWW/f39S12fn5+vrKwsSVJ0dHSJ7eeWT5w4IUnKyMhQVFTUeef77bqsrCw5HA4NHz5cfn5+Jf6lp6fr4MGDV3RtAACUB1NWAFRo50a4f3/T5dGjR0tsj4mJueCNmb9dFxYWJovFojFjxqhr167ntY2MjHRX2QAAXDYCOYAKrWHDhqpRo4beeecd3X///a71K1askL+/v1q2bClJatmypebOnatTp065pq189NFHrhF0SbLb7brjjju0c+dOTZo0ybsXAgBAKQjkACo0q9WqcePGaciQIYqKilLHjh315ZdfaurUqRo2bJiqV68uSRo2bJjmzJmjDh066Nlnn1VWVpYmTJjg2n7O9OnT1aZNG/Xq1UsJCQkKDw/XoUOHtHHjRg0YMEDx8fEGrhIAUJUxhxxAhTd48GDNnTtX69at03333ad58+YpKSlJ06ZNc7WJiYnR+vXrlZeXpx49emjq1KmaM2eOYmNjSxyrVatW+uyzz3TmzBkNGDBAHTt21MSJExUUFKT69et7+9IAAJDF6XQ6TRcBAAAAVFWMkAMAAAAGEcgBAAAAgwjkAAAAgEEEcgAAAMAgAjkAAABgEIEcAAAAMIhADgAAABhEIAcAAAAMIpADAAAABhHIAQAAAIMI5AAAAIBB/w+syawByayd6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ggplot: (1442454920029)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boxplot_metrics({\"PerturbNet\":results_PerturbNet, \"Training Sample\":results_sample,\n",
    "                \"KNN\": results_KNN }, metric_key = \"r2_large\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Perturbnet_scvi0.7",
   "language": "python",
   "name": "perturbnet_scvi0.7"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}