same_model_different_test_set_atoms-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing for individual systems.\n",
    "In the following block of codes we start with csv file of a particular system ('Ag2O') that contains the \n",
    "information of compound name. We then shuffle and split the selected rows into 80/20 train/test dataset. We will \n",
    "get the parameters for best performing model using different kernels (we used laplacian kernel in our case) in \n",
    "Kernel Ridge Regression (KRR). We then shuffle and get new train/test data using 80/20 split and use that best \n",
    "performing model to test the new train/test data. We repeat that step 20 times. We first start by importing the \n",
    "required modules. If you don't have some or all of the modules below try \"pip install pandas\" and so on for \n",
    "all missing modules in the terminal.\n",
    "We end up with a train and test-set MAE corresponding to various test-sets across the same system. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np \n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.kernel_ridge import KernelRidge\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import make_scorer\n",
    "#\n",
    "import tarfile\n",
    "import re\n",
    "import sys\n",
    "import random \n",
    "import operator \n",
    "import os\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mae(y_true, y_pred):\n",
    "    return np.mean(abs(y_true-y_pred))\n",
    "\n",
    "def root_mean_squared_error(y_true, y_pred):\n",
    "    return np.sqrt(np.mean((y_true-y_pred)**2))\n",
    "\n",
    "def string_number_separator(string):\n",
    "    import re\n",
    "    temp = re.compile(\"([a-zA-Z]+)([0-9]+)\")\n",
    "    res = temp.match(string).groups()\n",
    "    return res\n",
    "\n",
    "def list_of_atoms_and_numbers(string):\n",
    "    atoms=[]\n",
    "    atoms_number=[]\n",
    "    str1=''\n",
    "    while len(str1)<len(string):\n",
    "        string1=str(string)[len(str1):]\n",
    "        x,y=string_number_separator(string1)\n",
    "        str1=str1+str(x)+str(y)\n",
    "        atoms.append(str(x))\n",
    "        atoms_number.append(int(y))\n",
    "    return atoms, atoms_number\n",
    "\n",
    "def input_column_creater(atom_list,prop):\n",
    "    orbs=['1s','2s','2p','3s','3p','3d','4s','4p','4d','5s','5p','5d','6s','6p']\n",
    "    orbs_generated=[]\n",
    "    for elements in atom_list:\n",
    "        for i in range (len(orbs)):\n",
    "            new_str=str(orbs[i]+str(elements))\n",
    "            orbs_generated.append(new_str)\n",
    "    for j in range (len(prop)):\n",
    "        orbs_generated.append(str(prop[j]))\n",
    "    return orbs_generated\n",
    "\n",
    "def train_model(x_train, y_train, x_test,y_test, nfold, nthread, in_alpha, in_gamma, kernel=None, rseed=None):\n",
    "\n",
    "    if rseed:\n",
    "        random.seed(rseed)\n",
    "        np.random.seed(rseed)\n",
    "\n",
    "    neg_root_mean_squared_error = make_scorer(root_mean_squared_error, greater_is_better=False)\n",
    "    clf = GridSearchCV(KernelRidge(kernel=kernel), cv=nfold, n_jobs=nthread, verbose=1, scoring=neg_root_mean_squared_error, param_grid={\"alpha\":in_alpha, \"gamma\": in_gamma})\n",
    "    clf.fit(x_train, y_train)\n",
    "    y_train_pred = clf.predict(x_train)\n",
    "    y_test_pred = clf.predict(x_test)\n",
    "    train_err = mae(y_train_pred, y_train)\n",
    "    test_err = mae(y_test_pred, y_test)\n",
    "    return clf,train_err,test_err\n",
    "\n",
    "def predict_value(clf,x_test,y_test):\n",
    "    y_test_pred=clf.predict(x_test)\n",
    "    test_err=mae(y_test_pred,y_test)\n",
    "    return test_err\n",
    "\n",
    "def modified_pettifor(in_atom):\n",
    "\n",
    "    mod_scale = {\n",
    "                \"He\": 1,\n",
    "                \"Ne\": 2,\n",
    "                \"Ar\": 3,\n",
    "                \"Kr\": 4,\n",
    "                \"Xe\": 5,\n",
    "                \"Rn\": 6,\n",
    "                \"Fr\": 7,\n",
    "                \"Cs\": 8,\n",
    "                \"Rb\": 9,\n",
    "                \"K\": 10,\n",
    "                \"Na\": 11,\n",
    "                \"Li\": 12,\n",
    "                \"Ra\": 13,\n",
    "                \"Ba\": 14,\n",
    "                \"Sr\": 15,\n",
    "                \"Ca\": 16,\n",
    "                \"Eu\": 17,\n",
    "                \"Yb\": 18,\n",
    "                \"Lu\": 19,\n",
    "                \"Tm\": 20,\n",
    "                \"Y\": 21,\n",
    "                \"Er\": 22,\n",
    "                \"Ho\": 23,\n",
    "                \"Dy\": 24,\n",
    "                \"Tb\": 25,\n",
    "                \"Gd\": 26,\n",
    "                \"Sm\": 27,\n",
    "                \"Pm\": 28,\n",
    "                \"Nd\": 29,\n",
    "                \"Pr\": 30,\n",
    "                \"Ce\": 31,\n",
    "                \"La\": 32,\n",
    "                \"Ac\": 33,\n",
    "                \"Th\": 34,\n",
    "                \"Pa\": 35,\n",
    "                \"U\": 36,\n",
    "                \"Np\": 37,\n",
    "                \"Pu\": 38,\n",
    "                \"Am\": 39,\n",
    "                \"Cm\": 40,\n",
    "                \"Bk\": 41,\n",
    "                \"Cf\": 42,\n",
    "                \"Es\": 43,\n",
    "                \"Fm\": 44,\n",
    "                \"Md\": 45,\n",
    "                \"No\": 46,\n",
    "                \"Lr\": 47,\n",
    "                \"Sc\": 48,\n",
    "                \"Zr\": 49,\n",
    "                \"Hf\": 50,\n",
    "                \"Ti\": 51,\n",
    "                \"Ta\": 52,\n",
    "                \"Nb\": 53,\n",
    "                \"V\": 54,\n",
    "                \"Cr\": 55,\n",
    "                \"Mo\": 56,\n",
    "                \"W\": 57,\n",
    "                \"Re\": 58,\n",
    "                \"Tc\": 59,\n",
    "                \"Os\": 60,\n",
    "                \"Ru\": 61,\n",
    "                \"Ir\": 62,\n",
    "                \"Rh\": 63,\n",
    "                \"Pt\": 64,\n",
    "                \"Pd\": 65,\n",
    "                \"Au\": 66,\n",
    "                \"Ag\": 67,\n",
    "                \"Cu\": 68,\n",
    "                \"Ni\": 69,\n",
    "                \"Co\": 70,\n",
    "                \"Fe\": 71,\n",
    "                \"Mn\": 72,\n",
    "                \"Mg\": 73,\n",
    "                \"Zn\": 74,\n",
    "                \"Cd\": 75,\n",
    "                \"Hg\": 76,\n",
    "                \"Be\": 77,\n",
    "                \"Al\": 78,\n",
    "                \"Ga\": 79,\n",
    "                \"In\": 80,\n",
    "                \"Tl\": 81,\n",
    "                \"Pb\": 82,\n",
    "                \"Sn\": 83,\n",
    "                \"Ge\": 84,\n",
    "                \"Si\": 85,\n",
    "                \"B\": 86,\n",
    "                \"C\": 87,\n",
    "                \"N\": 88,\n",
    "                \"P\": 89,\n",
    "                \"As\": 90,\n",
    "                \"Sb\": 91,\n",
    "                \"Bi\": 92,\n",
    "                \"Po\": 93,\n",
    "                \"Te\": 94,\n",
    "                \"Se\": 95,\n",
    "                \"S\": 96,\n",
    "                \"O\": 97,\n",
    "                \"At\": 98,\n",
    "                \"I\": 99,\n",
    "                \"Br\": 100,\n",
    "                \"Cl\": 101,\n",
    "                \"F\": 102,\n",
    "                \"H\": 103\n",
    "            }\n",
    "\n",
    "    return mod_scale[in_atom]/float(max(mod_scale.values()))\n",
    "\n",
    "def compound_pettifor(compound):\n",
    "    cpd=compound.split('-')[0]\n",
    "    in_atom,num_atom=list_of_atoms_and_numbers(cpd)\n",
    "    if len(in_atom)==1:\n",
    "        norm_pettifor=modified_pettifor(in_atom[0])\n",
    "        return norm_pettifor\n",
    "    else:\n",
    "        norm_pettifor=0\n",
    "        for i in range (len(in_atom)):\n",
    "            tot=sum(num_atom)\n",
    "            norm_pettifor+=(num_atom[i]/tot)*modified_pettifor(in_atom[i])\n",
    "        return norm_pettifor\n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def best_fit(csv_path):\n",
    "    df = pd.read_csv(csv_path)\n",
    "    df['pet'] = [ compound_pettifor(compound) for compound in df['compound'] ]\n",
    "    tmp_feats = ['pet', '1s', '2s', '2p', '3s', '3p', '3d', '4s', '4p', '4d', '5s', '5p', '5d', '6s', '6p']\n",
    "    keep_cols = tmp_feats + ['EF_PBE']\n",
    "    # print(keep_cols)\n",
    "    # keep_cols = columns\n",
    "    y_feat = 'eigshift(eV)'\n",
    "    #y_feat='EF_HSE'\n",
    "    x_data = df[keep_cols].to_numpy()\n",
    "    y_data = df[y_feat].to_numpy()\n",
    "\n",
    "    print(np.shape(x_data)) #x_data)\n",
    "    print(np.shape(y_data)) #y_data)\n",
    "\n",
    "\n",
    "    test_size = 0.20 \n",
    "    nfold = 5\n",
    "    nthread = 2\n",
    "    ker = 'laplacian'\n",
    "    #ker = 'rbf'\n",
    "    alphas = [0.0000001, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0]\n",
    "    gammas = [0.0001, 0.001, 0.01, 0.1, 1.0 , 10.0 , 100.0]\n",
    "    number = random.randint(1, 1000)\n",
    "    print(\"random seed\", number)\n",
    "    rseed = random.seed(number)\n",
    "    #rseed=300\n",
    "    test_err=[]\n",
    "    x_train, x_test, y_train, y_test = train_test_split(x_data, y_data, test_size=test_size, shuffle=True)\n",
    "    clf,train_error,test_error= train_model(x_train, y_train, x_test, y_test, nfold, nthread, in_alpha=alphas, in_gamma=gammas, kernel=ker, rseed=rseed)\n",
    "    test_err.append(test_error)\n",
    "    print(train_error)\n",
    "    for i in range (25):\n",
    "        x_train, x_test, y_train, y_test = train_test_split(x_data, y_data, test_size=test_size, shuffle=True)\n",
    "        clf.fit(x_train,y_train)\n",
    "        test_mae=predict_value(clf,x_test,y_test)\n",
    "        test_err.append(test_mae)\n",
    "    return test_err,clf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cwd=os.getcwd()\n",
    "atom_path=os.path.join(cwd,'atom_sum_compound_csv_files/')\n",
    "concat_path=os.path.join(cwd,'special_concatenations_compounds_atom_sum/')\n",
    "listdir_cwd=os.listdir(cwd)\n",
    "#print(listdir)\n",
    "check=['atom_sum_compound_csv_files','special_concatenations_compounds_atom_sum']\n",
    "for entries in check:\n",
    "    if entries not in listdir_cwd:\n",
    "        file=tarfile.open(str(entries)+'.tar.gz')\n",
    "        file.extractall()\n",
    "        file.close\n",
    "    else:\n",
    "        print ('%s %s'%(entries,'exists'))\n",
    "\n",
    "test_err,clf=best_fit(atom_path+'Ag2O1.csv')\n",
    "test_error=np.array(test_err)\n",
    "print(clf.best_params_)\n",
    "print(np.min(test_error),np.max(test_error),np.mean(test_error),np.std(test_error))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}