biodeg_ml_for_upload.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "506865ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "import chemprop\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import time\n",
    "import umap\n",
    "import xyz2mol\n",
    "from tqdm import tqdm as tq\n",
    "#from tqdm.notebook import tqdm_notebook as tqnb\n",
    "from rdkit import Chem, DataStructs\n",
    "from rdkit.Chem import MolToSmiles, MolFromSmiles, MACCSkeys, AllChem, Descriptors, rdDistGeom, rdMolAlign\n",
    "from rdkit.Chem.AtomPairs import Pairs\n",
    "from rdkit.Chem.SaltRemover import SaltRemover\n",
    "from sklearn import svm\n",
    "from sklearn import tree,manifold\n",
    "from sklearn.naive_bayes import GaussianNB, MultinomialNB\n",
    "from sklearn.linear_model import LinearRegression, LogisticRegression\n",
    "from sklearn.feature_selection import RFE, RFECV, SelectFromModel, VarianceThreshold\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier,RandomForestRegressor\n",
    "from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold, GridSearchCV, RandomizedSearchCV\n",
    "from sklearn.metrics import classification_report, balanced_accuracy_score, confusion_matrix\n",
    "from sklearn.preprocessing import StandardScaler, Binarizer\n",
    "from sklearn.decomposition import PCA, KernelPCA, TruncatedSVD\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.utils.fixes import loguniform\n",
    "from sklearn.utils import shuffle\n",
    "from rdkit.ML.Cluster import Murtagh,Butina,ClusterUtils\n",
    "\n",
    "from collections import OrderedDict\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8724c634",
   "metadata": {},
   "source": [
    "### Processing of CREST data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9cb3e84",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('smiles_part2.csv', low_memory=False, names= [\"SMILES\"])\n",
    "for g in df.SMILES:\n",
    "    try:\n",
    "        m = MolFromSmiles(g)\n",
    "        m2 = Chem.AddHs(m)\n",
    "        cids = AllChem.EmbedMultipleConfs(m2, numConfs=4)\n",
    "        rmslist = []\n",
    "        AllChem.AlignMolConformers(m2, RMSlist=rmslist)\n",
    "        print(g)\n",
    "    except:\n",
    "        None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "853e5123",
   "metadata": {},
   "outputs": [],
   "source": [
    "crest = 'path_to_file'\n",
    "i = 0\n",
    "c_sizes = []\n",
    "filenames = []\n",
    "for fn in os.listdir(crest):\n",
    "    i += 1\n",
    "    if fn[-6:-4] == \"_0\":\n",
    "        c_sizes.append(min(100,i))\n",
    "        i = 0\n",
    "    if i>300: continue\n",
    "    filenames.append(fn)\n",
    "print(len(filenames), len(c_sizes))\n",
    "plt.hist(c_sizes,bins=50,range=(0,100))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "34d6893f",
   "metadata": {},
   "source": [
    "### Write to .sdf\n",
    "extract unique conformers from clustering on conformer IDs (cids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61a2ddf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "with Chem.SDWriter('cr.sdf') as w:\n",
    "    for z, m in enumerate(tq(da['SMILES'])):\n",
    "        try:\n",
    "            mh = Chem.AddHs(MolFromSmiles(m))\n",
    "            ps = rdDistGeom.ETKDGv3()\n",
    "            #ps.numThreads = 0\n",
    "            cids = rdDistGeom.EmbedMultipleConfs(mh,50,ps)\n",
    "            m3d = Chem.RemoveHs(mh)\n",
    "            dists = []\n",
    "            for i in range(len(cids)):\n",
    "                for j in range(i):\n",
    "                    dists.append(rdMolAlign.GetBestRMS(m3d,m3d,i,j))\n",
    "            clusts = Butina.ClusterData(dists, len(cids), 0.5, isDistData=True, reordering=True)\n",
    "            for mol in [x[0] for x in clusts]:\n",
    "                # WRITING FUNCTION!!\n",
    "                w.write(m3d,confId=mol)\n",
    "        except:\n",
    "            None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "433c4b8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "with Chem.SDWriter('ccm.sdf') as w:\n",
    "    for z, m in enumerate(tq(da['SMILES'])):\n",
    "        try:\n",
    "            cids = []\n",
    "            mh = Chem.AddHs(MolFromSmiles(m))\n",
    "            for i in range(50):\n",
    "                AllChem.EmbedMolecule(mh)\n",
    "                AllChem.MMFFOptimizeMolecule(mh)\n",
    "                mol = Chem.RemoveHs(mh)\n",
    "                cids.append(mol)\n",
    "            dists = []\n",
    "            timeout = time.time() + 20\n",
    "            iq = len(cids)\n",
    "            for i in range(len(cids)):\n",
    "                if time.time() > timeout:\n",
    "                    iq = i\n",
    "                    print(len(cids),iq)\n",
    "                    break\n",
    "                for j in range(i):\n",
    "                    dists.append(rdMolAlign.GetBestRMS(cids[i],cids[j]))#, numThreads=0))\n",
    "            clusts = Butina.ClusterData(dists, iq, 0.5, isDistData=True, reordering=True)\n",
    "            for mol in [cids[x[0]] for x in clusts]:\n",
    "                # WRITING FUNCTION!!\n",
    "                w.write(mol)\n",
    "            cids = []\n",
    "        except:\n",
    "            None\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0df3f19",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "m_sizes = []\n",
    "j = 1\n",
    "k = 0\n",
    "mols = []\n",
    "cids = []\n",
    "\n",
    "with Chem.SDWriter('cr.sdf') as w:\n",
    "    for z, fn in enumerate(tq(filenames)):\n",
    "        #if z<17545: continue\n",
    "        if fn[-6:-4] == \"_0\":\n",
    "            timeout = time.time() + 30\n",
    "            dists = []\n",
    "            iq = len(cids)\n",
    "            for i in range(len(cids)):\n",
    "                if time.time() > timeout:\n",
    "                    iq = i\n",
    "                    print(len(cids),iq)\n",
    "                    break\n",
    "                for j in range(i):\n",
    "                    dists.append(rdMolAlign.GetBestRMS(cids[i],cids[j]))#, numThreads=0))\n",
    "            clusts = Butina.ClusterData(dists, iq, 0.5, isDistData=True, reordering=True)\n",
    "            for mol in [cids[x[0]] for x in clusts]:\n",
    "                # WRITING FUNCTION!!\n",
    "                w.write(mol)\n",
    "            cids = []\n",
    "        try:\n",
    "            filename = os.path.join(crest, fn)#, \"crest_conformers.xyz\")\n",
    "            atoms, charge, xyz_coordinates = xyz2mol.read_xyz_file(filename)\n",
    "            cids.append(Chem.RemoveHs(xyz2mol.xyz2mol(atoms, xyz_coordinates, charge=charge)[0]))\n",
    "        except:\n",
    "            None\n",
    "    # k+=1\n",
    "    # if k>=min(50,c_sizes[j]):\n",
    "    #     k = 0\n",
    "    #     j += 1\n",
    "    #     m_sizes.append(sum(atoms))\n",
    "        #m_sizes.append(len([x for x in atoms if x != 1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da0aa184",
   "metadata": {},
   "outputs": [],
   "source": [
    "#reference code for current filter clustering\n",
    "maxd = -100\n",
    "for j in range(0,5):\n",
    "    for i in range(j,len(cids)):\n",
    "        d1 = rdMolAlign.AlignMol(m3d,m3d,prbCid=cids[i],refCid=cids[j])\n",
    "        d2 = rdMolAlign.GetBestRMS(m3d,m3d,prbId=cids[i],refId=cids[j])\n",
    "        delt = d1-d2\n",
    "        if delt<-1e-5:\n",
    "            print(f', {i}, {delt}')\n",
    "        if delt>maxd:\n",
    "            maxd = delt\n",
    "            maxi = i\n",
    "            maxj = j\n",
    "d1 = rdMolAlign.AlignMol(m3d,m3d,prbCid=cids[maxi],refCid=cids[maxj])\n",
    "d2 = rdMolAlign.GetBestRMS(m3d,m3d,prbId=cids[maxi],refId=cids[maxj])\n",
    "print(maxi,maxd,d1,d2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16f18936",
   "metadata": {},
   "source": [
    "\n",
    "### *Par generation. Extract biodeg targets, PubChem, and RDKit methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7b9cbe42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Degradation data (dr). For each scheme, also include other schemes as applicable.\n",
    "da = pd.read_excel('biodeg_target.xls',sheet_name=\"Sheet1\")\n",
    "dr = pd.read_excel('degradation_data.xls')\n",
    "substructs = pd.read_excel('parameters.xls')\n",
    "#da = pd.concat([da,dp],axis=1)\n",
    "#pc_data = da.values[:,3:13]\n",
    "\n",
    "f = []#\"HBondDonorCount\", \"RotatableBondCount\", \"HeavyAtomCount\", \\\n",
    "       #\"ZStericQuadrupole3D\", \"FeatureAnionCount3D\", \"FeatureCationCount3D\", \"FeatureRingCount3D\"]\n",
    "\n",
    "for j in da:\n",
    "    # \"MolecularWeight\", \"ExactMass\", \"MonoisotopicMass\",\"Volume3D\", \"XStericQuadrupole3D\", \"YStericQuadrupole3D\", \"ZStericQuadrupole3D\"\n",
    "    if j in [\"SMILES\",\"Target_INT\", \"Target\",\"Volume3D\",\"MolecularWeight\", \"MonoisotopicMass\", \"UndefinedAtomStereoCount\", \"UndefinedBondStereoCount\"]: continue\n",
    "    f.append(j)\n",
    "f=[]\n",
    "Xd = []\n",
    "MolMethods = ['BalabanJ',\n",
    " 'BertzCT',\n",
    " #'Chi0',\n",
    " #'Chi0n',\n",
    " #'Chi0v',\n",
    " 'Chi1',\n",
    " #'Chi1n',\n",
    " 'Chi1v',\n",
    " 'Chi2n',\n",
    " 'Chi2v',\n",
    " 'Chi3n',\n",
    " 'Chi3v',\n",
    " 'Chi4n',\n",
    " 'Chi4v',\n",
    " #'ChiNn_',\n",
    " #'ChiNv_',\n",
    " 'HallKierAlpha',\n",
    " #'Ipc',\n",
    " 'Kappa1',\n",
    " 'Kappa2',\n",
    " 'Kappa3']\n",
    "_3Descs = [\n",
    " 'CalcAsphericity',\n",
    " #'CalcCoulombMat',\n",
    " #'CalcCrippenDescriptors',\n",
    " #'CalcEEMcharges',\n",
    " 'CalcEccentricity',#\n",
    " 'CalcExactMolWt',\n",
    " 'CalcFractionCSP3',\n",
    " 'CalcInertialShapeFactor',#\n",
    " 'CalcLabuteASA',\n",
    " 'CalcNPR1',\n",
    " 'CalcNPR2',\n",
    " 'CalcPBF',\n",
    " 'CalcPMI1',\n",
    " 'CalcPMI2',\n",
    " 'CalcPMI3',\n",
    " 'CalcPhi',\n",
    " 'CalcRadiusOfGyration',\n",
    " 'CalcSpherocityIndex',\n",
    " 'CalcTPSA',\n",
    "]\n",
    "multi_3Descs=[#'CalcCoulombMat',\n",
    " #'CalcCrippenDescriptors',\n",
    " #'CalcEEMcharges',\n",
    " #'CalcGETAWAY',\n",
    " 'CalcMORSE',\n",
    " 'CalcRDF',\n",
    " 'CalcWHIM',\n",
    "]\n",
    "GETAWAY_entries = [\"ITH\",\"ISH\",\"HIC\",\"HGM\",\"H0u\",\"H1u\",\"H2u\",\"H3u\",\"H4u\",\"H5u\",\"H6u\",\"H7u\",\"H8u\",\"HTu\",\n",
    "\"HATS0u\",\"HATS1u\",\"HATS2u\",\"HATS3u\",\"HATS4u\",\"HATS5u\",\"HATS6u\",\"HATS7u\",\"HATS8u\",\"HATSu\",\"H0m\",\"H1m\",\"H2m\",\"H3m\",\"H4m\",\"H5m\",\n",
    "\"H6m\",\"H7m\",\"H8m\",\"HTm\",\"HATS0m\",\"HATS1m\",\"HATS2m\",\"HATS3m\",\"HATS4m\",\"HATS5m\",\"HATS6m\",\"HATS7m\",\"HATS8m\",\"HATSm\",\"H0v\",\"H1v\",\n",
    "\"H2v\",\"H3v\",\"H4v\",\"H5v\",\"H6v\",\"H7v\",\"H8v\",\"HTv\",\"HATS0v\",\"HATS1v\",\"HATS2v\",\"HATS3v\",\"HATS4v\",\"HATS5v\",\"HATS6v\",\"HATS7v\",\"HATS8v\",\n",
    "\"HATSv\",\"H0e\",\"H1e\",\"H2e\",\"H3e\",\"H4e\",\"H5e\",\"H6e\",\"H7e\",\"H8e\",\"HTe\",\"HATS0e\",\"HATS1e\",\"HATS2e\",\"HATS3e\",\"HATS4e\",\"HATS5e\",\"HATS6e\",\n",
    "\"HATS7e\",\"HATS8e\",\"HATSe\",\"H0p\",\"H1p\",\"H2p\",\"H3p\",\"H4p\",\"H5p\",\"H6p\",\"H7p\",\"H8p\",\"HTp\",\"HATS0p\",\"HATS1p\",\"HATS2p\",\"HATS3p\",\"HATS4p\",\n",
    "\"HATS5p\",\"HATS6p\",\"HATS7p\",\"HATS8p\",\"HATSp\",\"H0i\",\"H1i\",\"H2i\",\"H3i\",\"H4i\",\"H5i\",\"H6i\",\"H7i\",\"H8i\",\"HTi\",\"HATS0i\",\"HATS1i\",\"HATS2i\",\n",
    "\"HATS3i\",\"HATS4i\",\"HATS5i\",\"HATS6i\",\"HATS7i\",\"HATS8i\",\"HATSi\",\"H0s\",\"H1s\",\"H2s\",\"H3s\",\"H4s\",\"H5s\",\"H6s\",\"H7s\",\"H8s\",\"HTs\",\"HATS0s\",\n",
    "\"HATS1s\",\"HATS2s\",\"HATS3s\",\"HATS4s\",\"HATS5s\",\"HATS6s\",\"HATS7s\",\"HATS8s\",\"HATSs\",\"RCON\",\"RARS\",\"REIG\",\"R1u\",\"R2u\",\"R3u\",\"R4u\",\"R5u\",\n",
    "\"R6u\",\"R7u\",\"R8u\",\"RTu\",\"R1u+\",\"R2u+\",\"R3u+\",\"R4u+\",\"R5u+\",\"R6u+\",\"R7u+\",\"R8u+\",\"RTu+\",\"R1m\",\"R2m\",\"R3m\",\"R4m\",\"R5m\",\"R6m\",\"R7m\",\n",
    "\"R8m\",\"RTm\",\"R1m+\",\"R2m+\",\"R3m+\",\"R4m+\",\"R5m+\",\"R6m+\",\"R7m+\",\"R8m+\",\"RTm+\",\"R1v\",\"R2v\",\"R3v\",\"R4v\",\"R5v\",\"R6v\",\"R7v\",\"R8v\",\"RTv\",\n",
    "\"R1v+\",\"R2v+\",\"R3v+\",\"R4v+\",\"R5v+\",\"R6v+\",\"R7v+\",\"R8v+\",\"RTv+\",\"R1e\",\"R2e\",\"R3e\",\"R4e\",\"R5e\",\"R6e\",\"R7e\",\"R8e\",\"RTe\",\"R1e+\",\"R2e+\",\n",
    "\"R3e+\",\"R4e+\",\"R5e+\",\"R6e+\",\"R7e+\",\"R8e+\",\"RTe+\",\"R1p\",\"R2p\",\"R3p\",\"R4p\",\"R5p\",\"R6p\",\"R7p\",\"R8p\",\"RTp\",\"R1p+\",\"R2p+\",\"R3p+\",\"R4p+\",\n",
    "\"R5p+\",\"R6p+\",\"R7p+\",\"R8p+\",\"RTp+\",\"R1i\",\"R2i\",\"R3i\",\"R4i\",\"R5i\",\"R6i\",\"R7i\",\"R8i\",\"RTi\",\"R1i+\",\"R2i+\",\"R3i+\",\"R4i+\",\"R5i+\",\"R6i+\",\n",
    "\"R7i+\",\"R8i+\",\"RTi+\",\"R1s\",\"R2s\",\"R3s\",\"R4s\",\"R5s\",\"R6s\",\"R7s\",\"R8s\",\"RTs\",\"R1s+\",\"R2s+\",\"R3s+\",\"R4s+\",\"R5s+\",\"R6s+\",\"R7s+\",\"R8s+\",\"RTs+\"]\n",
    "#print(len(GETAWAY_entries))\n",
    "prs=[]\n",
    "morserdfsuffix = ['u', 'm', 'v', 'e', 'p', 'i', 's']\n",
    "fnames = {\"MolMethods\": MolMethods, \"_3Descs\": _3Descs, \n",
    "\"CalcGETAWAY\": GETAWAY_entries,\n",
    "\"CalcMORSE\": [\"MORSE_%s\" %(i%32) + morserdfsuffix[i//32] for i in range(224)],\n",
    "\"CalcRDF\": [\"RDF_%s\" %(5*(2+i%30)) + morserdfsuffix[i//30] for i in range(210)],\n",
    "\"CalcWHIM\": [\"WHIM_%s\" %i for i in range(114)],#provisional\n",
    "\"MACCS\":\n",
    "    ['0'] + [i[1][0]+'>'+str(i[1][1]) for i in Chem.MACCSkeys.smartsPatts.items()],\n",
    "\"GCT\": substructs.Pa.tolist(),}\n",
    "def fix_CCS(CCS):\n",
    "    \n",
    "    CCS = MolToSmiles(MolFromSmiles(CCS))\n",
    "    #CCS = CCS.replace('{t}','')\n",
    "    #CCS = CCS.replace('{c}','')\n",
    "    CCS = CCS.replace('Na{+}','Na+')\n",
    "    #CCS = CCS.replace('P{+}','[PH3+]')\n",
    "    #CCS = CCS.replace('N{+}','[NH3+]')\n",
    "    #CCS = CCS.replace('N{-}','[NH-]')        \n",
    "    CCS = CCS.replace('Ca{2+}','Ca++')        \n",
    "    CCS = CCS.replace('Cu{2+}','Cu++')\n",
    "    CCS = CCS.replace('Zn{2+}','Zn++')\n",
    "    CCS = CCS.replace('Ba{2+}','Ba++')\n",
    "    CCS = CCS.replace('O{-}','[O-]')  \n",
    "    CCS = CCS.replace('S{-}','[S-]')  \n",
    "    CCS = CCS.replace('Cl{-}','[Cl-]')        \n",
    "    CCS = CCS.replace('Br{-}','[Br-]')\n",
    "    #CCS = CCS.replace('n{+}','[nH+]')\n",
    "    CCS = CCS.replace('K','[K]')\n",
    "    \n",
    "    #CCS = CCS.replace('S{-}','[S-]')\n",
    "    CCS = CCS.replace('_','')\n",
    "    #CCS = CCS.replace('{','[')\n",
    "    #CCS = CCS.replace('}',']')    \n",
    "    return CCS\n",
    "def getdeg(dur = '28Days', classMode = 60, dr=dr):\n",
    "\n",
    "    ##Gets all molecules according to ClassMode. Doesn't classify yet, but includes other classModes that can be classified.\n",
    "\n",
    "    deg = {}\n",
    "    for i,row in dr[dr.Duration==dur].iterrows():\n",
    "        try:\n",
    "            CCS = fix_CCS(row.SMILES)\n",
    "            deg[CCS] = float(row.Data.rstrip('%'))\n",
    "        except:\n",
    "            None\n",
    "    if classMode:\n",
    "        for i, row in dr[dr.Duration!=dur].iterrows():\n",
    "            try:\n",
    "                CCS = fix_CCS(row.SMILES)\n",
    "                doyg = float(row.Data.rstrip('%'))\n",
    "                if doyg >= classMode and dur==\"28Days\":\n",
    "                    #print(row.Duration)\n",
    "                    deg[CCS] = doyg\n",
    "                elif doyg < classMode and dur=='14Days':\n",
    "                    #print(row.Duration)\n",
    "                    deg[CCS] = doyg\n",
    "            except:\n",
    "                None\n",
    "    print(sum(np.array(list(deg.values()))<classMode), sum(np.array(list(deg.values()))>=classMode),len(deg))\n",
    "    return deg\n",
    "def CalcFromsdf(path):\n",
    "    calc = {\"SMILES\":[], \"MolMethods\":[], \"_3Descs\": [], \"CalcMORSE\": [], \"CalcRDF\": [], \"CalcWHIM\": [], \"MACCS\": [], \"GCT\": []}\n",
    "    f=[]\n",
    "    #f = ['XLogP', 'ExactMass', 'RotatableBondCount','BondStereoCount', 'YStericQuadrupole3D', 'ZStericQuadrupole3D', \"ConformerModelRMSD3D\", 'EffectiveRotorCount3D', 'ConformerCount3D']\n",
    "    #f = [\"HBondDonorCount\", \"HeavyAtomCount\", \"FeatureDonorCount3D\", \"FeatureCationCount3D\", 'XLogP', \"FeatureRingCount3D\",\"FeatureHydrophobeCount3D\",\"ConformerModelRMSD3D\",] #ExactMass\", \"TPSA\", \"YStericQuadrupole3D\"]\n",
    "    with Chem.SDMolSupplier(path) as wdd:\n",
    "        for mol in tq(wdd):\n",
    "            #if i>10: break\n",
    "            try:\n",
    "                #need a .values + concatenation of da and dp\n",
    "                Q=[]\n",
    "                CCS = MolToSmiles(mol)\n",
    "                CCS = CCS.replace('/','')\n",
    "                CCS = CCS.replace('\\\\','')\n",
    "                CCS = CCS.replace('@','')\n",
    "                if f:\n",
    "                    a_list = da[da.SMILES==CCS]\n",
    "                for j in f:\n",
    "                    Q += [float(a_list[j])]\n",
    "                '''dps = dp[dp.SMILES==CCS]\n",
    "                    #Y.append(Descriptors.NumValenceElectrons(MolFromSmiles(CCS)))\n",
    "                    for j in prs:\n",
    "                        Q += [float(dps[j])]'''\n",
    "                ad = Chem.FindMolChiralCenters(mol,force=True,includeUnassigned=True)\n",
    "                #if ad:\n",
    "                    #print(CCS, ad)\n",
    "                calc[\"SMILES\"] += [CCS]\n",
    "                ph = []\n",
    "                for k in MolMethods:\n",
    "                    MolM = getattr(Chem.GraphDescriptors,k)\n",
    "                    ph += [MolM(mol)]\n",
    "                calc[\"MolMethods\"] += [ph]\n",
    "                ph = []\n",
    "                for k in _3Descs:\n",
    "                    MolM = getattr(Chem.rdMolDescriptors,k)\n",
    "                    ph += [MolM(mol)]\n",
    "                calc[\"_3Descs\"] += [ph]\n",
    "                for k in multi_3Descs:\n",
    "                    pass\n",
    "                    MolM = getattr(Chem.rdMolDescriptors,k)\n",
    "                    #print(k,len(MolM(mol)))\n",
    "                    calc[k] += [MolM(mol)]\n",
    "                ph = [k for k in MACCSkeys.GenMACCSKeys(Chem.MolFromSmiles(CCS))]\n",
    "                calc[\"MACCS\"] += [ph]\n",
    "                #par = Chem.RDKFingerprint(Chem.MolFromSmiles(CCS))\n",
    "                #pairFps = Pairs.GetAtomPairFingerprint(Chem.MolFromSmiles(CCS))\n",
    "                #print(pairFps.GetNonzeroElements())\n",
    "                ph = []\n",
    "                for k in fnames[\"GCT\"]:\n",
    "                    ph += [len(mol.GetSubstructMatches(Chem.MolFromSmarts(k)))]\n",
    "                calc[\"GCT\"] += [ph]\n",
    "            except:\n",
    "                None\n",
    "    \n",
    "    #all fingerprints to be SVD\n",
    "    '''        \n",
    "    svd = TruncatedSVD(n_components=40)\n",
    "    svd.fit(Xd)\n",
    "    Xd  = svd.transform(Xd)\n",
    "    for i in range(len(X)):\n",
    "        X[i] = list(Xd[i]) + X[i]\n",
    "    '''\n",
    "    for d in calc:\n",
    "        print(d,len(calc[d]))\n",
    "    return calc\n",
    "def CalcFromDB(path):\n",
    "    calc = {\"SMILES\":[], \"MolMethods\":[], \"_3Descs\": [], \"CalcMORSE\": [], \"CalcRDF\": [], \"CalcWHIM\": [], \"MACCS\": [], \"GCT\": []}\n",
    "    from ast import literal_eval\n",
    "    dcc = pd.read_excel(path, skiprows = 0)\n",
    "    dcc = dcc.drop(columns=[\"Unnamed: 0\"])\n",
    "    calc['SMILES'] = dcc.SMILES.to_list()\n",
    "    for p, r in dcc.iteritems():\n",
    "        if p=='SMILES': continue\n",
    "        if p=='MACCS': continue\n",
    "        calc[p] = [literal_eval(q) for q in r]\n",
    "        print(p,len(calc[p]))\n",
    "    #for p in dcc:\n",
    "        #print(pd.eval(p))\n",
    "    #calc = dcc.to_dict('records')\n",
    "    #print(pd.eval(calc[0][\"CalcMORSE\"]))\n",
    "    print(calc[\"MolMethods\"][1][9])\n",
    "    #fnames[\"MolMethods\"][9]\n",
    "    #print(dcc.SMILES.to_list())\n",
    "    return calc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d13e0143",
   "metadata": {},
   "source": [
    "### Write calc to .xlsx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "828d320a",
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = CalcFromsdf(path=\"c7.sdf\") #calc - oldconformers\n",
    "zzmols = set(calc[\"SMILES\"])\n",
    "print(len(zzmols))\n",
    "zzmolconfs = []\n",
    "for u in zzmols:\n",
    "    zzmolconfs.append(calc[\"SMILES\"].count(u))\n",
    "plt.hist(zzmolconfs,bins=50)\n",
    "plt.xlabel('Conformers used')\n",
    "plt.ylabel(\"No. of molecules\")\n",
    "calc['CalcWHIM']\n",
    "dcc = pd.DataFrame(data=calc)\n",
    "dcc.to_excel(\"c7.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cae20f",
   "metadata": {},
   "source": [
    "### SMILES reading, RDKit conformer generation (to .sdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ca7ee84",
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = pd.read_csv('Canonicalized_SMILES.csv')\n",
    "deg_old = {}\n",
    "for c in dx['SMILES']:\n",
    "    c =  MolToSmiles(MolFromSmiles(c))\n",
    "    deg_old[c]=0\n",
    "_new = []\n",
    "for i in deg_new:\n",
    "    try:\n",
    "        deg_old[i]\n",
    "    except:\n",
    "        _new.append(i)\n",
    "deg_old,len(_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8c5af1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for k in dir(Chem.GraphDescriptors):\n",
    "# ww.append(k)\n",
    "\n",
    "with Chem.SDWriter('c9.sdf') as w:\n",
    "    for m in tq(deg):\n",
    "        try:\n",
    "            for a in range(1):\n",
    "                mol = Chem.MolFromSmiles(m)\n",
    "                mol = Chem.AddHs(mol)\n",
    "                AllChem.EmbedMolecule(mol)\n",
    "                AllChem.MMFFOptimizeMolecule(mol)\n",
    "                mol = Chem.RemoveHs(mol)\n",
    "                # display(mol)\n",
    "                w.write(mol)\n",
    "                # previously our classes were 1:2, balance the former\n",
    "                if (deg[m] > classMode):\n",
    "                    mol = Chem.MolFromSmiles(m)\n",
    "                    mol = Chem.AddHs(mol)\n",
    "                    AllChem.EmbedMolecule(mol)\n",
    "                    AllChem.MMFFOptimizeMolecule(mol)\n",
    "                    mol = Chem.RemoveHs(mol)\n",
    "                    # display(mol)\n",
    "                    w.write(mol)\n",
    "        except:\n",
    "            None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f52c2312",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#test mols\n",
    "with Chem.SDMolSupplier('c8.sdf') as w:\n",
    "    k=0\n",
    "    for a in w:\n",
    "        print(Chem.MolToSmiles(a))\n",
    "        try:\n",
    "            print(deg[Chem.MolToSmiles(a)])\n",
    "        except:\n",
    "            None\n",
    "        display(a)\n",
    "        k=k+1\n",
    "        if (k>20): break"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e0e625e2",
   "metadata": {},
   "source": [
    "\n",
    "### Clustering (OLD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a01186f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fps={}\n",
    "with Chem.SDMolSupplier('c8.sdf') as wdd:\n",
    "    for mol in tq(wdd):\n",
    "        try:\n",
    "            CCS = MolToSmiles(mol)\n",
    "            fps[CCS] = AllChem.GetMorganFingerprintAsBitVect(mol,2,nBits=1024)\n",
    "            #fps[CCS] = MACCSkeys.GenMACCSKeys(mol)\n",
    "        except:\n",
    "            None\n",
    "    sim = [[1-DataStructs.FingerprintSimilarity(x, y) for y in list(fps.values())] for x in list(fps.values())]\n",
    "dist_sim=[]\n",
    "for i in range(len(sim)):\n",
    "    dist_sim += sim[i][0:i]\n",
    "print(len(fps),len(sim))\n",
    "dist_sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acc407cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "clust = Butina.ClusterData(dist_sim,len(sim),0.6,1)\n",
    "print([len(k) for k in clust])\n",
    "for u in clust:\n",
    "    for v in u:\n",
    "        print(list(fps.keys())[v])\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91429fcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sim = np.array(sim)\n",
    "#plt.plot(sim[0:12].T)\n",
    "from rdkit.ML import Cluster\n",
    "clust = Murtagh.ClusterData(dist_sim,len(sim),Murtagh.WARDS,1)\n",
    "for r in clust[0].GetChildren():\n",
    "    print(r.Print())\n",
    "ClusterUtils.GetNodesDownToCentroids(clust[0], above=1)\n",
    "clusting=[]\n",
    "for r in ClusterUtils.SplitIntoNClusters(clust[0],10):\n",
    "    clusting.append([k.GetData() for k in r.GetPoints()])\n",
    "for k in clusting:\n",
    "    print(len(k))\n",
    "for u in clusting:\n",
    "    for v in u:\n",
    "        print(list(fps.keys())[v])\n",
    "    print(\"\\n\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cf87129e",
   "metadata": {},
   "source": [
    "### *DATA LOADER - .xlsx to calc[]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71125750",
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = CalcFromDB(path=\"calc7.xlsx\") #calc - oldconformers\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a63fb10",
   "metadata": {},
   "outputs": [],
   "source": [
    "zzmols = set(calc[\"SMILES\"])\n",
    "print(len(zzmols))\n",
    "zzmolconfs = []\n",
    "for u in zzmols:\n",
    "    zzmolconfs.append(calc[\"SMILES\"].count(u))\n",
    "plt.hist(zzmolconfs,bins=50)\n",
    "plt.xlabel('Conformers used')\n",
    "plt.ylabel(\"No. of molecules\")\n",
    "calc['CalcWHIM']\n",
    "plt.hist(zzmolconfs,bins=27)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418501ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_dr(dr):\n",
    "    deg_new={}\n",
    "    for i,row in dr.iterrows():\n",
    "        try:\n",
    "            CCS = row.SMILES.replace('/','')\n",
    "            CCS = CCS.replace('\\\\','')\n",
    "            CCS = CCS.replace('@','')\n",
    "            CCS = MolToSmiles(MolFromSmiles(CCS))\n",
    "            deg_new[CCS] = float(row.Data.rstrip('%'))\n",
    "        except:\n",
    "            None\n",
    "    return deg_new\n",
    "def check_zzmols(zzmols):\n",
    "    zzmol = {}\n",
    "    for CC in zzmols:\n",
    "        CCS = CC.replace('/','')\n",
    "        CCS = CCS.replace('\\\\','')\n",
    "        CCS = CCS.replace('@','')\n",
    "        CCS = MolToSmiles(MolFromSmiles(CCS))\n",
    "        try:\n",
    "            print(CC,CCS,zzmol[CCS])\n",
    "        except:\n",
    "            zzmol[CCS]=0\n",
    "    print(len(zzmol), len(zzmols))\n",
    "    return zzmol\n",
    "    k=0\n",
    "def check_deg(deg,zzmol):\n",
    "    for z,CC in enumerate(deg.keys()):\n",
    "        CCS = CC.replace('/','')\n",
    "        CCS = CCS.replace('\\\\','')\n",
    "        CCS = CCS.replace('@','')\n",
    "        CCS = MolToSmiles(MolFromSmiles(CCS))\n",
    "        assert(CC==CCS)\n",
    "        try:\n",
    "            zzmol[CCS]\n",
    "        except:\n",
    "            print(CCS)\n",
    "            k += 1\n",
    "    return k"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d22ab9f9",
   "metadata": {},
   "source": [
    "### Set up the x and y from calc[]\n",
    "- Get the pre-filtered matrix (Mansouri 2013) ('fB') parameters\n",
    "- fix CCS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb267f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp={'CalcWHIM': [False, False, False, False, False, False, False, False, True, True, False, False, False, False, False, False, False, False, False, True, False, False, False, False, False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, True, True, False, False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False], 'CalcRDF': [False, False, True, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, True, True, True, False, True, False, False, False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, True, True, False, False, False, False, False, False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False, False, True, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, True, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, False, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False], 'CalcMORSE': [False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False, False, False, False, True, True, False, False, False, True, False, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, True, False, True, True, False, False, False, True, False, False, False, True, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, False, False, True, False, False, False, True, False, False, False, False, False, False, False, False, False, True, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False, False, False, False, False, False, False, True, False, False, True, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, False, True, False]}\n",
    "getFeatureMask = {'CalcWHIM': [8, 9, 19, 28, 41, 42, 48, 66],\n",
    "'CalcRDF': [2, 5, 31, 34, 35, 36, 38, 45, 62, 64, 65, 75, 87, 92, 94, 121, 123, 124, 155, 156, 158, 182, 183, 188],\n",
    "'CalcMORSE': [2, 32, 39, 40, 44, 48, 62, 64, 66, 67, 71, 75, 77, 108, 109, 112, 116, 126, 128, 162, 163, 192, 200, 203, 222]\n",
    "}\n",
    "for k in getFeatureMask:\n",
    "    #gamma = [x for x,i in enumerate(sp[k]) if i]\n",
    "    #print(k,gamma,len(fnames[k]))\n",
    "    sp[k] = [True if i in getFeatureMask[k] else False for i in range(len(fnames[k]))]\n",
    "deg_new={}\n",
    "rej=[]\n",
    "for i,row in dr.iterrows():\n",
    "        try:\n",
    "            CCS = fix_CCS(row.SMILES)\n",
    "            CCS = MolToSmiles(MolFromSmiles(CCS))\n",
    "            deg_new[CCS] = float(row.Data.rstrip('%'))\n",
    "        except:\n",
    "            rej.append([CCS,row.SMILES])\n",
    "len(deg_new),len(rej)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39338ab6",
   "metadata": {},
   "outputs": [],
   "source": [
    "classMode = 60\n",
    "\n",
    "deg = getdeg('28Days', classMode, dr)\n",
    "X = []\n",
    "Y = []\n",
    "g = []\n",
    "SMILEScount = 1\n",
    "samecount = 1\n",
    "last_CCS = 0\n",
    "rej_CCS = set()\n",
    "\n",
    "elimratio=sum(np.array(list(deg.values()))>classMode)/sum(np.array(list(deg.values()))<classMode)\n",
    "for i, CCS in enumerate(tq(calc[\"SMILES\"])):\n",
    "    CCS = CCS.replace('/','')\n",
    "    CCS = CCS.replace('\\\\','')\n",
    "    CCS = CCS.replace('@','')\n",
    "    #m1 = MolFromSmiles(CCS)\n",
    "    #m1 = SaltRemover().StripMol(m1)\n",
    "    fA = [\"GCT\",\"MolMethods\",\"_3Descs\"] #GETAWAY_entries\n",
    "    fB = [\"CalcMORSE\", \"CalcRDF\", \"CalcWHIM\"]\n",
    "    CCS = MolToSmiles(MolFromSmiles(CCS))\n",
    "    try:\n",
    "        Q=[]\n",
    "        #Y.append(calc[\"MolMethods\"][i][1])#float(da.XLogP[i]))#YStericQuadrupole3D\n",
    "        #untagged CCS do not run below     \n",
    "        if 1:\n",
    "            if last_CCS==CCS:\n",
    "                samecount += 1 #Q = X[-1]\n",
    "            else:\n",
    "                samecount = 1\n",
    "                SMILEScount += 1\n",
    "            if samecount > 5:\n",
    "                #continue\n",
    "                samecount = 1\n",
    "        last_CCS = CCS\n",
    "        if deg[CCS]<classMode and np.random.rand()>0.67:#and samecount>1 :\n",
    "            continue\n",
    "            Y.append(deg[CCS])\n",
    "            g.append(CCS)\n",
    "            X.append(Q)\n",
    "        for k in fA:\n",
    "            Q += calc[k][i]\n",
    "        for k in fB:\n",
    "            Q += np.array(calc[k][i])[sp[k]].tolist()\n",
    "        Y.append(deg[CCS])\n",
    "        g.append(CCS)\n",
    "        X.append(Q)\n",
    "    except:\n",
    "        rej_CCS.add(CCS)\n",
    "\n",
    "g = [{v: k for k, v in enumerate(OrderedDict.fromkeys(g))} [n] for n in g]\n",
    "f = []\n",
    "for k in fA:\n",
    "    print(k)\n",
    "    f += fnames[k]\n",
    "\n",
    "for k in fB:\n",
    "    print(k)\n",
    "    f += np.array(fnames[k])[sp[k]].tolist()\n",
    "\n",
    "sr = VarianceThreshold(1e-2)\n",
    "X = sr.fit_transform(X)\n",
    "f = sr.get_feature_names_out(f)\n",
    "\n",
    "plt.hist(Y,bins=60)#,range=(10,110))\n",
    "plt.xticks(np.arange(0,max(Y),10))\n",
    "plt.xlabel(\"Biochemical Oxygen Demand\")\n",
    "plt.ylabel(\"Number of compounds\")\n",
    "plt.show()\n",
    "#archive\n",
    "X,Y,f,g = np.array(X),np.array(Y),np.array(f),np.array(g)\n",
    "arcX, arcY, arcf, arcg = X,Y,f,g\n",
    "\n",
    "X.shape,Y.shape,Y.mean(),g[:100],sum(Y>classMode)/sum(Y<classMode),SMILEScount,len(set(g)),len(rej_CCS)\n",
    "#pd.DataFrame(X,index=g,columns=f).head(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6661cbf8",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def InitData(classMode, cors_th, nm_x=True):\n",
    "    mask = None\n",
    "    g = arcg\n",
    "    if classMode:\n",
    "        Y = np.array([min(1, int((a)/classMode)) for a in arcY])\n",
    "    else:\n",
    "        Y = arcY\n",
    "    #correlation filter\n",
    "    if cors_th:\n",
    "        cors = abs(np.corrcoef(np.transpose(arcX)))\n",
    "        to_rm = set()                                \n",
    "        for a in cors:  # eliminate correlations\n",
    "            if (a > cors_th).sum() > 1:\n",
    "                for c in np.nonzero(a > cors_th):\n",
    "                    for b in c:\n",
    "                        if cors[b].sum() > a.sum():\n",
    "                            to_rm.add(b)\n",
    "        mask = [i not in to_rm for i in range(len(arcf))]\n",
    "        X = arcX[:, mask]\n",
    "        f = arcf[mask]\n",
    "    else:\n",
    "        X = arcX\n",
    "        f = arcf\n",
    "    if nm_x:\n",
    "        scaler = StandardScaler().fit(X)\n",
    "        X = scaler.transform(X)\n",
    "    return X, Y, f, mask\n",
    "\n",
    "\n",
    "X, Y, f, mask = InitData(60, cors_th=0.97)\n",
    "X, Y, g = shuffle(X, Y, g)\n",
    "# re-number g\n",
    "#g = [{v: k for k, v in enumerate(OrderedDict.fromkeys(g))}[n] for n in tq(g)]\n",
    "#g = np.array(g)\n",
    "# X(f*g), Y(g)\n",
    "X.shape, Y.shape, g.shape, X.mean(), Y.mean(), len(f), g[:100], #Y[0:100], X[101:102], f"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "659f9cde",
   "metadata": {},
   "source": [
    "### Chemprop\n",
    "- ONLY .csv files from this point are produced by Chemprop!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5b1eda",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install chemprop==1.6.1.dev0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e79380",
   "metadata": {},
   "outputs": [],
   "source": [
    "##only for generating the csv\n",
    "classMode = 60\n",
    "deg = getdeg('28Days', classMode, dr)\n",
    "feas = []\n",
    "for CCS in deg:\n",
    "    \n",
    "    #CCS = CCS.replace('/','')\n",
    "    #CCS = CCS.replace('\\\\','')\n",
    "    #CCS = CCS.replace('@','')\n",
    "    CSSS = fix_CCS(CCS)\n",
    "    #try:\n",
    "    #v = calc[calc['SMILES'] == CCS]\n",
    "    #print(v)\n",
    "    if(deg[CCS]<60):\n",
    "        deg[CCS] = 0\n",
    "    else:\n",
    "        deg[CCS] = 1\n",
    "    # except:\n",
    "    #     None\n",
    "dcc = pd.DataFrame.from_dict(data=deg,orient=\"index\")\n",
    "dcc.to_csv(\"cr_r.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f95bcc4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [\n",
    "    '--data_path', 'cr1.csv',\n",
    "    '--dataset_type', 'regression', #<----------------------------\n",
    "    '--metric', 'r2',\n",
    "    '--num_folds', '5',\n",
    "    '--save_dir', 'cr_checkpoints'\n",
    "]\n",
    "\n",
    "args = chemprop.args.TrainArgs().parse_args(arguments)\n",
    "mean_score, std_score = chemprop.train.cross_validate(args=args, train_func=chemprop.train.run_training)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb24974a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "arguments = [\n",
    "    '--config_save_path', 'hp.json',\n",
    "    '--data_path', 'cr1.csv',\n",
    "    '--dataset_type', 'classification',\n",
    "    '--metrics', 'f1'\n",
    "]\n",
    "args = chemprop.args.HyperoptArgs().parse_args(arguments)\n",
    "chemprop.hyperparameter_optimization.hyperopt(args=args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c460b958",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['morgan', 'morgan_count', 'rdkit_2d', 'rdkit_2d_normalized']"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chemprop.features.features_generators.get_available_features_generators()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2b25a1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#import chemprop\n",
    "arguments = [\n",
    "    '--config_path', 'hp.json',\n",
    "    '--data_path', 'cr_r.csv',\n",
    "    '--dataset_type', 'classification',\n",
    "    '--features_generator', 'rdkit_2d_normalized',\n",
    "    '--no_features_scaling',\n",
    "    '--metric', 'accuracy',\n",
    "    '--quiet',\n",
    "    '--class_balance',\n",
    "    '--num_folds', '5',\n",
    "    '--save_dir', 'cu_checkpoints'\n",
    "]\n",
    "args = chemprop.args.TrainArgs().parse_args(arguments)\n",
    "#args.features_generator = ['morgan_count', 'rdkit_2d_normalized']\n",
    "args.extra_metrics = ['balanced_accuracy',  'recall', 'precision', 'f1','auc']\n",
    "#print(args.extra_metrics)\n",
    "mean_score, std_score = chemprop.train.cross_validate(args=args, train_func=chemprop.train.run_training)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48ddbeae",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [\n",
    "    '--test_path', 'cza.csv',\n",
    "    '--preds_path', 'czp.csv',\n",
    "    '--checkpoint_dir', 'cu_checkpoints'\n",
    "]\n",
    "args = chemprop.args.PredictArgs().parse_args(arguments)\n",
    "preds = chemprop.train.make_predictions(args=args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a668a70",
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = [\n",
    "    '--data_path', 'ca.csv',\n",
    "    '--checkpoint_dir', 'cd_checkpoints'\n",
    "]\n",
    "args = chemprop.args.InterpretArgs().parse_args(arguments)\n",
    "#interprets = chemprop.train.evaluate_predictions\n",
    "chemprop.interpret()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "03a28eca",
   "metadata": {},
   "source": [
    "## Training and Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "87a70289",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import svm\n",
    "from sklearn import tree\n",
    "from sklearn.naive_bayes import GaussianNB, MultinomialNB\n",
    "from sklearn.linear_model import LinearRegression, LogisticRegression\n",
    "from sklearn.feature_selection import RFE, RFECV, SelectFromModel\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier,RandomForestRegressor\n",
    "from sklearn.model_selection import train_test_split, cross_val_score, GroupKFold, StratifiedGroupKFold, GroupShuffleSplit #we don't need GroupStrat for now\n",
    "from sklearn.metrics import classification_report, balanced_accuracy_score, confusion_matrix\n",
    "#import seaborn as sns\n",
    "def udr(embedding,Y):\n",
    "    fig, axs = plt.subplots(1,embedding.shape[1]-1)\n",
    "    for p in range(embedding.shape[1]-1):\n",
    "        axs[p].scatter(embedding[:, 0],embedding[:, p+1],c=Y,marker=0)\n",
    "    plt.gca().set_aspect('equal', 'datalim')\n",
    "    plt.show()\n",
    "    return #embedding\n",
    "    #plt.title('UMAP projection of the Penguin dataset', fontsize=24)\n",
    "def udr_old(X, Y=None, min_dist=0.1, n_neighbors=15, n_components=2):\n",
    "    reducer = umap.UMAP(a=None, angular_rp_forest=False, b=None,\n",
    "         force_approximation_algorithm=False, init='spectral', learning_rate=1.0,\n",
    "         local_connectivity=1.0, low_memory=False, metric='euclidean',\n",
    "         metric_kwds=None, min_dist=min_dist, n_components=n_components, n_epochs=None,\n",
    "         n_neighbors=n_neighbors, negative_sample_rate=5, output_metric='euclidean',\n",
    "         output_metric_kwds=None, random_state=None, repulsion_strength=1.0,\n",
    "         set_op_mix_ratio=1.0, spread=1.0, target_metric='categorical',\n",
    "         target_metric_kwds=None, target_n_neighbors=-1, target_weight=0.5,\n",
    "         transform_queue_size=4.0, transform_seed=None, unique=False, verbose=False)\n",
    "    Xw = StandardScaler().fit_transform(X)\n",
    "    embedding = reducer.fit_transform(Xw,Y)\n",
    "    embedding.shape\n",
    "    plt.scatter(embedding[:, 0],embedding[:, 1],c=Y,marker=0)\n",
    "    plt.gca().set_aspect('equal', 'datalim')\n",
    "    plt.show()\n",
    "    return embedding\n",
    "    #plt.title('UMAP projection of the Penguin dataset', fontsize=24)\n",
    "def tsn(X, min_dist=12, n_neighbors=30, n_components=2):\n",
    "    reducer = manifold.TSNE(n_components=n_components, perplexity=n_neighbors, early_exaggeration=min_dist, learning_rate='warn',\n",
    "                  n_iter=1000, n_iter_without_progress=300, min_grad_norm=1e-07, metric='euclidean',\n",
    "                  metric_params=None, init='warn', verbose=0, random_state=None, method='barnes_hut',\n",
    "                  angle=0.5, n_jobs=None, square_distances='deprecated')\n",
    "    Xw = StandardScaler().fit_transform(X)\n",
    "    embedding = reducer.fit_transform(Xw,Y)\n",
    "    #embedding.shape\n",
    "    plt.scatter(embedding[:, 0],embedding[:, 1],c=Y,marker=0)\n",
    "    plt.gca().set_aspect('equal', 'datalim')\n",
    "    plt.show()\n",
    "    return embedding\n",
    "def tts(X,Y,g,test_size=0.25,random_state=0):\n",
    "    train_indx, test_indx = next(GroupShuffleSplit(n_splits=1, test_size=test_size, random_state=random_state).split(X,Y,groups=g))\n",
    "    #np.unique(g[train_indx]), np.unique(g[test_indx])\n",
    "    #I continue to keep W for legacy purposes\n",
    "    return X[train_indx], X[test_indx], Y[train_indx], Y[test_indx], g[train_indx], g[test_indx]\n",
    "#udr_old(X,Y,0.1,15,4)\n",
    "#tsn(X,12,30,2)\n",
    "def confusion(X_test,X_train,Y_test,Y_train):\n",
    "    Y1 = sklearn_model.predict(X_test)\n",
    "    print(confusion_matrix(Y1,Y_test))\n",
    "    print(classification_report(Y1,Y_test))\n",
    "    print(balanced_accuracy_score(Y1,Y_test), sklearn_model.score(X_test,Y_test))\n",
    "\n",
    "    Y1 = sklearn_model.predict(X_train)\n",
    "    print(confusion_matrix(Y1,Y_train))\n",
    "    print(classification_report(Y1,Y_train))\n",
    "    print(balanced_accuracy_score(Y1,Y_train), sklearn_model.score(X_train,Y_train))\n",
    "confusion(X_test,X_train,Y_test,Y_train)\n",
    "\n",
    "#print(sklearn_model.best_estimator_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d8c1a89",
   "metadata": {},
   "outputs": [],
   "source": [
    "X,Y,f=InitData(classMode)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1fc3bafd",
   "metadata": {},
   "source": [
    "### Feature Reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5d95f90",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#dataset = dc.data.NumpyDataset(X, Y)\n",
    "\n",
    "scaler = StandardScaler().fit(X)\n",
    "X = scaler.transform(X)\n",
    "sklearn_model = LogisticRegression(max_iter=500000,class_weight='balanced')\n",
    "sklearn_model = svm.SVC(kernel='linear')#, class_weight='balanced')\n",
    "#sklearn_model = RandomForestClassifier(n_estimators=50,n_jobs=-1,max_features=\"sqrt\",random_state=0)\n",
    "rfecv = RFECV(\n",
    "    estimator = sklearn_model,\n",
    "    step = 2,\n",
    "    cv = 5,\n",
    "    scoring = \"balanced_accuracy\",\n",
    "    #importance_getter = 'coef_',\n",
    "    min_features_to_select = 1,\n",
    "    n_jobs = -1,\n",
    "    verbose = True\n",
    ")\n",
    "rfecv.fit(X, Y)\n",
    "\n",
    "print(\"Optimal number of features : %d\" % rfecv.n_features_)\n",
    "\n",
    "# Plot number of features VS. cross-validation scores\n",
    "plt.figure()\n",
    "plt.title(\"Random Forest\")\n",
    "plt.xlabel(\"Number of features selected\")\n",
    "plt.ylabel(\"Cross validation score (accuracy)\")\n",
    "plt.plot(\n",
    "    range(1, rfecv.step*len(rfecv.grid_scores_) + 1,rfecv.step),\n",
    "    rfecv.grid_scores_,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0d7c915",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['[CH3,CH4]' '*-C(=O)-*' '*-[OR0]-*' '[O]=[C]-[O]-*' '[N]=*'\n",
      " '*-[CHR](-*)-*' '*-[CH2]-[OH]' 'n' '*-C(=[!O&!S])-*' '*-[CR](=O)-*'\n",
      " '[R][N](=O)O' '*-[S]-*' '[R0][C](=O)[OH]' '*-S=O' '*-[CH]-[OH]'\n",
      " '*-[CH0R](-*)(-*)-*' '[F][R0]' '[Br][R0]' '[Br][R]'\n",
      " '[CH2][CH2][C](=O)[OH]' '[NH2][CH]C=O' 'BalabanJ' 'BertzCT' 'Chi1'\n",
      " 'Chi1v' 'Chi2v' 'Chi3n' 'Chi3v' 'Chi4n' 'HallKierAlpha' 'Kappa1' 'Kappa2'\n",
      " 'CalcAsphericity' 'CalcFractionCSP3' 'CalcNPR2' 'CalcPBF' 'CalcPMI2'\n",
      " 'CalcPhi' 'MORSE_2u' 'MORSE_0m' 'MORSE_7m' 'MORSE_12m' 'MORSE_30m'\n",
      " 'MORSE_3v' 'MORSE_7v' 'MORSE_12e' 'MORSE_30e' 'MORSE_0p' 'MORSE_3i'\n",
      " 'MORSE_0s' 'RDF_30m' 'RDF_40m' 'RDF_85m' 'RDF_20v' 'RDF_30v' 'RDF_145v'\n",
      " 'RDF_15p' 'RDF_25p' 'RDF_30p' 'RDF_35i' 'RDF_40i' 'RDF_50i' 'RDF_20s'\n",
      " 'RDF_50s' 'WHIM_66'] 65\n"
     ]
    }
   ],
   "source": [
    "sklearn_model = LogisticRegression(max_iter=500000,class_weight='balanced')\n",
    "\n",
    "#sklearn_model = RandomForestClassifier(n_estimators=20,n_jobs=-1,max_features=\"sqrt\",random_state=0)\n",
    "#selector = SelectFromModel(sklearn_model, max_features=50, threshold=-np.inf)# threshold = \"0.5*mean\")\n",
    "selector = RFE(sklearn_model, n_features_to_select = 65)\n",
    "\n",
    "selector = selector.fit(X, Y)\n",
    "#for j in range(len(f)):\n",
    "    #print(selector.ranking_[j])\n",
    "X = selector.transform(X)\n",
    "f = selector.get_feature_names_out(f)\n",
    "print(f,len(f))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed3e1ca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "rfecv = RFE(sklearn_model, n_features_to_select = 64)\n",
    "rfecv = rfecv.fit(X, Y)\n",
    "print(f, rfecv.get_support())\n",
    "sp1 = rfecv.get_support()\n",
    "#mask using the original set,\n",
    "\n",
    "print(mask)\n",
    "ii = 0\n",
    "for i, r in enumerate(mask):\n",
    "    if(r==True):\n",
    "        mask[i] = sp1[ii]\n",
    "        ii+=1\n",
    "#sp[\"CalcMORSE\"] = mask\n",
    "#print(mask)\n",
    "sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fa20115",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CalcWHIM 8 114\n",
      "CalcRDF 24 210\n",
      "CalcMORSE 25 224\n"
     ]
    }
   ],
   "source": [
    "for k in sp:\n",
    "    print(k,np.count_nonzero(sp[k]),np.size(sp[k]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "7c3fc51c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['MORSE_2u', 'MORSE_0m', 'MORSE_7m', 'MORSE_8m', 'MORSE_12m',\n",
       "       'MORSE_16m', 'MORSE_30m', 'MORSE_0v', 'MORSE_2v', 'MORSE_3v',\n",
       "       'MORSE_7v', 'MORSE_11v', 'MORSE_13v', 'MORSE_12e', 'MORSE_13e',\n",
       "       'MORSE_16e', 'MORSE_20e', 'MORSE_30e', 'MORSE_0p', 'MORSE_2i',\n",
       "       'MORSE_3i', 'MORSE_0s', 'MORSE_8s', 'MORSE_11s', 'MORSE_30s'],\n",
       "      dtype='<U9')"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calca = np.array(calc[\"CalcMORSE\"][1])\n",
    "calca[sp[\"CalcMORSE\"]]\n",
    "ff1 = np.array(fnames[\"CalcMORSE\"])\n",
    "ff1[sp[\"CalcMORSE\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da419a42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((14110, 129), (14110,))"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array(X).shape,np.array(Y).shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b1f6d609",
   "metadata": {},
   "source": [
    "#### Logistic Reg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d4ad0a2c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7864407354715832 [0.77405338 0.80136561 0.74736189 0.82929857 0.78012422]\n",
      "0.7868289391041138 [0.77418432 0.80525078 0.74540259 0.82931994 0.77998707]\n",
      "0.7864407354715832 [0.77405338 0.80136561 0.74736189 0.82929857 0.78012422]\n",
      "0.7909222466771666 [0.77434383 0.82015922 0.74716516 0.82948938 0.78345365]\n",
      "0.7858538795502281 [0.77410198 0.79946815 0.74698547 0.82927858 0.77943523]\n",
      "0.7805370498405397 [0.76962025 0.78378378 0.72629455 0.83097726 0.7920094 ]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1e779ff04f0>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAAsTAAALEwEAmpwYAAAS8ElEQVR4nO3df4xdZX7f8fcH/4C6KwVTZrcbzNpOoAJ2SyC6cQUUeeuULts2EDZpNchqo0SR/2hYqUSR4opW3bK7UuImYkWKunKrRFHlxWHdTWsUEWuzGCmrmq7H1A5rZg3GBGPXZWfT0q5rJQjy7R/3uFwPg+fOzB2PJ8/7JR3NOc/znMPzZaTzmXOeeyFVhSSpPVcs9QQkSUvDAJCkRhkAktQoA0CSGmUASFKjVi71BObi2muvrQ0bNiz1NCRpWTl06ND3qmpsevuyCoANGzYwMTGx1NOQpGUlyesztfsKSJIaZQBIUqMMAElqlAEgSY0yACSpUQaAtAC7dsGGDXDFFf2fu3Yt9Yyk4S2rj4FKl5Ndu2DbNjh3rn/8+uv9Y4CtW5duXtKwfAKQ5umRR967+Z937ly/XVoODABpnk6enFu7dLkxAKR5+tjH5tYuXW4MAGmevvhFWLPmwrY1a/rt0nJgAEjztHUr7NwJ69dD0v+5c6cLwFo+/BSQtABbt3rD1/LlE4AkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0aKgCS3JvkWJLjSbbP0P9YksPd9nKStwb6diQ5mmQyyeNJ0rWvTrKzG/+dJD81sqokSbOa9YtgSVYATwD3AKeAg0n2VtVL58dU1cMD4z8L3N7t3wncBdzadX8T2Aw8BzwCfLeq/lqSK4BrRlGQJGk4w3wTeBNwvKpOACTZDdwPvPQB4x8E/mW3X8BVwGogwCrgza7v54CbAKrqz4HvzWP+kqR5GuYV0HXAGwPHp7q290myHtgIPAtQVQeA/cCZbttXVZNJru5O+XySF5J8NclHPuCa25JMJJmYmpoapiZJ0hBGvQg8DuypqncBktwA3Aysox8aW5LcTf/JYx3wX6rqR4EDwK/NdMGq2llVvarqjY2NjXi6ktSuYQLgNHD9wPG6rm0m48CTA8cPAM9X1dmqOgs8A9wB/AlwDvhaN+6rwI/OYd6SpAUaJgAOAjcm2ZhkNf2b/N7pg5LcBKyl/9f8eSeBzUlWJllFfwF4sqoKeBr4ZDfux/ngNQVJ0iKYdRG4qt5J8hCwD1gB/GZVHU3yKDBRVefDYBzY3d3cz9sDbAFepL8g/PtV9XTX98vAf0jyJWAK+NlRFCRJGk4uvF9f3nq9Xk1MTCz1NCRpWUlyqKp609v9JrAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRQwVAknuTHEtyPMn2GfofS3K4215O8tZA344kR5NMJnk8Sbr257prnj/vwyOrSpI0q5WzDUiyAngCuAc4BRxMsreqXjo/pqoeHhj/WeD2bv9O4C7g1q77m8Bm4LnueGtVTSy8DEnSXA3zBLAJOF5VJ6rqbWA3cP9Fxj8IPNntF3AVsBq4ElgFvDn/6UqSRmWYALgOeGPg+FTX9j5J1gMbgWcBquoAsB840237qmpy4JTf6l7//Ivzr4ZmuOa2JBNJJqampoaYriRpGKNeBB4H9lTVuwBJbgBuBtbRD40tSe7uxm6tqr8O3N1t/2imC1bVzqrqVVVvbGxsxNOVpHYNEwCngesHjtd1bTMZ573XPwAPAM9X1dmqOgs8A9wBUFWnu5/fB75C/1WTJOkSGSYADgI3JtmYZDX9m/ze6YOS3ASsBQ4MNJ8ENidZmWQV/QXgye742u68VcDfB769sFIkSXMxawBU1TvAQ8A+YBJ4qqqOJnk0yX0DQ8eB3VVVA217gFeBF4EjwJGqepr+gvC+JH8EHKb/RPHvRlCPJGlIufB+fXnr9Xo1MeGnRiVpLpIcqqre9Ha/CSxJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjRoqAJLcm+RYkuNJts/Q/1iSw932cpK3Bvp2JDmaZDLJ40ky7dy9Sb694EokSXOycrYBSVYATwD3AKeAg0n2VtVL58dU1cMD4z8L3N7t3wncBdzadX8T2Aw81/V/Bjg7ikIkSXMzzBPAJuB4VZ2oqreB3cD9Fxn/IPBkt1/AVcBq4EpgFfAmQJIPAb8IfGF+U5ckLcQwAXAd8MbA8amu7X2SrAc2As8CVNUBYD9wptv2VdVkN/zzwK8D5y72D0+yLclEkompqakhpitJGsaoF4HHgT1V9S5AkhuAm4F19ENjS5K7k9wG/HBV/e5sF6yqnVXVq6re2NjYiKcrSe2adQ0AOA1cP3C8rmubyTjwCwPHDwDPV9VZgCTPAHcA3wd6Sf64m8OHkzxXVZ+c0+wlSfM2zBPAQeDGJBuTrKZ/k987fVCSm4C1wIGB5pPA5iQrk6yivwA8WVX/tqp+sKo2AH8TeNmbvyRdWrMGQFW9AzwE7AMmgaeq6miSR5PcNzB0HNhdVTXQtgd4FXgROAIcqaqnRzZ7SdK85cL79eWt1+vVxMTEUk9DkpaVJIeqqje93W8CS1KjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGjVUACS5N8mxJMeTbJ+h/7Ekh7vt5SRvDfTtSHI0yWSSx5Oka//9JEe6vi8nWTGyqiRJs5o1ALob8xPAp4FbgAeT3DI4pqoerqrbquo24DeAr3Xn3gncBdwKfAL4MWBzd9o/rKof6drHgH8wioIkScMZ5glgE3C8qk5U1dvAbuD+i4x/EHiy2y/gKmA1cCWwCngToKr+TzdmZddfc569JGnehgmA64A3Bo5PdW3vk2Q9sBF4FqCqDgD7gTPdtq+qJgfG7wO+C3wf2PMB19yWZCLJxNTU1BDTlSQNY9SLwOPAnqp6FyDJDcDNwDr6obElyd3nB1fVp4CP0n862DLTBatqZ1X1qqo3NjY24ulKUruGCYDTwPUDx+u6tpmM897rH4AHgOer6mxVnQWeAe4YPKGq/hT4z1z8tZIkacSGCYCDwI1JNiZZTf8mv3f6oCQ3AWuBAwPNJ4HNSVYmWUV/AXgyyYeSfLQ7byXw94DvLKwUSdJczBoAVfUO8BCwD5gEnqqqo0keTXLfwNBxYHdVDS7m7gFeBV4EjgBHqupp4C8De5P8EXCY/jrAl0dQjyRpSLnwfn156/V6NTExsdTTkKRlJcmhqupNb/ebwJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUqKECIMm9SY4lOZ5k+wz9jyU53G0vJ3lroG9HkqNJJpM8nr41SX4vyXe6vl8ZYU2SpCGsnG1AkhXAE8A9wCngYJK9VfXS+TFV9fDA+M8Ct3f7dwJ3Abd23d8ENgPfAn6tqvYnWQ18I8mnq+qZ0ZQlSZrNME8Am4DjVXWiqt4GdgP3X2T8g8CT3X4BVwGrgSuBVcCbVXWuqvYDdNd8AVg3vxIkSfMxTABcB7wxcHyqa3ufJOuBjcCzAFV1ANgPnOm2fVU1Oe2cq4GfAL7xAdfclmQiycTU1NQQ05UkDWPUi8DjwJ6qehcgyQ3AzfT/ur8O2JLk7vODk6yk/7TweFWdmOmCVbWzqnpV1RsbGxvxdCWpXcMEwGng+oHjdV3bTMZ57/UPwAPA81V1tqrOAs8Adwz07wReqaovDT1jSdJIDBMAB4Ebk2zsFmzHgb3TByW5CVgLHBhoPglsTrIyySr6C8CT3fgvAD8A/NMFVSBJmpdZA6Cq3gEeAvbRv3k/VVVHkzya5L6BoePA7qqqgbY9wKvAi8AR4EhVPZ1kHfAIcAvwQvfx0Z8fTUmSpGHkwvv15a3X69XExMRST0OSlpUkh6qqN73dbwJLUqMMAElqlAEgSY0yACSpUQaAJDXKAJCkRhkAktQoA0CSGmUASFKjDABJapQBIEmNMgAkqVEGgCQ1ygCQpEYZAJLUKANAkhplAEhSowwASWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaNVQAJLk3ybEkx5Nsn6H/sSSHu+3lJG8N9O1IcjTJZJLHk6Rr/2KSN5KcHVk1kqShzRoASVYATwCfBm4BHkxyy+CYqnq4qm6rqtuA3wC+1p17J3AXcCvwCeDHgM3daU8Dm0ZThiRproZ5AtgEHK+qE1X1NrAbuP8i4x8Enuz2C7gKWA1cCawC3gSoquer6sx8Jy5JWphhAuA64I2B41Nd2/skWQ9sBJ4FqKoDwH7gTLftq6rJuUwwybYkE0kmpqam5nKqJOkiRr0IPA7sqap3AZLcANwMrKMfGluS3D2XC1bVzqrqVVVvbGxsxNOVpHYNEwCngesHjtd1bTMZ573XPwAPAM9X1dmqOgs8A9wxn4lKkkZrmAA4CNyYZGOS1fRv8nunD0pyE7AWODDQfBLYnGRlklX0F4Dn9ApIkrQ4Zg2AqnoHeAjYR//m/VRVHU3yaJL7BoaOA7urqgba9gCvAi8CR4AjVfU0/P+Ph54C1iQ5leRzI6lIkjSUXHi/vrz1er2amJhY6mlI0rKS5FBV9aa3+01gSWqUASBJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIC/HaLvhPG+ArV/R/vrZrqWckDW3lUk9AWrZe2wXf2gbvnusfn3u9fwywcevSzUsakk8A0nwdeeS9m/95757rt0vLgAEgzde5k3Nrly4zBoA0X2s+Nrd26TJjAEjz9SNfhBVrLmxbsabfLi0DBoA0Xxu3wqadsGY9kP7PTTtdANay4aeApIXYuNUbvpYtnwAkqVEGgCQ1ygCQpEYZAJLUKANAkhq1rP6n8EmmgNeXeh5zdC3wvaWexCVmzW2w5uVjfVWNTW9cVgGwHCWZqKreUs/jUrLmNljz8ucrIElqlAEgSY0yABbfzqWewBKw5jZY8zLnGoAkNconAElqlAEgSY0yAEYgyTVJvp7kle7n2g8Y9zPdmFeS/MwM/XuTfHvxZ7xwC6k5yZokv5fkO0mOJvmVSzv7uUlyb5JjSY4n2T5D/5VJfqfr/69JNgz0/bOu/ViST13SiS/AfGtOck+SQ0le7H5uueSTn4eF/I67/o8lOZvkly7ZpEehqtwWuAE7gO3d/nbgV2cYcw1wovu5tttfO9D/GeArwLeXup7FrhlYA/ytbsxq4A+BTy91TR9Q5wrgVeCHurkeAW6ZNuafAF/u9seB3+n2b+nGXwls7K6zYqlrWuSabwd+sNv/BHB6qetZzHoH+vcAXwV+aanrmcvmE8Bo3A/8drf/28BPzjDmU8DXq+p/VtX/Ar4O3AuQ5EPALwJfWPypjsy8a66qc1W1H6Cq3gZeANYt/pTnZRNwvKpOdHPdTb/2QYP/LvYAP54kXfvuqvqzqnoNON5d73I375qr6r9V1X/v2o8CfynJlZdk1vO3kN8xSX4SeI1+vcuKATAaH6mqM93+/wA+MsOY64A3Bo5PdW0Anwd+HTi3aDMcvYXWDECSq4GfAL6xCHMchVlrGBxTVe8A/xv4K0OeezlaSM2Dfgp4oar+bJHmOSrzrrf74+2XgX91CeY5cv4fwYaU5A+AvzpD1yODB1VVSYb+bG2S24AfrqqHp79XXGqLVfPA9VcCTwKPV9WJ+c1Sl6MkHwd+Ffg7Sz2XRfY54LGqOts9ECwrBsCQqupvf1BfkjeTfLSqziT5KPDdGYadBj45cLwOeA64A+gl+WP6v48PJ3muqj7JElvEms/bCbxSVV9a+GwXzWng+oHjdV3bTGNOdaH2A8CfDHnu5WghNZNkHfC7wD+uqlcXf7oLtpB6/wbw00l2AFcDf57kT6vq3yz6rEdhqRch/iJswL/mwgXRHTOMuYb+e8K13fYacM20MRtYPovAC6qZ/nrHfwSuWOpaZqlzJf3F6428t0D48WljfoELFwif6vY/zoWLwCdYHovAC6n56m78Z5a6jktR77Qxn2OZLQIv+QT+Imz0331+A3gF+IOBm1wP+PcD436O/kLgceBnZ7jOcgqAeddM/y+sAiaBw93280td00Vq/bvAy/Q/KfJI1/YocF+3fxX9T4AcB74F/NDAuY905x3jMv2k0yhrBv458H8Hfq+HgQ8vdT2L+TseuMayCwD/UxCS1Cg/BSRJjTIAJKlRBoAkNcoAkKRGGQCS1CgDQJIaZQBIUqP+H7vYbdimuhf1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X_train , X_test , Y_train, Y_test, W_train, W_test = tts(X,Y,g,test_size=0.25,random_state=0)\n",
    "sklearn_model = LogisticRegression(max_iter=500000, class_weight='balanced')#solver='saga')\n",
    "sklearn_model = tree.DecisionTreeClassifier(max_depth=44, min_samples_leaf=2, random_state = 0)\n",
    "sklearn_model = RandomForestClassifier(n_estimators=50,n_jobs=-1,max_features=\"sqrt\",random_state=0)\n",
    "sklearn_model = svm.SVC(class_weight='balanced')#solver='saga')\n",
    "\n",
    "n_comps = []\n",
    "Qw=[]\n",
    "Qx=[]\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5))\n",
    "print(Q.mean(),Q)\n",
    "plt.plot(0,Q.mean(),marker=\"o\",color=\"blue\")\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='balanced_accuracy')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='recall_weighted')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='precision_weighted')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='f1_weighted')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='f1')\n",
    "print(Q.mean(),Q)\n",
    "plt.plot(0,Q.mean(),marker=\"o\",color=\"orange\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b88342e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9597448618001417 [0.96172927 0.9588944  0.95747697 0.95712261 0.96350106]\n",
      "0.9595610226533188 [0.96178522 0.95859807 0.95735907 0.95669654 0.96336622]\n",
      "0.9597448618001417 [0.96172927 0.9588944  0.95747697 0.95712261 0.96350106]\n",
      "0.960589886286859 [0.96246186 0.95994013 0.95824171 0.95800213 0.96430359]\n",
      "0.959718539334682 [0.96171635 0.95885855 0.95745403 0.95708354 0.96348023]\n",
      "0.9609066392925805 [0.96236934 0.96032832 0.95856354 0.9588015  0.96447051]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1bef1ac8d30>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVXUlEQVR4nO3df4ydVX7f8fcHRmaLlqxZ20sTDLbJQot3RZbV4KaNkKk37cJK4YcXVVCiwpaGPxqqtIQIkPPH1ltEWJMSVUuVmoKKIzZAKGUdpRVLDFarKssyrDEsds1OTACbtHgJXpW6WWrn2z/uGbjMM6yvZ+54POb9kkb3ec45z5nzZaT53Oc5d0yqCkmS+p0w1wuQJB17DAdJUofhIEnqMBwkSR2GgySpY2SuFzAMixcvruXLl8/1MiRpXnnuued+WFVLpuo7LsJh+fLljI2NzfUyJGleSfLqh/X5WEmS1GE4SJI6DAdJUofhIEnqMBwkSR2GgzQbXnkQHl8O3zyh9/rKg3O9IumIHBcfZZWOKa88CN+9AQ4d6J0feLV3DrDimrlbl3QEvHOQhm37uveDYcKhA712aZ4wHKRhO/DakbVLxyDDQRq2k888snbpGGQ4SMP2c7fDiSd/sO3Ek3vt0jxhOEjDtuIaWLURTl4GpPe6aqOb0ZpX/LSSNBtWXGMYaF7zzkGS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgYKhyQXJ9mVZDzJrVP0L0uyJckLSbYmWdrXd2aSbyfZmWRHkuWtfUWSZ9qcDydZ0Df+6STb2nxfGlKtkqQBHTYckpwI3ANcAqwErk6yctKwu4BNVXUesB64o69vE7Chqs4FVgFvtvY7gbur6tPA28D1rf03gUeq6nzgKuDfTqcwSdL0DXLnsAoYr6rdVfUu8BBw2aQxK4Gn2vHTE/0tREaq6kmAqnqnqg4kCbAGeLRd8wBweTsu4Kfa8SeAN460KEnSzAwSDqcDr/ed72lt/bYDa9vxFcApSRYB5wD7kzzWHhNtaHcii4D9VXVwijm/Cvxykj3Afwb+2VSLSnJDkrEkY/v27RugDEnSoIa1IX0zsDrJNmA1sBc4RO+fBL+w9V8AnAVcd5i5rgb+Q1UtBb4E/F6SzjqramNVjVbV6JIlS4ZUhiQJBguHvcAZfedLW9t7quqNqlrb9gnWtbb99O4Inm+PpA4CjwOfB94CFiYZmWLO64FH2hx/AnwMWHzElUmSpm2QcHgWOLt9umgBvU3izf0Dkizue3d/G3B/37ULk0y8tV8D7Kiqorc3cWVrvxb4Vjt+DfhCm/dceuHgcyNJOooOGw7tHf+NwBPATnqfJHopyfokl7ZhFwG7krwMnAbc3q49RO+R0pYkLwIB7m3X3ALclGSc3h7Efa3914FfSbId+H3guhYmkqSjJMfD793R0dEaGxub62VI0ryS5LmqGp2qz7+QliR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1DBQOSS5OsivJeJJbp+hflmRLkheSbE2ytK/vzCTfTrIzyY4ky1v7iiTPtDkfTrKg75p/0Ma+lOSbQ6hTknQEDhsOSU4E7gEuAVYCVydZOWnYXcCmqjoPWA/c0de3CdhQVecCq4A3W/udwN1V9WngbeD69v3OBm4DfqGqPgP88+mVJkmarkHuHFYB41W1u6reBR4CLps0ZiXwVDt+eqK/hchIVT0JUFXvVNWBJAHWAI+2ax4ALm/HvwLcU1Vvt2smwkSSdJQMEg6nA6/3ne9pbf22A2vb8RXAKUkWAecA+5M8lmRbkg3tTmQRsL+qDk4x5znAOUn+e5LvJLl4qkUluSHJWJKxffv2DVCGJGlQw9qQvhlYnWQbsBrYCxwCRoALW/8FwFnAdYeZawQ4G7gIuBq4N8nCyYOqamNVjVbV6JIlS4ZThSQJGCwc9gJn9J0vbW3vqao3qmptVZ0PrGtt++ndETzfHkkdBB4HPg+8BSxMMjLFnHuAzVX1/6rqFeBlemEhSTpKBgmHZ4Gz26eLFgBXAZv7ByRZnGRirtuA+/uuXZhk4q39GmBHVRW9vYkrW/u1wLfa8eP07hpIspjeY6bdR1aWJGkmDhsO7R3/jcATwE7gkap6Kcn6JJe2YRcBu5K8DJwG3N6uPUTvkdKWJC8CAe5t19wC3JRknN4exH2t/QngrSQ76AXIb1TVWzOuVJI0sPTexM9vo6OjNTY2NtfLkKR5JclzVTU6VZ9/IS1J6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6hgoHJJcnGRXkvEkt07RvyzJliQvJNmaZGlf35lJvp1kZ5IdSZa39hVJnmlzPpxkwaQ5v5ykkozOsEZJ0hE6bDgkORG4B7gEWAlcnWTlpGF3AZuq6jxgPXBHX98mYENVnQusAt5s7XcCd1fVp4G3gev7vucpwK8Bz0ynKEnSzAxy57AKGK+q3VX1LvAQcNmkMSuBp9rx0xP9LURGqupJgKp6p6oOJAmwBni0XfMAcHnffF+jFx5/ecQVSZJmbJBwOB14ve98T2vrtx1Y246vAE5Jsgg4B9if5LEk25JsaHcii4D9VXVw8pxJPg+cUVV/9JMWleSGJGNJxvbt2zdAGZKkQQ1rQ/pmYHWSbcBqYC9wCBgBLmz9FwBnAdd92CRJTgD+NfDrh/uGVbWxqkaranTJkiUzLkCS9L5BwmEvcEbf+dLW9p6qeqOq1lbV+cC61raf3h3B8+2R1EHgceDzwFvAwiQjk+Y8BfgssDXJnwE/D2x2U1qSjq5BwuFZ4Oz26aIFwFXA5v4BSRa3d/0AtwH39127MMnEW/s1wI6qKnp7E1e29muBb1XVj6pqcVUtr6rlwHeAS6tqbJr1SZKm4bDh0N7x3wg8AewEHqmql5KsT3JpG3YRsCvJy8BpwO3t2kP0HiltSfIiEODeds0twE1JxuntQdw3tKokSTOS3pv4+W10dLTGxry5kKQjkeS5qprysb1/IS1J6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6hgoHJJcnGRXkvEkt07RvyzJliQvJNmaZGlf35lJvp1kZ5IdSZa39hVJnmlzPpxkQWu/qY17oc25bEi1SpIGdNhwSHIicA9wCbASuDrJyknD7gI2VdV5wHrgjr6+TcCGqjoXWAW82drvBO6uqk8DbwPXt/ZtwGib61Hg69MpTJI0fYPcOawCxqtqd1W9CzwEXDZpzErgqXb89ER/C5GRqnoSoKreqaoDSQKsoffLH+AB4PI25umqOtDavwO8dxciSTo6BgmH04HX+873tLZ+24G17fgK4JQki4BzgP1JHkuyLcmGdieyCNhfVQd/wpzQu5v4L1MtKskNScaSjO3bt2+AMiRJgxrWhvTNwOok24DVwF7gEDACXNj6LwDOAq4bZMIkvwyMAhum6q+qjVU1WlWjS5YsmXEBkqT3jQwwZi9wRt/50tb2nqp6g3bnkOTjwJeran+SPcDzVbW79T0O/DxwP7AwyUi7e/jAnEl+EVgHrK6qH0+zNknSNA1y5/AscHb7dNEC4Cpgc/+AJIuTTMx1G71f/hPXLkwy8dZ+DbCjqore3sSVrf1a4FttrvOBfwdcWlUTm9eSpKPosOHQ3tnfCDwB7AQeqaqXkqxPcmkbdhGwK8nLwGnA7e3aQ/QeKW1J8iIQ4N52zS3ATUnG6e1B3NfaNwAfB/4gyfNJPhBEkqTZl96b+PltdHS0xsbG5noZkjSvJHmuqkan6vMvpCVJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHQOFQ5KLk+xKMp7k1in6lyXZkuSFJFuTLO3rOzPJt5PsTLIjyfLWviLJM23Oh5MsaO0ntfPx1r98OKVKkgZ12HBIciJwD3AJsBK4OsnKScPuAjZV1XnAeuCOvr5NwIaqOhdYBbzZ2u8E7q6qTwNvA9e39uuBt1v73W2cJOkoGuTOYRUwXlW7q+pd4CHgskljVgJPteOnJ/pbiIxU1ZMAVfVOVR1IEmAN8Gi75gHg8nZ8WTun9X+hjZckHSWDhMPpwOt953taW7/twNp2fAVwSpJFwDnA/iSPJdmWZEO7E1kE7K+qg1PM+d73a/0/auMlSUfJsDakbwZWJ9kGrAb2AoeAEeDC1n8BcBZw3TC+YZIbkowlGdu3b98wppQkNYOEw17gjL7zpa3tPVX1RlWtrarzgXWtbT+9O4Ln2yOpg8DjwOeBt4CFSUammPO979f6P9HGf0BVbayq0aoaXbJkyQBlSJIGNUg4PAuc3T5dtAC4CtjcPyDJ4iQTc90G3N937cIkE7+91wA7qqro7U1c2dqvBb7Vjje3c1r/U228JOkoOWw4tHf8NwJPADuBR6rqpSTrk1zahl0E7EryMnAacHu79hC9R0pbkrwIBLi3XXMLcFOScXp7Cve19vuARa39JqDz0VlJ0uzK8fCmfHR0tMbGxuZ6GZI0ryR5rqpGp+rzL6QlSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR0DhUOSi5PsSjKe5NYp+pcl2ZLkhSRbkyzt6zuU5Pn2tbmvfU2S7yX5fpIHkoy09k8k+cMk25O8lOQrwyhUkjS4w4ZDkhOBe4BLgJXA1UlWThp2F7Cpqs4D1gN39PX936r6XPu6tM15AvAAcFVVfRZ4Fbi2jf9VYEdV/RxwEfDbSRZMt0BJ0pEb5M5hFTBeVbur6l3gIeCySWNWAk+146en6J9sEfBuVb3czp8EvtyOCzglSYCPA38BHBxgnZKkIRkkHE4HXu8739Pa+m0H1rbjK+j9cl/Uzj+WZCzJd5Jc3tp+CIwkGW3nVwJntONvAOcCbwAvAr9WVX81YD2SpCEY1ob0zcDqJNuA1cBe4FDrW1ZVo8A/BH4nyc9WVQFXAXcn+S7wv/vGfxF4HvgZ4HPAN5L81ORvmOSGFjpj+/btG1IZkiQYLBz28v67eoClre09VfVGVa2tqvOBda1tf3vd2153A1uB89v5n1TVhVW1CvivwMQjpq8Aj1XPOPAK8DcnL6qqNlbVaFWNLlmyZMByJUmDGCQcngXOTrKibQxfBWzuH5BkcdtkBrgNuL+1n5rkpIkxwC8AO9r5p9rrScAtwO+2618DvtD6TgP+BrB7ugVKko7cYcOhqg4CNwJPADuBR6rqpSTrk1zahl0E7EryMnAacHtrPxcYS7Kd3kb1b1XVjtb3G0l2Ai8Af1hVExvaXwP+TpIXgS3ALVX1w5kWKkkaXHqP/+e30dHRGhsbm+tlSNK8kuS5tifc4V9IS5I6DAdpFjz4ICxfDiec0Ht98MG5XpF0ZEbmegHS8ebBB+GGG+DAgd75q6/2zgGuuWbu1iUdCe8cpCFbt+79YJhw4ECvXZovDAdpyF577cjapWOR4SAN2ZlnHlm7dCwyHKQhu/12OPnkD7adfHKvXZovDAdpyK65BjZuhGXLIOm9btzoZrTmFz+tJM2Ca64xDDS/eecgSeowHCRJHYaDJKnDcJAkdRgOkqSO4+Kf7E6yD3h1rtcxDYvp/f+0P0o+ajV/1OoFa55PllXVlP8rzeMiHOarJGMf9m+pH68+ajV/1OoFaz5e+FhJktRhOEiSOgyHubVxrhcwBz5qNX/U6gVrPi645yBJ6vDOQZLUYThIkjoMh1mW5JNJnkzyg/Z66oeMu7aN+UGSa6fo35zk+7O/4pmZSb1JTk7yR0n+R5KXkvzW0V39kUlycZJdScaT3DpF/0lJHm79zyRZ3td3W2vfleSLR3XhMzDdmpP8vSTPJXmxva456oufppn8nFv/mUneSXLzUVv0MFSVX7P4BXwduLUd3wrcOcWYTwK72+up7fjUvv61wDeB7891PbNZL3Ay8HfbmAXAfwMumeuaPqTOE4E/Bc5qa90OrJw05p8Cv9uOrwIebscr2/iTgBVtnhPnuqZZrvl84Gfa8WeBvXNdz2zX3Nf/KPAHwM1zXc+RfHnnMPsuAx5oxw8Al08x5ovAk1X1F1X1NvAkcDFAko8DNwH/avaXOhTTrreqDlTV0wBV9S7wPWDp7C95WlYB41W1u631IXq19+v/b/Eo8IUkae0PVdWPq+oVYLzNd6ybds1Vta2q3mjtLwF/LclJR2XVMzOTnzNJLgdeoVfzvGI4zL7TqurP2/H/BE6bYszpwOt953taG8DXgN8GDszaCodrpvUCkGQh8EvAlllY4zActob+MVV1EPgRsGjAa49FM6m535eB71XVj2dpncM07ZrbG7tbgH95FNY5dP6f4IYgyR8Df32KrnX9J1VVSQb+7HCSzwE/W1X/YvJzzLk0W/X2zT8C/D7wb6pq9/RWqWNRks8AdwJ/f67XchR8Fbi7qt5pNxLziuEwBFX1ix/Wl+R/JfnpqvrzJD8NvDnFsL3ARX3nS4GtwN8GRpP8Gb2f1aeSbK2qi5hDs1jvhI3AD6rqd2a+2lmzFzij73xpa5tqzJ4WeJ8A3hrw2mPRTGomyVLgPwH/qKr+dPaXOxQzqflvAVcm+TqwEPirJH9ZVd+Y9VUPw1xvehzvX8AGPrhB+/UpxnyS3nPJU9vXK8AnJ41ZzvzYkJ5RvfT2Vv4jcMJc13KYOkfobaSv4P2Nys9MGvOrfHCj8pF2/Bk+uCG9m/mxIT2Tmhe28Wvnuo6jVfOkMV9lnm1Iz/kCjvcves9btwA/AP6475fgKPDv+8b9Y3obk+PAV6aYZ76Ew7TrpfeurICdwPPt65/MdU0/odYvAS/T+zTLuta2Hri0HX+M3qdUxoHvAmf1XbuuXbeLY/QTWcOsGfhN4P/0/VyfBz411/XM9s+5b455Fw7+8xmSpA4/rSRJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjr+PyBWreZMmOajAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "X_train , X_test , Y_train, Y_test, W_train, W_test = tts(X,Y,g,test_size=0.25,random_state=0)\n",
    "sklearn_model = LogisticRegression(max_iter=500000, class_weight='balanced')#solver='saga')\n",
    "sklearn_model = tree.DecisionTreeClassifier(max_depth=44, min_samples_leaf=2, random_state = 0)\n",
    "sklearn_model = RandomForestClassifier(n_estimators=50,n_jobs=-1,max_features=\"sqrt\",random_state=0)\n",
    "sklearn_model = svm.SVC(class_weight='balanced')#solver='saga')\n",
    "\n",
    "n_comps = []\n",
    "Qw=[]\n",
    "Qx=[]\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5))\n",
    "print(Q.mean(),Q)\n",
    "plt.plot(0,Q.mean(),marker=\"o\",color=\"blue\")\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='balanced_accuracy')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='recall_weighted')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='precision_weighted')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='f1_weighted')\n",
    "print(Q.mean(),Q)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='f1')\n",
    "print(Q.mean(),Q)\n",
    "plt.plot(0,Q.mean(),marker=\"o\",color=\"orange\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94e70f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train , X_test , Y_train, Y_test, W_train, W_test = tts(X,Y,g,test_size=0.05,random_state=5)\n",
    "sklearn_model = LogisticRegression(max_iter=500000, class_weight='balanced')#solver='saga')\n",
    "sklearn_model.fit(X_train,Y_train)\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5))# scoring='balanced_accuracy')\n",
    "print(Q.mean(),Q,np.std(Q),sklearn_model.score(X_test,Y_test),sklearn_model.score(X_train,Y_train))\n",
    "plt.plot(0,Q.mean(),marker=\"o\",color=\"blue\")\n",
    "Q=cross_val_score(sklearn_model, X,Y, groups = g, cv = GroupKFold(5), scoring='f1')\n",
    "print(Q.mean(),Q,np.std(Q),sklearn_model.score(X_test,Y_test),sklearn_model.score(X_train,Y_train))\n",
    "plt.plot(0,Q.mean(),marker=\"o\",color=\"orange\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "334c8d45",
   "metadata": {},
   "source": [
    "#### Feature visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b261b147",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "coefs = []\n",
    "ba = []\n",
    "for i in tq(range(40)):\n",
    "    X_train , X_test , Y_train, Y_test, W_train, W_test = tts(X,Y,g,test_size=0.05,random_state=i)\n",
    "    sklearn_model = LogisticRegression(max_iter=500000, class_weight='balanced')#solver='saga')\n",
    "    cross_val_score(sklearn_model, X_train,Y_train, groups = g, cv = GroupKFold(5))\n",
    "    coefs.append(sklearn_model.coef_[0])\n",
    "    ba.append([sklearn_model.score(X_test,Y_test),sklearn_model.score(X_train,Y_train)])\n",
    "\n",
    "#print(pd.DataFrame(coefs,columns=f).T.to_string())\n",
    "coefsa = np.array(coefs)\n",
    "plt.figure(figsize=(5,30))\n",
    "plt.grid()\n",
    "print(pd.DataFrame([np.mean(coefsa,0),np.std(coefsa,0)],columns=f).T.to_string())\n",
    "plt.barh(y=np.flip(f),width=np.mean(np.flip(coefsa),0),xerr=np.std(np.flip(coefsa),0))\n",
    "\n",
    "#coefsa = np.flip(coefsa)\n",
    "#print(pd.DataFrame([np.mean(coefsa,0),np.std(coefsa,0)],columns=f).T.to_string())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0112f58b",
   "metadata": {},
   "source": [
    "## Post-Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a038e2f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "coefs = []\n",
    "ba = []\n",
    "for i in tq(range(40)):\n",
    "    X_train , X_test , Y_train, Y_test, g_train, g_test = tts(X,Y,g,test_size=0.2,random_state=i)\n",
    "    sklearn_model = LogisticRegression(max_iter=500000, class_weight='balanced')#solver='saga')\n",
    "    sklearn_model.fit(X_train,Y_train)#, groups = g_train, cv = GroupKFold(5))\n",
    "    #coefs.append(sklearn_model.coef_[0])\n",
    "    ba.append([sklearn_model.score(X_test,Y_test),sklearn_model.score(X_train,Y_train)])\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.grid()\n",
    "from matplotlib import pyplot as plt\n",
    "from itertools import cycle\n",
    "import numpy as np\n",
    "colors = cycle(['#368195', '#4198AF', '#91C3D5', '#BBD7E3'])\n",
    "bars=plt.bar(x=[i for i in range(len(ba[0]))],height=np.mean(ba,0),yerr=np.std(ba,0),color=[next(colors) for _ in range(len(ba[0]))])\n",
    "plt.yticks(np.arange(0,1,0.1))\n",
    "plt.ylabel('accuracy')\n",
    "plt.legend(handles=bars, labels=['test CREST', 'train CREST', 'test MMFF', 'train MMFF'],loc = 'center right',bbox_to_anchor=(1.2, 0.5))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6148fbcc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "coefs = []\n",
    "ba = []\n",
    "for i in tq(range(40)):\n",
    "    X_train , X_test , Y_train, Y_test, W_train, W_test = tts(X,Y,g,test_size=0.05,random_state=i)\n",
    "    sklearn_model = LogisticRegression(max_iter=500000, class_weight='balanced')#solver='saga')\n",
    "    cross_val_score(sklearn_model, X_train,Y_train, groups = g, cv = GroupKFold(5))\n",
    "    coefs.append(sklearn_model.coef_[0])\n",
    "    ba.append([sklearn_model.score(X_test,Y_test),sklearn_model.score(X_train,Y_train)])\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.grid()\n",
    "\n",
    "colors = cycle(['#368195', '#4198AF', '#91C3D5', '#BBD7E3'])\n",
    "bars=plt.bar(x=[i for i in range(len(ba[0]))],height=np.mean(ba,0),yerr=np.std(ba,0),color=[next(colors) for _ in range(len(ba[0]))])\n",
    "plt.yticks(np.arange(0,1,0.1))\n",
    "plt.ylabel('accuracy')\n",
    "plt.legend(handles=bars, labels=['test CREST', 'train CREST', 'test MMFF', 'train MMFF'],loc = 'center right',bbox_to_anchor=(1.2, 0.5))\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "my-rdkit-env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "7f57446601f641341b8bf34e3127a7399b7cc86b59253d3e5183adc8a9149cb7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}