Solvents_predictor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solvent selection using Chem_VAE-generated molecular descriptors - with UDM for filehandling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Workflow:\n",
    "1. Create 2 pandas df from my 2 UDM files. One with the CAS and SMILES of all solvents used, and another with the CAS, d.e. and conv. Use these two UDM files to create two pandas dfs.\n",
    "2. Use Chem_VAE to generate an array of molecular discriptors in the latent space for each solvent (1x196)\n",
    "3. 196 discriptors is too many to train a GP model on - do PCA for dimensionality reduction to arrive at fewer discriptors\n",
    "4. Then a random_forest model with the descriptors as x and yields as y is used to predict the optimal next solvent to test\n",
    "5. Finally, a UDM file with the CAS and SMILES of this solvent is produced"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Import the functions needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using standarized functions? True\n",
      "Standarization: estimating mu and std values ...done!\n"
     ]
    }
   ],
   "source": [
    "#Normal imports required\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from numpy import shape\n",
    "import xml.etree.ElementTree  as ET\n",
    "from lxml.builder import ElementMaker # lxml only !\n",
    "from lxml import etree\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "import scikitplot as skplt\n",
    "import datetime\n",
    "from lxml import etree\n",
    "from io import StringIO\n",
    "from urllib.request import urlopen\n",
    "E = ElementMaker()\n",
    "\n",
    "#My functions defined in UDM_converters.py\n",
    "from UDM_converters import udm_to_df\n",
    "from UDM_converters import CAS_to_SMILES\n",
    "from UDM_converters import prediction_udm_file\n",
    "from UDM_converters import UDM_verifier\n",
    "\n",
    "#Imports required for the VAE\n",
    "#from https://github.com/aspuru-guzik-group/chemical_vae\n",
    "\n",
    "# tensorflow backend\n",
    "from os import environ\n",
    "environ['KERAS_BACKEND'] = 'tensorflow'\n",
    "# vae stuff\n",
    "import sys\n",
    "sys.path.insert(1, '/Users/dsw46/OneDrive - University of Cambridge/\\\n",
    "Cambridge/udm_project/UDM_github/UDM/chemical_vae_master')\n",
    "from chemvae.vae_utils import VAEUtils\n",
    "from chemvae import mol_utils as mu\n",
    "# import scientific py\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "# rdkit stuff\n",
    "from rdkit.Chem import AllChem as Chem\n",
    "from rdkit.Chem import PandasTools\n",
    "# plotting stuff\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from IPython.display import SVG, display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "%matplotlib inline\n",
    "\n",
    "#load the vae model\n",
    "vae = VAEUtils(directory='chemical_vae_master/models/zinc_properties')\n",
    "\n",
    "#ML imports\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import explained_variance_score\n",
    "from sklearn.metrics import median_absolute_error\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn import preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Define the functions needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SMILES_to_VAE(dataframe):\n",
    "    \"\"\"\n",
    "    This function takes in a dataframe of molecules (e.g. solvents) that have smiles strings under the heading\n",
    "    \"MOLSTRUCTURE\" and creates an array with the VAE generated molecular descriptors.\n",
    "    \"\"\"\n",
    "    VAE_array = np.zeros([len(dataframe['MOLSTRUCTURE']), 196])\n",
    "    for i in range(len(dataframe['MOLSTRUCTURE'])):\n",
    "        try:\n",
    "            smiles_1=mu.canon_smiles(dataframe['MOLSTRUCTURE'][i])\n",
    "            X_1 = vae.smiles_to_hot(smiles_1,canonize_smiles=True)\n",
    "            z_1 = vae.encode(X_1)\n",
    "            VAE_array[i, :] = z_1\n",
    "        except:\n",
    "            VAE_array[i, :] = 0\n",
    "            print(\"Converting molecule in index \", i, \"didn't work\")\n",
    "    normalized_VAE_array = preprocessing.normalize(VAE_array)\n",
    "    return(normalized_VAE_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 4 unconvertable \"bad SMILES\" in the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_PCA(array, num_components):\n",
    "    '''\n",
    "    PCA takes in an array of size (x,y) with x components, with each component being described by y dimensions.\n",
    "    In this notebook each component is a solvent, and the y-values are generated by the Chem_VAE.\n",
    "    num_components is the number of components you want in your output PCA.\n",
    "    The function returns an array of size (x, z), where z=num_components\n",
    "    Code copied from Kobi's summit\n",
    "    '''\n",
    "    pca = PCA(n_components=num_components)\n",
    "    pcs_practical = pca.fit_transform(array)\n",
    "    print(f\"{round(pca.explained_variance_ratio_.sum()*100)}% of variance\\\n",
    "    is explained by {num_components} principal components\")\n",
    "    return(pcs_practical)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This block will get you a graph showing the variance explained as a funciton of PCA dimensions\n",
    "\n",
    "# pca = PCA(random_state=1)\n",
    "# pca.fit(VAE_list)\n",
    "# skplt.decomposition.plot_pca_component_variance(pca)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above plot shows how much variance is explained depending on the number of principle components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now I can move on to training my surrogate model: A random forest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next step is to link the cas numbers, such that the latent vectors are linked to the correct conversions/d.e.\n",
    "\n",
    "I need to make sure the input dataset (the latent vectors, as 'X') is linked to the corresponding output data (conversion as 'y'). \n",
    "NB: Since solvent choice only has very little influence on the d.e., and the number of solvent experiments is relatively low, we will only use conversion in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def CAS_matching(results_df, all_solvents_df, PCA_all_solvents):\n",
    "    \"\"\"\n",
    "    This function returns an array of Chem VAE mol descriptors reduced by PCA, such that the entries are in the same\n",
    "    order as for the results dataframe (by matching CAS numbers).\n",
    "    \n",
    "    - results_df: The pandas df containing results, i.e. conversion and d.e. Must have 'CAS' as a column.\n",
    "    - all_solvents_df: The pandas df containing all solvents (with CAS numbers). Must have 'CAS' as a column.\n",
    "    - PCA_all_solvents: An array containing the result of applying PCA to the VAE mol descriptors.\n",
    "    \"\"\"\n",
    "    X_temp=np.array([PCA_all_solvents[0]])\n",
    "    for i in range(len(results_df)):\n",
    "        for j in range(len(all_solvents_df)):\n",
    "            if results_df['CAS'][i] == all_solvents_df['CAS'][j]:\n",
    "                a = np.array([PCA_all_solvents[j]])\n",
    "                X_temp = np.concatenate([X_temp,a])\n",
    "    X_temp=np.delete(X_temp,0,0)\n",
    "    X=X_temp\n",
    "    return(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now that we finally have the data in the right shape, let's try implementing random forest regressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def surrogate_model(X, y, PCA_all_solvents, all_solvents_df, model_choice):\n",
    "    \"\"\"\n",
    "    This function implements the random_forest_regressor to predict the optimal solvent for the experiment.\n",
    "    - X: PCA reduced descriptors for each solvent used in the experiment, in same order as y\n",
    "    - y: conversion for each solvent used in the experiment\n",
    "    - PCA_all_solvents: array with the PCA reduced mol descriptors of all solvents\n",
    "    - conv_max: Is a list with 2 numbers (Predicted max conversion achievable, index of this solvent)\n",
    "    - best_solvent: CAS number of the predicted best solvent\n",
    "    \"\"\"\n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state = 0)\n",
    "    \n",
    "    if model_choice == 'random_forest_regressor':\n",
    "        regr = RandomForestRegressor(criterion= 'mae', random_state = 0)\n",
    "    elif model_choice == 'linear_regressor':\n",
    "        regr = LinearRegression()\n",
    "    elif model_choice == 'gaussian_process_regressor':\n",
    "        regr = GaussianProcessRegressor(random_state = 0)\n",
    "    else:\n",
    "        print('That is not a valid surrogate model, please try again')\n",
    "        quit()\n",
    "    \n",
    "    # fit the model\n",
    "    regr.fit(X_train, y_train)\n",
    "\n",
    "    # we can now use it like any other estimator\n",
    "    print('R^2 error on the training data (best is 1.0) is: ', r2_score(y_train, regr.predict(X_train)))\n",
    "    print('R^2 error on the test data (best is 1.0) is: ', r2_score(y_test, regr.predict(X_test)))\n",
    "    \n",
    "    # predicting the best solvent for conversion using the mol descriptors reduced by PCA\n",
    "    predictions = regr.predict(PCA_all_solvents)\n",
    "    conv_max = max_of_array(predictions)\n",
    "    print('Predicted max achievable conversion is ', conv_max[0])\n",
    "    best_solvent = all_solvents_df['CAS'][conv_max[1]]\n",
    "    return(best_solvent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assuming I'd built a model with good predictive power, I'd now like to take the whole solvent dataset into account, and try to predict the best solvent.\n",
    "NB: Since we're just training the model to predict conversion (1D) and not \"conv. and d.e.\" (2D), we don't need a pareto front"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def max_of_array(array):\n",
    "    overall_max = 0\n",
    "    for i in range(len(array)):\n",
    "        temp_max = array[i]\n",
    "        if temp_max > overall_max:\n",
    "            overall_max = temp_max\n",
    "            overall_max_index = i\n",
    "    return(overall_max, overall_max_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now create a UDM file containing just that cas number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Define the solvent predictor function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Solvent_predictor(UDM_schema, solvent_collection, solvent_results, num_PCA_components, output_UDM_file, model_choice):\n",
    "    \"\"\"\n",
    "    This function outputs a UDM file with the predicted best solvent to experiment with next.\n",
    "    \n",
    "    It takes the following inputs:\n",
    "    - solvent_collection: A UDM file containing all solvents to be considered in the experiment. It has colums \n",
    "        with the molecular structure (column name: 'MOLSTRUCTURE', a SMILES string), molecular weight \n",
    "        (name: 'MW'), name (name: 'NAME'), and CAS number (name: 'CAS'). Other properties can also be added in \n",
    "        other columns.\n",
    "    - solvent results: A UDM file containing the solvents used in an experiment which has the same columns as\n",
    "        \"solvent_collection\", and also the columns \"conv\" and \"d.e.\" for conversion and diastereomeric excess,\n",
    "        respectively.\n",
    "    - num_PCA_components: An integer. The number of principle components in the PCA.\n",
    "    - model_choice: Can either be 'random_forest_regressor', 'gaussian_progress_regressor', or 'linear_regressor'\n",
    "    \"\"\"\n",
    "    #Verify the UDM files\n",
    "    UDM_verifier(solvent_collection, UDM_schema)\n",
    "    UDM_verifier(solvent_results, UDM_schema)\n",
    "    \n",
    "    #Import the UDM files\n",
    "    full_all_solvents_df=udm_to_df(solvent_collection)\n",
    "    full_results_df=udm_to_df(solvent_results)\n",
    "\n",
    "    #Extract only the CAS numbers and SMILES strings from the all_solvents_df\n",
    "    #in this instance, the MOLSTRUCTURE is a SMILES string\n",
    "    all_solvents_df = full_all_solvents_df[['CAS', 'MOLSTRUCTURE']]\n",
    "\n",
    "    #Extract only the CAS numbers, conversions and d.e. from the full_results_df\n",
    "    results_df = full_results_df[['CAS', 'conv', 'd.e.']]\n",
    "    \n",
    "    \n",
    "    #Create an array with the VAE mol descriptors for each solvent\n",
    "    VAE_array = SMILES_to_VAE(all_solvents_df)\n",
    "    \n",
    "    #Do PCA on the VAE array\n",
    "    PCA_array = my_PCA(VAE_array, num_PCA_components)\n",
    "    \n",
    "    #Define X and y for the ML model\n",
    "    y = results_df['conv'].copy()\n",
    "    X = CAS_matching(results_df, all_solvents_df, PCA_array)\n",
    "    \n",
    "    # Use the model to predict the best solvent\n",
    "    CAS_best_solvent = surrogate_model(X, y, PCA_array, all_solvents_df, model_choice)\n",
    "    \n",
    "    # Create a UDM file with the result\n",
    "    prediction_udm_file(CAS_best_solvent, output_UDM_file)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XML well formed, syntax ok.\n",
      "XML valid, schema validation ok.\n",
      "XML well formed, syntax ok.\n",
      "XML valid, schema validation ok.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RDKit WARNING: [17:07:49] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:07:49] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:07:49] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:07:49] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:07:49] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:07:49] WARNING: not removing hydrogen atom without neighbors\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERROR: Check chars file. Bad SMILES: CN1CCCCC1.[Cl-].[H+]                                                                                                    \n",
      "Converting molecule in index  345 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: NNc1ccccc1.[Cl-].[H+]                                                                                                   \n",
      "Converting molecule in index  372 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: C1CCNC1.[Cl-].[H+]                                                                                                      \n",
      "Converting molecule in index  393 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: C[Si](C)(C)C                                                                                                            \n",
      "Converting molecule in index  402 didn't work\n",
      "76.0% of variance    is explained by 11 principal components\n",
      "R^2 error on the training data (best is 1.0) is:  0.5622878122505801\n",
      "R^2 error on the test data (best is 1.0) is:  -0.5550496426487648\n",
      "Predicted max achievable conversion is  128.1454150998271\n"
     ]
    }
   ],
   "source": [
    "#Linear regression\n",
    "Solvent_predictor('udm_6_0_0.xsd', 'data/all_solvents.xml', 'data/solvents_results.xml', 11, 'predicted_best_solvent.xml', 'linear_regressor')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XML well formed, syntax ok.\n",
      "XML valid, schema validation ok.\n",
      "XML well formed, syntax ok.\n",
      "XML valid, schema validation ok.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RDKit WARNING: [17:08:32] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:32] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:32] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:32] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:32] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:32] WARNING: not removing hydrogen atom without neighbors\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERROR: Check chars file. Bad SMILES: CN1CCCCC1.[Cl-].[H+]                                                                                                    \n",
      "Converting molecule in index  345 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: NNc1ccccc1.[Cl-].[H+]                                                                                                   \n",
      "Converting molecule in index  372 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: C1CCNC1.[Cl-].[H+]                                                                                                      \n",
      "Converting molecule in index  393 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: C[Si](C)(C)C                                                                                                            \n",
      "Converting molecule in index  402 didn't work\n",
      "76.0% of variance    is explained by 11 principal components\n",
      "R^2 error on the training data (best is 1.0) is:  0.9999999999999983\n",
      "R^2 error on the test data (best is 1.0) is:  -4.692608645446754\n",
      "Predicted max achievable conversion is  170.05527863782982\n"
     ]
    }
   ],
   "source": [
    "#gaussian process\n",
    "Solvent_predictor('udm_6_0_0.xsd', 'data/all_solvents.xml', 'data/solvents_results.xml', 11, 'predicted_best_solvent.xml', 'gaussian_process_regressor')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XML well formed, syntax ok.\n",
      "XML valid, schema validation ok.\n",
      "XML well formed, syntax ok.\n",
      "XML valid, schema validation ok.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "RDKit WARNING: [17:08:43] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:43] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:43] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:43] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:43] WARNING: not removing hydrogen atom without neighbors\n",
      "RDKit WARNING: [17:08:43] WARNING: not removing hydrogen atom without neighbors\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ERROR: Check chars file. Bad SMILES: CN1CCCCC1.[Cl-].[H+]                                                                                                    \n",
      "Converting molecule in index  345 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: NNc1ccccc1.[Cl-].[H+]                                                                                                   \n",
      "Converting molecule in index  372 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: C1CCNC1.[Cl-].[H+]                                                                                                      \n",
      "Converting molecule in index  393 didn't work\n",
      "ERROR: Check chars file. Bad SMILES: C[Si](C)(C)C                                                                                                            \n",
      "Converting molecule in index  402 didn't work\n",
      "76.0% of variance    is explained by 11 principal components\n",
      "R^2 error on the training data (best is 1.0) is:  0.8661882223329469\n",
      "R^2 error on the test data (best is 1.0) is:  -0.2598825857030702\n",
      "Predicted max achievable conversion is  92.67\n"
     ]
    }
   ],
   "source": [
    "#random forest regressor\n",
    "Solvent_predictor('udm_6_0_0.xsd', 'data/all_solvents.xml', 'data/solvents_results.xml', 11, 'predicted_best_solvent.xml', 'random_forest_regressor')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "chemvae",
   "language": "python",
   "name": "chemvae"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}