Notebooks/gpr-sklearn-linearprobes-PCA-git.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0551acf8",
   "metadata": {},
   "source": [
    "#### AI/ML directed high-throughput virtual screening of linear probes\n",
    "\n",
    "#### Siva Dasetty, Maximilian Topel | The Ferguson Lab, Pritzker School of Molecular Engineering, U.Chicago"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68f5ed2d",
   "metadata": {},
   "source": [
    "## imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "031fc96c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# STANDARD IMPORTS -- SEE environment.yml.\n",
    "\n",
    "import os\n",
    "import os.path\n",
    "import sys\n",
    "import importlib\n",
    "import yaml\n",
    "\n",
    "import subprocess\n",
    "\n",
    "from time import time\n",
    "from datetime import datetime\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "\n",
    "import scipy\n",
    "\n",
    "import rdkit\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "from rdkit import DataStructs\n",
    "from rdkit.Chem.Draw import IPythonConsole\n",
    "from rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions\n",
    "from rdkit.Chem import Draw\n",
    "from rdkit.Chem.Draw import rdMolDraw2D\n",
    "from rdkit.Chem import Descriptors\n",
    "from rdkit.Chem.Fragments import *\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "from matplotlib import rcParams\n",
    "import matplotlib as mpl\n",
    "from matplotlib  import cm\n",
    "from PIL import Image\n",
    "from matplotlib import colors\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.offsetbox import (TextArea, DrawingArea, OffsetImage,\n",
    "                                  AnnotationBbox)\n",
    "from matplotlib.cbook import get_sample_data\n",
    "\n",
    "import multiprocessing\n",
    "\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "# Path to the codes directory\n",
    "sys.path.append(\"../Codes/\")\n",
    "\n",
    "# For getting one hot encoding vectors to obtain latent space representation of probes\n",
    "import utils_activeLearning_Oliver as ut \n",
    "\n",
    "# For GPR training\n",
    "import utils_functions as utils_functions\n",
    "\n",
    "#import utils_figures as utils_figures\n",
    "\n",
    "importlib.reload(ut)\n",
    "importlib.reload(utils_functions)\n",
    "#importlib.reload(utils_figures)\n",
    "\n",
    "# For converting SMILES to SELFIES\n",
    "import selfies as sf\n",
    "\n",
    "# For getting latent representation of probe\n",
    "from oliver_vae.VAE import VAE_encode, VAE_decode\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c9af439",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "10a01605",
   "metadata": {},
   "source": [
    "## plot settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b70b55ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 20})\n",
    "\n",
    "# time stamp\n",
    "now = datetime.now()\n",
    "now = str(now)[:10]\n",
    "\n",
    "# Uncomment to save figures in the below directory.\n",
    "#figures_dir = './Figures/{}'.format(str(now)).replace(':','-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1199d0a5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b6ebc79e",
   "metadata": {},
   "source": [
    "## read all search molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "21484b0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read all potential linear probes\n",
    "trainingfolder = '../VAE_data'\n",
    "all_potential_linear_probes = pd.read_csv(trainingfolder + '/SMILES_linearProbes_All_Properties.csv')\n",
    "\n",
    "# extract smiles\n",
    "all_potential_linear_probesList = all_potential_linear_probes.smiles\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8611bd6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7a7611ed",
   "metadata": {},
   "source": [
    "## generate dictionary of smiles and index in search set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fa641c69",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_smiles_all_potential_linear_probes = dict(zip(all_potential_linear_probes.smiles, all_potential_linear_probes.index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ce4dd91",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cdf90e45",
   "metadata": {},
   "source": [
    "## read saved PCA arrays of search set and molecules in ZINC trainings et"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6efb7cac",
   "metadata": {},
   "outputs": [],
   "source": [
    "## read saved PCA arrays\n",
    "projected_all_search_linear_probes_array = np.load('../VAE_data/PCA_AllSearchLinearProbes.npy')\n",
    "\n",
    "## read ZINC training data latent space.\n",
    "projected_zinc_set = np.load('../VAE_data/PCA_ZINC_TrainingSet.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca1cbdb8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "06b61c88",
   "metadata": {},
   "source": [
    "## read completed dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "475a673f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Method</th>\n",
       "      <th>Analyte</th>\n",
       "      <th>lprobe</th>\n",
       "      <th>lflourines</th>\n",
       "      <th>lhydrogen</th>\n",
       "      <th>hhydrogen</th>\n",
       "      <th>hcarbons</th>\n",
       "      <th>Probe</th>\n",
       "      <th>Avg</th>\n",
       "      <th>SD</th>\n",
       "      <th>smiles</th>\n",
       "      <th>formula</th>\n",
       "      <th>selectivity</th>\n",
       "      <th>selectivityError</th>\n",
       "      <th>sensitivity</th>\n",
       "      <th>sensitivityError</th>\n",
       "      <th>smilesValid</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>C4/F0</td>\n",
       "      <td>-0.501825</td>\n",
       "      <td>0.534080</td>\n",
       "      <td>CCCC</td>\n",
       "      <td>C4H10</td>\n",
       "      <td>0.124415</td>\n",
       "      <td>0.541591</td>\n",
       "      <td>-0.501825</td>\n",
       "      <td>0.534080</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>14</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>C6/F0</td>\n",
       "      <td>-1.819519</td>\n",
       "      <td>0.640062</td>\n",
       "      <td>CCCCCC</td>\n",
       "      <td>C6H14</td>\n",
       "      <td>-0.040476</td>\n",
       "      <td>0.679777</td>\n",
       "      <td>-1.819519</td>\n",
       "      <td>0.640062</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>18</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>C8/F0</td>\n",
       "      <td>-2.620010</td>\n",
       "      <td>0.264206</td>\n",
       "      <td>CCCCCCCC</td>\n",
       "      <td>C8H18</td>\n",
       "      <td>1.196455</td>\n",
       "      <td>0.290748</td>\n",
       "      <td>-2.620010</td>\n",
       "      <td>0.264206</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "      <td>22</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>C10/F0</td>\n",
       "      <td>-3.852597</td>\n",
       "      <td>0.212905</td>\n",
       "      <td>CCCCCCCCCC</td>\n",
       "      <td>C10H22</td>\n",
       "      <td>0.595865</td>\n",
       "      <td>0.240479</td>\n",
       "      <td>-3.852597</td>\n",
       "      <td>0.212905</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>12</td>\n",
       "      <td>0</td>\n",
       "      <td>26</td>\n",
       "      <td>-1</td>\n",
       "      <td>-1</td>\n",
       "      <td>C12/F0</td>\n",
       "      <td>-4.290535</td>\n",
       "      <td>0.607122</td>\n",
       "      <td>CCCCCCCCCCCC</td>\n",
       "      <td>C12H26</td>\n",
       "      <td>1.676133</td>\n",
       "      <td>0.666850</td>\n",
       "      <td>-4.290535</td>\n",
       "      <td>0.607122</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>514</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "      <td>20</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>C9Br0/NH2</td>\n",
       "      <td>-4.578463</td>\n",
       "      <td>0.088387</td>\n",
       "      <td>CCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(B...</td>\n",
       "      <td>C9H9Br12N</td>\n",
       "      <td>0.693603</td>\n",
       "      <td>0.232540</td>\n",
       "      <td>-4.578463</td>\n",
       "      <td>0.088387</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>519</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>C3F3/N(CH3)3</td>\n",
       "      <td>-1.359146</td>\n",
       "      <td>0.327960</td>\n",
       "      <td>C[N+](C)(C)CCC(F)(F)F</td>\n",
       "      <td>C6H13F3N</td>\n",
       "      <td>-0.042718</td>\n",
       "      <td>0.408764</td>\n",
       "      <td>-1.359146</td>\n",
       "      <td>0.327960</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>522</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>12</td>\n",
       "      <td>25</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>C12Cl25/P(CH3)2</td>\n",
       "      <td>-7.310449</td>\n",
       "      <td>0.631328</td>\n",
       "      <td>CP(C)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)...</td>\n",
       "      <td>C14H6Cl25P</td>\n",
       "      <td>5.527516</td>\n",
       "      <td>1.097425</td>\n",
       "      <td>-7.310449</td>\n",
       "      <td>0.631328</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>523</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>7</td>\n",
       "      <td>15</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>C7F15/N(CH3)2</td>\n",
       "      <td>-4.964294</td>\n",
       "      <td>0.862647</td>\n",
       "      <td>CN(C)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F...</td>\n",
       "      <td>C9H6F15N</td>\n",
       "      <td>-0.463810</td>\n",
       "      <td>0.935841</td>\n",
       "      <td>-4.964294</td>\n",
       "      <td>0.862647</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>524</th>\n",
       "      <td>PBMD-4CVS-GPU</td>\n",
       "      <td>PFOSD</td>\n",
       "      <td>7</td>\n",
       "      <td>15</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>C7F15/N(CH3)3</td>\n",
       "      <td>-3.640009</td>\n",
       "      <td>0.307166</td>\n",
       "      <td>C[N+](C)(C)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)...</td>\n",
       "      <td>C10H9F15N</td>\n",
       "      <td>-0.278788</td>\n",
       "      <td>0.579390</td>\n",
       "      <td>-3.640009</td>\n",
       "      <td>0.307166</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>252 rows × 17 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            Method Analyte  lprobe  lflourines  lhydrogen  hhydrogen   \n",
       "0    PBMD-4CVS-GPU   PFOSD       4           0         10         -1  \\\n",
       "1    PBMD-4CVS-GPU   PFOSD       6           0         14         -1   \n",
       "2    PBMD-4CVS-GPU   PFOSD       8           0         18         -1   \n",
       "3    PBMD-4CVS-GPU   PFOSD      10           0         22         -1   \n",
       "4    PBMD-4CVS-GPU   PFOSD      12           0         26         -1   \n",
       "..             ...     ...     ...         ...        ...        ...   \n",
       "514  PBMD-4CVS-GPU   PFOSD       9           0         20          2   \n",
       "519  PBMD-4CVS-GPU   PFOSD       3           3         11          0   \n",
       "522  PBMD-4CVS-GPU   PFOSD      12          25          5          0   \n",
       "523  PBMD-4CVS-GPU   PFOSD       7          15          5          0   \n",
       "524  PBMD-4CVS-GPU   PFOSD       7          15          7          0   \n",
       "\n",
       "     hcarbons            Probe       Avg        SD   \n",
       "0          -1            C4/F0 -0.501825  0.534080  \\\n",
       "1          -1            C6/F0 -1.819519  0.640062   \n",
       "2          -1            C8/F0 -2.620010  0.264206   \n",
       "3          -1           C10/F0 -3.852597  0.212905   \n",
       "4          -1           C12/F0 -4.290535  0.607122   \n",
       "..        ...              ...       ...       ...   \n",
       "514         0        C9Br0/NH2 -4.578463  0.088387   \n",
       "519         3     C3F3/N(CH3)3 -1.359146  0.327960   \n",
       "522         2  C12Cl25/P(CH3)2 -7.310449  0.631328   \n",
       "523         2    C7F15/N(CH3)2 -4.964294  0.862647   \n",
       "524         3    C7F15/N(CH3)3 -3.640009  0.307166   \n",
       "\n",
       "                                                smiles     formula   \n",
       "0                                                 CCCC       C4H10  \\\n",
       "1                                               CCCCCC       C6H14   \n",
       "2                                             CCCCCCCC       C8H18   \n",
       "3                                           CCCCCCCCCC      C10H22   \n",
       "4                                         CCCCCCCCCCCC      C12H26   \n",
       "..                                                 ...         ...   \n",
       "514  CCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(B...   C9H9Br12N   \n",
       "519                              C[N+](C)(C)CCC(F)(F)F    C6H13F3N   \n",
       "522  CP(C)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)...  C14H6Cl25P   \n",
       "523  CN(C)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F...    C9H6F15N   \n",
       "524  C[N+](C)(C)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)...   C10H9F15N   \n",
       "\n",
       "     selectivity  selectivityError  sensitivity  sensitivityError  smilesValid  \n",
       "0       0.124415          0.541591    -0.501825          0.534080         True  \n",
       "1      -0.040476          0.679777    -1.819519          0.640062         True  \n",
       "2       1.196455          0.290748    -2.620010          0.264206         True  \n",
       "3       0.595865          0.240479    -3.852597          0.212905         True  \n",
       "4       1.676133          0.666850    -4.290535          0.607122         True  \n",
       "..           ...               ...          ...               ...          ...  \n",
       "514     0.693603          0.232540    -4.578463          0.088387         True  \n",
       "519    -0.042718          0.408764    -1.359146          0.327960         True  \n",
       "522     5.527516          1.097425    -7.310449          0.631328         True  \n",
       "523    -0.463810          0.935841    -4.964294          0.862647         True  \n",
       "524    -0.278788          0.579390    -3.640009          0.307166         True  \n",
       "\n",
       "[252 rows x 17 columns]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Load canonicalized SMILES of all linear probes\n",
    "\n",
    "simdir = '../Data/GPR_training_data/'\n",
    "\n",
    "# READ RUNNING LIST OF ALL PROBES EXPLORED\n",
    "deltaGfile = simdir + '/wsmiles_deltaG-linearProbesBlockArea-b3-pbmd-4cvs.txt'\n",
    "deltaGdata = pd.read_csv(deltaGfile, delim_whitespace=True, names=[\"Method\", \"Analyte\", \"lprobe\", \n",
    "                                                                   \"lflourines\", \"lhydrogen\", \"hhydrogen\", \n",
    "                                                                   \"hcarbons\", \"Probe\", \n",
    "                                                                   \"Avg\", \"SD\", \"smiles\", \"formula\"], comment=\"#\")\n",
    "\n",
    "#print(deltaGdata.Probe.values)\n",
    "\n",
    "sdsfldata = deltaGdata[ (deltaGdata.Analyte == \"SDS\")\n",
    "                      & (deltaGdata.Method == \"PBMD-4CVS-GPU\") ]\n",
    "#sdsfldata = sdsfldata.sort_values(by=['lprobe', 'lflourines'])\n",
    "\n",
    "pfosfldata = deltaGdata[ (deltaGdata.Analyte == \"PFOSD\") \n",
    "                       & (deltaGdata.Method == \"PBMD-4CVS-GPU\") ]\n",
    "\n",
    "pfosfldataCopy = pfosfldata.copy(deep=False)\n",
    "#pfosfldata = pfosfldata.sort_values(by=['lprobe', 'lflourines'])\n",
    "\n",
    "# for ratio measure of selectivity\n",
    "#yflerror = (pfosfldata.Avg.values/sdsfldata.Avg.values)*np.sqrt( (sdsfldata.SD.values/sdsfldata.Avg.values)**2 +  (pfosfldata.SD.values/pfosfldata.Avg.values)**2 )\n",
    "#pfosfldataCopy['selectivity'] = pfosfldata.Avg.values/sdsfldata.Avg.values\n",
    "\n",
    "# for delta delta G measure of selectivity\n",
    "pfosfldataCopy['selectivity'] = pfosfldata.Avg.values - sdsfldata.Avg.values\n",
    "yflerror = np.sqrt( (sdsfldata.SD.values)**2 +  (pfosfldata.SD.values)**2 )\n",
    "pfosfldataCopy['selectivityError'] = yflerror\n",
    "\n",
    "# for sensitivity\n",
    "pfosfldataCopy['sensitivity'] = pfosfldata.Avg.values\n",
    "pfosfldataCopy['sensitivityError'] = pfosfldata.SD.values\n",
    "\n",
    "# filter probes using alphabet list extracted from the training set molecules\n",
    "alphabet_list = pd.read_csv('../Data/VAE_data/alphabet-list.csv')['alphabets'].tolist()\n",
    "#print(alphabet_list)\n",
    "\n",
    "alpha_exists = ut.check_molecules_alphabets_exists(pfosfldataCopy)\n",
    "\n",
    "#print(\"\\n probes alphabet list\")\n",
    "#print(alpha_exists)\n",
    "#print(np.shape(alpha_exists))\n",
    "\n",
    "pfosfldataCopy['smilesValid'] = np.array(alpha_exists)\n",
    "\n",
    "pfosfldataCopy = pfosfldataCopy[ (pfosfldataCopy.smilesValid == True) ]\n",
    "\n",
    "# filter for C < 15\n",
    "pfosfldataCopy = pfosfldataCopy[ pfosfldataCopy['lprobe'] < 15 ] # MAX LENGTH IN DESIGN SPACE: 14\n",
    "\n",
    "# sort by ncarbons; if required.\n",
    "#pfosfldataCopy = pfosfldataCopy.sort_values(by=['lprobe']) #, 'lflourines'])\n",
    "\n",
    "pfosfldataCopy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d427ee3f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c3c31972",
   "metadata": {},
   "source": [
    "## find indices of completed probes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "152d5347",
   "metadata": {},
   "outputs": [],
   "source": [
    "index_completed_probes = []\n",
    "for key in pfosfldataCopy.smiles.values.tolist():\n",
    "    index_completed_probes.append( dict_smiles_all_potential_linear_probes.get(key) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34375414",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f08c46b3",
   "metadata": {},
   "source": [
    "## separate PCs and smiles of completed and unexplored probes using indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ce26bdf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "pc_completed = projected_all_search_linear_probes_array[index_completed_probes]\n",
    "pc_unexplored = np.delete(projected_all_search_linear_probes_array, index_completed_probes, axis=0)\n",
    "\n",
    "smiles_completed = all_potential_linear_probes.iloc[index_completed_probes]\n",
    "smiles_unexplored = all_potential_linear_probes[~all_potential_linear_probes.index.isin(index_completed_probes)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aa99eb7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d822761a",
   "metadata": {},
   "source": [
    "# check smiles extracted and read are the same in the dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0eebfba5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "# two identical dataframes\n",
    "df1 = smiles_completed[['smiles']].reset_index(drop=True)\n",
    "df2 = pfosfldataCopy[['smiles']].reset_index(drop=True)\n",
    "\n",
    "# check if both are equal\n",
    "print(df1.equals(df2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fba5530",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b080bb0b",
   "metadata": {},
   "source": [
    "## set device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d88ed85e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Planning to run on 3 GPUs.\n"
     ]
    }
   ],
   "source": [
    "torch.cuda.empty_cache()\n",
    "n_devices = torch.cuda.device_count()\n",
    "print('Planning to run on {} GPUs.'.format(n_devices-1)) # ONLY IF BOTORCH is used; BELOW IS SIMPLE CPU IMPLEMENTATION."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "203126a3",
   "metadata": {},
   "source": [
    "## prepare data for training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e1bda238",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-9.435065792882614\n",
      "3.100481351149976\n",
      "-0.9333558445990184\n",
      "9.794834454266525\n"
     ]
    }
   ],
   "source": [
    "npcs = 5 #select number of PCs\n",
    "pcListNames=[]\n",
    "for (item1, item2) in zip(['pc']*npcs, np.arange(0,npcs)):\n",
    "    pcListNames.append(item1+str(item2))\n",
    "\n",
    "trainData = pd.DataFrame( pc_completed[:,np.arange(0,npcs)], columns=pcListNames, dtype = 'float' )\n",
    "\n",
    "ysd = pfosfldataCopy.selectivityError.values\n",
    "y = pfosfldataCopy.selectivity.values\n",
    "\n",
    "trainData['f_obj_selectivity'] = y * -1 # for maximizing \n",
    "trainData['f_obj_error_selectivity'] = ysd\n",
    "\n",
    "ysd = pfosfldataCopy.sensitivityError.values\n",
    "y = pfosfldataCopy.sensitivity.values\n",
    "\n",
    "trainData['f_obj_sensitivity'] = y * -1 # for maximizing\n",
    "trainData['f_obj_error_sensitivity'] = ysd\n",
    "\n",
    "print( trainData.f_obj_selectivity.min() )\n",
    "print( trainData.f_obj_selectivity.max() )\n",
    "\n",
    "print( trainData.f_obj_sensitivity.min() )\n",
    "print( trainData.f_obj_sensitivity.max() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c27ef34",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b7cb58b4",
   "metadata": {},
   "source": [
    "## multi-objective Bayesian optimization using random scalarizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "94697a38",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pc0</th>\n",
       "      <th>pc1</th>\n",
       "      <th>pc2</th>\n",
       "      <th>pc3</th>\n",
       "      <th>pc4</th>\n",
       "      <th>f_obj_selectivity</th>\n",
       "      <th>f_obj_error_selectivity</th>\n",
       "      <th>f_obj_sensitivity</th>\n",
       "      <th>f_obj_error_sensitivity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-49.518326</td>\n",
       "      <td>9.005636</td>\n",
       "      <td>-0.347310</td>\n",
       "      <td>-3.799499</td>\n",
       "      <td>12.873461</td>\n",
       "      <td>-0.124415</td>\n",
       "      <td>0.541591</td>\n",
       "      <td>0.501825</td>\n",
       "      <td>0.534080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-49.171539</td>\n",
       "      <td>9.326957</td>\n",
       "      <td>0.078968</td>\n",
       "      <td>-3.641393</td>\n",
       "      <td>12.023334</td>\n",
       "      <td>0.040476</td>\n",
       "      <td>0.679777</td>\n",
       "      <td>1.819519</td>\n",
       "      <td>0.640062</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-49.430145</td>\n",
       "      <td>8.978824</td>\n",
       "      <td>0.882766</td>\n",
       "      <td>-3.796247</td>\n",
       "      <td>11.605146</td>\n",
       "      <td>-1.196455</td>\n",
       "      <td>0.290748</td>\n",
       "      <td>2.620010</td>\n",
       "      <td>0.264206</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-47.139896</td>\n",
       "      <td>7.172888</td>\n",
       "      <td>2.530454</td>\n",
       "      <td>-3.327389</td>\n",
       "      <td>9.268016</td>\n",
       "      <td>-0.595865</td>\n",
       "      <td>0.240479</td>\n",
       "      <td>3.852597</td>\n",
       "      <td>0.212905</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-44.073818</td>\n",
       "      <td>5.382994</td>\n",
       "      <td>4.331920</td>\n",
       "      <td>-1.006031</td>\n",
       "      <td>5.727169</td>\n",
       "      <td>-1.676133</td>\n",
       "      <td>0.666850</td>\n",
       "      <td>4.290535</td>\n",
       "      <td>0.607122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>247</th>\n",
       "      <td>3.287469</td>\n",
       "      <td>-5.507845</td>\n",
       "      <td>-3.342284</td>\n",
       "      <td>0.881792</td>\n",
       "      <td>1.544212</td>\n",
       "      <td>-0.693603</td>\n",
       "      <td>0.232540</td>\n",
       "      <td>4.578463</td>\n",
       "      <td>0.088387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>248</th>\n",
       "      <td>-34.426353</td>\n",
       "      <td>1.413089</td>\n",
       "      <td>1.311949</td>\n",
       "      <td>1.265554</td>\n",
       "      <td>-3.103536</td>\n",
       "      <td>0.042718</td>\n",
       "      <td>0.408764</td>\n",
       "      <td>1.359146</td>\n",
       "      <td>0.327960</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>249</th>\n",
       "      <td>67.572304</td>\n",
       "      <td>-1.427993</td>\n",
       "      <td>0.189692</td>\n",
       "      <td>6.350353</td>\n",
       "      <td>0.993674</td>\n",
       "      <td>-5.527516</td>\n",
       "      <td>1.097425</td>\n",
       "      <td>7.310449</td>\n",
       "      <td>0.631328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>250</th>\n",
       "      <td>20.535923</td>\n",
       "      <td>-3.621423</td>\n",
       "      <td>3.692337</td>\n",
       "      <td>-4.924819</td>\n",
       "      <td>-0.820699</td>\n",
       "      <td>0.463810</td>\n",
       "      <td>0.935841</td>\n",
       "      <td>4.964294</td>\n",
       "      <td>0.862647</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>251</th>\n",
       "      <td>25.844276</td>\n",
       "      <td>-3.530218</td>\n",
       "      <td>4.332822</td>\n",
       "      <td>-1.410431</td>\n",
       "      <td>-1.658447</td>\n",
       "      <td>0.278788</td>\n",
       "      <td>0.579390</td>\n",
       "      <td>3.640009</td>\n",
       "      <td>0.307166</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>252 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           pc0       pc1       pc2       pc3        pc4  f_obj_selectivity  \\\n",
       "0   -49.518326  9.005636 -0.347310 -3.799499  12.873461          -0.124415   \n",
       "1   -49.171539  9.326957  0.078968 -3.641393  12.023334           0.040476   \n",
       "2   -49.430145  8.978824  0.882766 -3.796247  11.605146          -1.196455   \n",
       "3   -47.139896  7.172888  2.530454 -3.327389   9.268016          -0.595865   \n",
       "4   -44.073818  5.382994  4.331920 -1.006031   5.727169          -1.676133   \n",
       "..         ...       ...       ...       ...        ...                ...   \n",
       "247   3.287469 -5.507845 -3.342284  0.881792   1.544212          -0.693603   \n",
       "248 -34.426353  1.413089  1.311949  1.265554  -3.103536           0.042718   \n",
       "249  67.572304 -1.427993  0.189692  6.350353   0.993674          -5.527516   \n",
       "250  20.535923 -3.621423  3.692337 -4.924819  -0.820699           0.463810   \n",
       "251  25.844276 -3.530218  4.332822 -1.410431  -1.658447           0.278788   \n",
       "\n",
       "     f_obj_error_selectivity  f_obj_sensitivity  f_obj_error_sensitivity  \n",
       "0                   0.541591           0.501825                 0.534080  \n",
       "1                   0.679777           1.819519                 0.640062  \n",
       "2                   0.290748           2.620010                 0.264206  \n",
       "3                   0.240479           3.852597                 0.212905  \n",
       "4                   0.666850           4.290535                 0.607122  \n",
       "..                       ...                ...                      ...  \n",
       "247                 0.232540           4.578463                 0.088387  \n",
       "248                 0.408764           1.359146                 0.327960  \n",
       "249                 1.097425           7.310449                 0.631328  \n",
       "250                 0.935841           4.964294                 0.862647  \n",
       "251                 0.579390           3.640009                 0.307166  \n",
       "\n",
       "[252 rows x 9 columns]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multi_trainData = trainData.copy()\n",
    "multi_trainData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "05d9a5ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy score for training data: 0.2072\n",
      "R2 score for training data: 0.9217\n",
      "[2.43601074e+01 2.00000000e+02 9.39990032e-01 6.70665475e-02\n",
      " 1.59722354e+00]\n",
      "Accuracy score for training data: 0.0697\n",
      "R2 score for training data: 0.9832\n",
      "[ 13.71718364 200.           0.34967344   5.48768455   0.59303081]\n",
      "Selected probe #1/20: CNC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(C)(Br)Br\n",
      "iteration 0\n",
      "Accuracy score for training data: 0.2073\n",
      "R2 score for training data: 0.9214\n",
      "[2.43713456e+01 2.00000000e+02 9.39338962e-01 6.70527741e-02\n",
      " 1.59647515e+00]\n",
      "Accuracy score for training data: 0.0694\n",
      "R2 score for training data: 0.9832\n",
      "[ 13.80633268 200.           0.34877534   5.57453071   0.59482758]\n",
      "Selected probe #2/20: CNCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 1\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9194\n",
      "[5.34595110e+00 6.32229673e-02 2.34007591e+00 4.24237973e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0694\n",
      "R2 score for training data: 0.9832\n",
      "[ 13.75148615 200.           0.35084656   5.57149915   0.59387511]\n",
      "Selected probe #3/20: CP(C)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(C)(Cl)Cl\n",
      "iteration 2\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9191\n",
      "[5.34302748e+00 6.31875922e-02 2.33534739e+00 4.24287275e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0692\n",
      "R2 score for training data: 0.9833\n",
      "[ 13.99436788 200.           0.34942656   5.68973963   0.59530911]\n",
      "Selected probe #4/20: CPCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 3\n",
      "Accuracy score for training data: 0.2115\n",
      "R2 score for training data: 0.9187\n",
      "[5.34028771e+00 6.31517278e-02 2.33066568e+00 4.24336065e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0690\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.13148227 200.           0.34947769   5.67431524   0.595291  ]\n",
      "Selected probe #5/20: PCCCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)Br\n",
      "iteration 4\n",
      "Accuracy score for training data: 0.2115\n",
      "R2 score for training data: 0.9184\n",
      "[5.33990294e+00 6.31193931e-02 2.32580149e+00 4.24409590e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0689\n",
      "R2 score for training data: 0.9834\n",
      "[ 14.14479602 200.           0.34877601   5.70285002   0.5987954 ]\n",
      "Selected probe #6/20: CN(C)CCCCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 5\n",
      "Accuracy score for training data: 0.2114\n",
      "R2 score for training data: 0.9181\n",
      "[5.33709595e+00 6.30861505e-02 2.32121535e+00 4.24466537e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0689\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.15832234 200.           0.35103385   5.66786083   0.59760431]\n",
      "Accuracy score for training data: 0.2114\n",
      "R2 score for training data: 0.9181\n",
      "[5.33708884e+00 6.30862223e-02 2.32121494e+00 4.24466872e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0689\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.15832239 200.           0.35103389   5.66785977   0.59760426]\n",
      "Selected probe #7/20: NCCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 7\n",
      "Accuracy score for training data: 0.2114\n",
      "R2 score for training data: 0.9178\n",
      "[5.33488235e+00 6.30515280e-02 2.31934579e+00 4.24553804e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0687\n",
      "R2 score for training data: 0.9834\n",
      "[ 14.15083412 200.           0.35061961   5.66632848   0.60073031]\n",
      "Selected probe #8/20: CCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(N)(Br)Br\n",
      "iteration 8\n",
      "Accuracy score for training data: 0.2106\n",
      "R2 score for training data: 0.9178\n",
      "[2.65143573e+01 2.00000000e+02 9.57167438e-01 6.72206524e-02\n",
      " 1.59961193e+00]\n",
      "Accuracy score for training data: 0.0687\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.14157237 200.           0.34967566   5.93379171   0.60145402]\n",
      "Selected probe #9/20: CCCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(P)(Br)Br\n",
      "iteration 9\n",
      "Accuracy score for training data: 0.2121\n",
      "R2 score for training data: 0.9168\n",
      "[5.39113434e+00 6.31112855e-02 2.31606452e+00 4.28871859e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0685\n",
      "R2 score for training data: 0.9834\n",
      "[ 14.20065171 200.           0.34852395   5.96163376   0.60498121]\n",
      "Selected probe #10/20: C[N+](C)(C)CCCCCCCCCC(F)(F)C(F)(F)F\n",
      "iteration 10\n",
      "Accuracy score for training data: 0.2119\n",
      "R2 score for training data: 0.9172\n",
      "[5.41890052e+00 6.30796414e-02 2.40041062e+00 4.40205290e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0685\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.17570749 200.           0.34835633   5.94631115   0.60456131]\n",
      "Accuracy score for training data: 0.2119\n",
      "R2 score for training data: 0.9172\n",
      "[5.41890151e+00 6.30796446e-02 2.40040923e+00 4.40205471e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0685\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.17570577 200.           0.34835633   5.94630967   0.60456134]\n",
      "Accuracy score for training data: 0.2119\n",
      "R2 score for training data: 0.9172\n",
      "[5.41890299e+00 6.30796338e-02 2.40040972e+00 4.40205173e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0685\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.17570823 200.           0.34835628   5.94631227   0.60456137]\n",
      "Selected probe #11/20: C[N+](C)(C)CCCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)Br\n",
      "iteration 13\n",
      "Accuracy score for training data: 0.2118\n",
      "R2 score for training data: 0.9169\n",
      "[5.41589678e+00 6.30453211e-02 2.39575158e+00 4.40230974e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0685\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.21717923 200.           0.34949425   5.9290396    0.6039387 ]\n",
      "Selected probe #12/20: CC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(N)(Br)Br\n",
      "iteration 14\n",
      "Accuracy score for training data: 0.2118\n",
      "R2 score for training data: 0.9166\n",
      "[5.41287653e+00 6.30140070e-02 2.39117711e+00 4.40262295e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0684\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.33162161 200.           0.34952689   5.9325937    0.60394663]\n",
      "Accuracy score for training data: 0.2118\n",
      "R2 score for training data: 0.9166\n",
      "[5.41287673e+00 6.30139784e-02 2.39117320e+00 4.40263514e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0684\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.33162241 200.           0.34952687   5.9325923    0.60394666]\n",
      "Accuracy score for training data: 0.2118\n",
      "R2 score for training data: 0.9166\n",
      "[5.41287919e+00 6.30140067e-02 2.39117842e+00 4.40262024e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0684\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.33162381 200.           0.34952687   5.93259306   0.60394666]\n",
      "Accuracy score for training data: 0.2118\n",
      "R2 score for training data: 0.9166\n",
      "[5.41287684e+00 6.30140448e-02 2.39117374e+00 4.40262129e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0684\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.33162494 200.           0.34952686   5.93259229   0.60394668]\n",
      "Selected probe #13/20: CCCCCCCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)NC\n",
      "iteration 18\n",
      "Accuracy score for training data: 0.2117\n",
      "R2 score for training data: 0.9163\n",
      "[5.41033318e+00 6.30308924e-02 2.38613885e+00 4.40304517e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.37451315 200.           0.34959463   5.94611495   0.60657789]\n",
      "Accuracy score for training data: 0.2117\n",
      "R2 score for training data: 0.9163\n",
      "[5.41031996e+00 6.30308535e-02 2.38613991e+00 4.40304297e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.37451195 200.           0.34959462   5.94611647   0.60657785]\n",
      "Selected probe #14/20: PCCCCCCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)Br\n",
      "iteration 20\n",
      "Accuracy score for training data: 0.2117\n",
      "R2 score for training data: 0.9160\n",
      "[5.40746772e+00 6.29975171e-02 2.38165161e+00 4.40344210e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.37100147 200.           0.34962288   6.01534901   0.60732884]\n",
      "Accuracy score for training data: 0.2117\n",
      "R2 score for training data: 0.9160\n",
      "[5.40746432e+00 6.29974629e-02 2.38165385e+00 4.40344477e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9833\n",
      "[ 14.37100166 200.           0.34962288   6.01534932   0.60732883]\n",
      "Selected probe #15/20: PCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 22\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9157\n",
      "[5.40464858e+00 6.29643939e-02 2.37721286e+00 4.40387488e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9832\n",
      "[ 14.38217356 200.           0.3489021    6.1474591    0.60813226]\n",
      "Selected probe #16/20: CN(C)CCCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 23\n",
      "Accuracy score for training data: 0.2117\n",
      "R2 score for training data: 0.9153\n",
      "[5.40637654e+00 6.29501713e-02 2.38912495e+00 4.40764195e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9832\n",
      "[ 14.44289049 200.           0.34920777   6.16453905   0.60895242]\n",
      "Selected probe #17/20: CPCCCCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 24\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9150\n",
      "[5.40324581e+00 6.29191979e-02 2.38597949e+00 4.40818645e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9832\n",
      "[ 14.51036363 200.           0.34805598   6.27441387   0.61085695]\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9150\n",
      "[5.40325828e+00 6.29190670e-02 2.38598085e+00 4.40818849e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9832\n",
      "[ 14.5103641  200.           0.34805593   6.27441369   0.61085703]\n",
      "Selected probe #18/20: CCCCCCCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)P(C)C\n",
      "iteration 26\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9147\n",
      "[5.40056959e+00 6.28869706e-02 2.38146250e+00 4.40869618e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9831\n",
      "[ 14.50659261 200.           0.34897159   6.25823816   0.61060125]\n",
      "Selected probe #19/20: CCCC(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(Br)(Br)C(N)(Br)Br\n",
      "iteration 27\n",
      "Accuracy score for training data: 0.2116\n",
      "R2 score for training data: 0.9143\n",
      "[5.40058671e+00 6.28838580e-02 2.37642115e+00 4.41243464e+00\n",
      " 2.00000000e+02]\n",
      "Accuracy score for training data: 0.0683\n",
      "R2 score for training data: 0.9831\n",
      "[ 14.60054735 200.           0.34914714   6.2516548    0.61052074]\n",
      "Selected probe #20/20: PCCC(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)C(Cl)(Cl)Cl\n",
      "iteration 28\n"
     ]
    }
   ],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "npcs=5\n",
    "\n",
    "# random scalarization\n",
    "niter=0\n",
    "q = 20 # number of candidates\n",
    "selected_probes = list()\n",
    "selected_index = list()\n",
    "\n",
    "while (len(selected_probes) < q):\n",
    "    \n",
    "    selectivityData = utils_functions.trainGPR(multi_trainData, 'selectivity', npcs, isotropic=0)\n",
    "    sensitivityData = utils_functions.trainGPR(multi_trainData, 'sensitivity', npcs, isotropic=0)\n",
    "    \n",
    "    gprModel_selectivity =  selectivityData['gprModel']\n",
    "    gprModel_sensitivity =  sensitivityData['gprModel']\n",
    "    \n",
    "    X_selectivity = selectivityData['X']\n",
    "    y_selectivity = selectivityData['y']\n",
    "    ysdplot_selectivity = selectivityData['ysdplot']\n",
    "    y_sensitivity = sensitivityData['y']\n",
    "    ysdplot_sensitivity = sensitivityData['ysdplot']\n",
    "    \n",
    "    theta = utils_functions.get_scalarization()\n",
    "    kappa = utils_functions.loguniform(low=np.log(1e-4), high=np.log(10000))\n",
    "\n",
    "    # search set\n",
    "    x = pc_unexplored[:,0:npcs]\n",
    "    \n",
    "    # EI acq [enable for EI] [NOT USED]\n",
    "    #ei_selectivity = utils_functions.expected_improvement(x, gprModel_selectivity, y_selectivity, greater_is_better=False, n_params=npcs)\n",
    "    #ei_sensitivity = utils_functions.expected_improvement(x, gprModel_sensitivity, y_sensitivity, greater_is_better=False, n_params=npcs)\n",
    "    \n",
    "    # predict values\n",
    "    y_pred_selectivity, y_pred_std_selectivity = gprModel_selectivity.predict(x,  return_std=True)\n",
    "    y_pred_sensitivity, y_pred_std_sensitivity = gprModel_sensitivity.predict(x,  return_std=True)\n",
    "    \n",
    "    # UCB acquisition function.\n",
    "    ucb_selectivity = utils_functions.ucb(y_pred_selectivity.ravel(), y_pred_std_selectivity, kappa)\n",
    "    ucb_sensitivity = utils_functions.ucb(y_pred_sensitivity.ravel(), y_pred_std_selectivity, kappa)\n",
    "        \n",
    "    # random scalarizations\n",
    "    aq = utils_functions.scalarize(ucb_selectivity, ucb_sensitivity, theta).flatten()\n",
    "    \n",
    "    selected_name = smiles_unexplored.smiles.values[aq.argmax()]\n",
    "    \n",
    "    #print(\"iter\")\n",
    "    #print(selected_name)\n",
    "    # update training set\n",
    "    \n",
    "    if selected_name not in selected_probes:\n",
    "        \n",
    "        selected_probes.append(selected_name)\n",
    "        selected_index.append(aq.argmax())\n",
    "        \n",
    "        print(f'Selected probe #{len(selected_probes)}/{q}: {selected_name}')\n",
    "        \n",
    "        #################################################################################################################\n",
    "        # ADD IDENTIFIED PROBE AND ITS PREDICTED RESPONSE TO TRAINING SET FOR NEXT ITERATION [KRIGING BELIEVER APPROACH]        \n",
    "        Xnew = pc_unexplored[aq.argmax()][0:npcs]\n",
    "        \n",
    "        # updated set\n",
    "        Xup = np.vstack((X_selectivity, Xnew))    \n",
    "        yup_selectivity = np.append(y_selectivity, y_pred_selectivity[aq.argmax()]).reshape(-1,1) \n",
    "        yup_sd_selectivity = np.append(ysdplot_selectivity, y_pred_std_selectivity[aq.argmax()]).reshape(-1,1)\n",
    "        yup_sensitivity = np.append(y_sensitivity, y_pred_sensitivity[aq.argmax()]).reshape(-1,1) \n",
    "        yup_sd_sensitivity = np.append(ysdplot_sensitivity, y_pred_std_sensitivity[aq.argmax()]).reshape(-1,1)\n",
    "\n",
    "        pcListNames=[]\n",
    "        for (item1, item2) in zip(['pc']*npcs, np.arange(0,npcs)):\n",
    "            pcListNames.append(item1+str(item2))\n",
    "        \n",
    "        columnNames = pcListNames + ['f_obj_selectivity', 'f_obj_error_selectivity',\n",
    "                                     'f_obj_sensitivity', 'f_obj_error_sensitivity'] \n",
    "        \n",
    "        # update data frame\n",
    "        multi_trainData = pd.DataFrame(np.column_stack((Xup, \n",
    "                                                        yup_selectivity,\n",
    "                                                        yup_sd_selectivity,\n",
    "                                                        yup_sensitivity,\n",
    "                                                        yup_sd_sensitivity)), columns=columnNames)\n",
    "\n",
    "        print(\"iteration\", niter)\n",
    "\n",
    "    if niter > 1000:\n",
    "        break\n",
    "        \n",
    "    niter += 1\n",
    "    \n",
    "selected_probes_df = pd.DataFrame(selected_probes, columns =['smiles']) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5480c3b1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d1488880",
   "metadata": {},
   "source": [
    "### Generate probe names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fec16a05",
   "metadata": {},
   "outputs": [],
   "source": [
    "probeNames = []\n",
    "for i in range(selected_probes_df.shape[0]):\n",
    "#for i in range(c8_12_hyd_fluor_amine_head_probes_df.shape[0]):\n",
    "#for i in range(c3_fluorinated_amine_head_probes_df.shape[0]):\n",
    "\n",
    "    smiles = selected_probes_df.smiles.values[i]\n",
    "    #smiles = c8_12_hyd_fluor_amine_head_probes_df.smiles.values[i]\n",
    "    #smiles = c3_fluorinated_amine_head_probes_df.smiles.values[i]\n",
    "    \n",
    "    label_i_c_count = smiles.count('C')\n",
    "    \n",
    "    label_i_f_count = smiles.count('F')\n",
    "    \n",
    "    label_i_cl_count = smiles.count('Cl')\n",
    "    \n",
    "    label_i_br_count = smiles.count('Br')\n",
    "    \n",
    "    if label_i_cl_count != 0:\n",
    "        label_i_c_count -= label_i_cl_count\n",
    "        \n",
    "    fstring = '-f'\n",
    "    if label_i_cl_count > 0:\n",
    "        label_i_f_count = label_i_cl_count\n",
    "        fstring = '-cl'\n",
    "        \n",
    "    if label_i_br_count > 0:\n",
    "        label_i_f_count = label_i_br_count\n",
    "        fstring = '-br'\n",
    "        \n",
    "        \n",
    "     \n",
    "    countP = -1\n",
    "    if smiles.count('P') > 0:\n",
    "        \n",
    "        countP = smiles.count('P')\n",
    "        smiles = smiles.replace('P', 'N')\n",
    "        \n",
    "    \n",
    "    rdkProbe = Chem.MolFromSmiles(smiles, sanitize=True)\n",
    "    \n",
    "    head_index = None\n",
    "    for atom in rdkProbe.GetAtoms():\n",
    "        if atom.GetAtomicNum() == 7:\n",
    "            head_index = atom.GetIdx()\n",
    "        \n",
    "    #print(head_index)\n",
    "    \n",
    "    label_i_head_group_h_count = -1\n",
    "    label_i_head_group_c_count = -1\n",
    "    \n",
    "    if fr_NH2(rdkProbe) == 1 and fr_quatN(rdkProbe) == 0 :\n",
    "        label_i_head_group_h_count = 2\n",
    "        label_i_head_group_c_count = 0\n",
    "        #print(probes[i], \"primary amine\")\n",
    "    elif fr_NH1(rdkProbe) == 1 and fr_quatN(rdkProbe) == 0:\n",
    "        label_i_head_group_h_count = 1\n",
    "        label_i_head_group_c_count = 1\n",
    "        #print(probes[i], \"secondary amine\")\n",
    "    elif fr_NH0(rdkProbe) == 1 and fr_quatN(rdkProbe) == 0:\n",
    "        label_i_head_group_h_count = 0\n",
    "        label_i_head_group_c_count = 2\n",
    "        #label_i_head_group = \"0\"\n",
    "        #print(probes[i], \"tertiary amine\")\n",
    "    elif fr_NH0(rdkProbe) == 1 and fr_quatN(rdkProbe) == 1:\n",
    "        label_i_head_group_h_count = 0\n",
    "        label_i_head_group_c_count = 3\n",
    "        #label_i_head_group = \"\"\n",
    "        #print(probes[i], \"quaternary amine\")\n",
    "    elif fr_NH2(rdkProbe) == 1 and fr_quatN(rdkProbe) == 1 :\n",
    "        label_i_head_group_h_count = 2\n",
    "        label_i_head_group_c_count = 1\n",
    "        #print(probes[i], \"secondary amine cation\") \n",
    "    elif fr_NH1(rdkProbe) == 1 and fr_quatN(rdkProbe) == 1:\n",
    "        label_i_head_group_h_count = 1\n",
    "        label_i_head_group_c_count = 2\n",
    "        #print(probes[i], \"tertiary amine cation\")\n",
    "    elif fr_NH0(rdkProbe) == 0 and fr_quatN(rdkProbe) == 1:\n",
    "        label_i_head_group_h_count = 3\n",
    "        label_i_head_group_c_count = 0\n",
    "\n",
    "    if label_i_head_group_c_count > 0:\n",
    "        label_i_c_count = label_i_c_count - label_i_head_group_c_count\n",
    "    \n",
    "    if countP == -1:\n",
    "    \n",
    "        if label_i_head_group_h_count == -1 and label_i_head_group_c_count == -1:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count)\n",
    "        elif label_i_head_group_h_count == 0:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count) + '-n' + str(head_index) + '_ch3_' + str(label_i_head_group_c_count)\n",
    "        elif label_i_head_group_c_count == 0:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count) + '-n' + str(head_index) + '_h' + str(label_i_head_group_h_count)\n",
    "        elif label_i_head_group_h_count > 0 and label_i_head_group_c_count > 0:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count) + '-n' + str(head_index) + '_h' + str(label_i_head_group_h_count) + '_ch3_' + str(label_i_head_group_c_count)\n",
    "    else:\n",
    "        \n",
    "        if label_i_head_group_h_count == -1 and label_i_head_group_c_count == -1:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count)\n",
    "        elif label_i_head_group_h_count == 0:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count) + '-p' + str(head_index) + '_ch3_' + str(label_i_head_group_c_count)\n",
    "        elif label_i_head_group_c_count == 0:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count) + '-p' + str(head_index) + '_h' + str(label_i_head_group_h_count)\n",
    "        elif label_i_head_group_h_count > 0 and label_i_head_group_c_count > 0:\n",
    "            label_i = 'c' + str(label_i_c_count) + fstring + str(label_i_f_count) + '-p' + str(head_index) + '_h' + str(label_i_head_group_h_count) + '_ch3_' + str(label_i_head_group_c_count)\n",
    "        \n",
    "    \n",
    "    probeName = label_i\n",
    "    \n",
    "    probeNames.append(probeName)\n",
    "    print(probeName)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe21d92a",
   "metadata": {},
   "source": [
    "### save probes to txt fie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "556b30c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_probes_df = pd.DataFrame(selected_probes, columns =['smiles']) \n",
    "selected_probes_df['name'] = probeNames\n",
    "\n",
    "\n",
    "target_dir = \"./\" #Output directory\n",
    "probe_dir = target_dir + '{}'.format(str(now)).replace(':','-')\n",
    "\n",
    "selected_probes_df.to_csv(probe_dir + \"-mobo-candidates.txt\", \n",
    "                          header=None, index=None, sep='\\t', mode='a')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d861fb83",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e739b172",
   "metadata": {},
   "source": [
    "### visualize the new sample probe, if needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1addc70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mol_with_atom_index( mol ):\n",
    "    atoms = mol.GetNumAtoms()\n",
    "    for idx in range( atoms ):\n",
    "        mol.GetAtomWithIdx( idx ).SetProp( 'molAtomMapNumber', str( mol.GetAtomWithIdx( idx ).GetIdx() ) )\n",
    "    return mol\n",
    "\n",
    "mol_with_atom_index( Chem.MolFromSmiles( selected_probes_df.smiles[0] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f654ab88",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}