krr_lap_concatenated_mod_pettifor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Best bit parameters and concatenated file\n",
    "\n",
    " The purpose of the following notebook is to perform grid search to find the best two parameters namely,\n",
    " alpha and gamma that are to be properly tuned in Kernel Ridge Regression using Laplacian Kernel. In addition \n",
    " we choose 10 percent (1 percent here ) of the whole randomly selected concatenated data and do 80/20 train/test \n",
    " split. As usual we start by importing required modules. Try installing using pip (eg. pip install scikit-learn ) \n",
    " if you don't have some or all of the following modules. Since most of the functions used here are used on other \n",
    " notebooks herewith, I will not describe these here. We get a pdf file plotting train/test error and the best\n",
    " fit parameter. Here we can choose kernel = 'laplacian' or 'rbf', etc.,\n"
   ]
  },
  {
   "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_test, y_test, x_train, y_train, nfold, nthread, in_alpha=np.logspace(-15, 5, 21, base=2), in_gamma=np.logspace(-15, 3, 19, base=2), 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",
    "\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",
    "\n",
    "    clf.fit(x_train, y_train)\n",
    "\n",
    "    y_train_pred = clf.predict(x_train)\n",
    "    y_test_pred = clf.predict(x_test)\n",
    "\n",
    "    return clf, clf.best_params_, y_train, y_test, y_train_pred, y_test_pred\n",
    "\n",
    "def scatter_comp(x1, y1, x2, y2, error_unit, xlabel, ylabel, plot_name):\n",
    "\n",
    "    import os, sys\n",
    "    import matplotlib as mpl\n",
    "    import matplotlib.pyplot as plt\n",
    "    import pandas as pd\n",
    "    import numpy as np\n",
    "\n",
    "    from sklearn.metrics import r2_score\n",
    "\n",
    "    axis_width = 2\n",
    "    mpl.rcParams['mathtext.default'] = 'regular'\n",
    "    mpl.rcParams['axes.linewidth'] = axis_width\n",
    "    mpl.rcParams['axes.unicode_minus']=False\n",
    "\n",
    "    fontsize = 24\n",
    "    alpha = 0.2\n",
    "    lw45 = 1.5\n",
    "    axis_width=1.5\n",
    "    figsize = (9.2,7)\n",
    "    label_font = 24\n",
    "    tick_font = 20\n",
    "    leg_font = 20\n",
    "    tick_len = 6\n",
    "    test_marker = 'o'  #test_mark = '^'\n",
    "    train_marker = 'o'\n",
    "    train_color = 'gray'\n",
    "    test_color = 'blue' #'red'\n",
    "    pt_alpha = 0.3\n",
    "    pt_lw = 0\n",
    "    pt_s = 80\n",
    "    htp=0.1\n",
    "    diag_color = 'black'\n",
    "    diag_lw = axis_width\n",
    "    diag_ls = '-'\n",
    "    text_font = 20\n",
    "    axis_width = 2\n",
    "    leg_fancy = False\n",
    "    leg_frame = True\n",
    "    leg_alpha = 1\n",
    "    leg_loc = 'upper right'\n",
    "    htp = 0.2\n",
    "    leg_color = 'white'\n",
    "    leg_edge_color = 'black'\n",
    "    leg_shadow = False\n",
    "    leg_lw = axis_width\n",
    "    s = 150\n",
    "    pt_lw = 0\n",
    "    decision_lw = axis_width\n",
    "    decision_ls = '--'\n",
    "    decision_color = 'black'\n",
    "    markersize = 80\n",
    "    alpha = 0.6\n",
    "    lw45 = 1.5\n",
    "    axis_width=1.5\n",
    "\n",
    "    train_err = mae(y1, x1)\n",
    "    test_err = mae(y2, x2)\n",
    "    print((\"MAE TRAIN ERROR\", train_err))\n",
    "    print((\"MAE TEST ERROR\", test_err))\n",
    "\n",
    "    fig = plt.figure(figsize=figsize)    \n",
    "    plt.scatter(x1, y1, color=train_color, marker=train_marker, s=markersize, alpha=alpha, lw=0)\n",
    "    plt.scatter(x2, y2, color=test_color, marker=test_marker, s=markersize, alpha=alpha, lw=0)\n",
    "    plt.legend(['train', 'test'], fontsize=fontsize-4, frameon=False, loc='lower right')\n",
    "\n",
    "    xmax = np.max(x1)\n",
    "    xmin = np.min(x1)\n",
    "    x45 = np.linspace(xmax*.98, xmin*1.02, num=100)\n",
    "    y45 = x45\n",
    "    plt.plot(x45, y45, color='black', lw=lw45, label='__nolegend__')\n",
    "\n",
    "    plt.ylabel(ylabel, fontsize=fontsize)\n",
    "    plt.xlabel(xlabel, fontsize=fontsize)\n",
    "\n",
    "    plt.xticks(fontsize=fontsize)\n",
    "    plt.yticks(fontsize=fontsize)\n",
    "    axes = plt.gca()\n",
    "    axes.set_title(axes.get_title()) #* 2)\n",
    "    axes.set_xlabel(axes.get_xlabel(), size=fontsize) #* 0.5) # fontname=\"Times New Roman\")\n",
    "    axes.set_ylabel(axes.get_ylabel(), size=fontsize)\n",
    "\n",
    "    fig_label = \"RMSE\"\n",
    "\n",
    "    plt.tick_params('both', length = 6, width = axis_width, which = 'major',right=True,top=True)\n",
    "\n",
    "    fig.savefig(plot_name, dpi=600, bbox_inches='tight')\n",
    "    plt.show()\n",
    "    return train_err,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": [
    "## use ker='rbf' if you want to use rbf kernel here. \n",
    "def best_fit(csv_path,ker='laplacian'):\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','EF_PBE','eigshift(eV)']\n",
    "    keep_cols = tmp_feats\n",
    "    x_data_ini=df[keep_cols].to_numpy()\n",
    "    number_of_rows=x_data_ini.shape[0]\n",
    "    random_indices=x_data_ini[np.random.randint(x_data_ini.shape[0], size=int(number_of_rows/100)), :] # 1 percent \n",
    "    x_data =random_indices[:, :-1]\n",
    "    y_data =random_indices[:, -1]\n",
    "\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",
    "\n",
    "    number = random.randint(1, 1000)\n",
    "    print(\"random seed\", number)\n",
    "    rseed = random.seed(number)\n",
    "\n",
    "    x_train, x_test, y_train, y_test = train_test_split(x_data, y_data, test_size=test_size, shuffle=True)\n",
    "\n",
    "    alphas = [0.0000001, 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.0] #np.logspace(-15, 5, 21, base=2) #  # #\n",
    "    gammas = [0.0001, 0.001, 0.01, 0.1, 1.0, 10.0, 100.0] #np.linspace(0.001,0.1, 20) #np.logspace(-15, 5, 21, base=2) \n",
    "\n",
    "    res_dict = { idx: {'pred_test': [], 'pred_train': [], \"best_params\":{'alpha': 0, 'gamma': 0}} for idx in range(len(gammas)) }\n",
    "    all_data_mae = []\n",
    "    for it, g in enumerate(gammas):\n",
    "        clf, clf.best_params_, y_train, y_test, y_train_pred, y_test_pred = train_model(x_test, y_test, x_train, y_train, nfold, nthread, in_alpha=alphas, in_gamma=[g], kernel=ker, rseed=rseed)\n",
    "        print(clf.best_params_)\n",
    "        print(np.shape(y_train_pred), np.shape(y_test_pred))\n",
    "        res_dict[it]['pred_train'].append( y_train_pred )\n",
    "        res_dict[it]['pred_test'].append( y_test_pred )\n",
    "        res_dict[it]['best_params']['alpha'] = clf.best_params_['alpha']\n",
    "        res_dict[it]['best_params']['gamma'] = clf.best_params_['gamma']\n",
    "        print(it, \"RMSE\", root_mean_squared_error(y_test, y_test_pred))\n",
    "        print(it, \"MAE\", np.mean(abs(y_test - y_test_pred)))\n",
    "        diff = abs(y_test - y_test_pred)\n",
    "        all_data_mae.append(np.mean(diff))\n",
    "\n",
    "    min_index, min_value = min(enumerate(all_data_mae), key=operator.itemgetter(1))\n",
    "\n",
    "    tmp_figname =str(csv_path.split('/')[-1])[:-4]+'_'.join([\"train_nsamples=\"+str(len(y_train)), \"ker=\"+ker, \"MAE=\"+str(min_value)])\n",
    "    figname = tmp_figname+\".pdf\"\n",
    "    xlabel = r\"DFT $\\Delta eigen$ (eV)\"\n",
    "    ylabel = r\"ML $\\Delta eigen}$ (eV)\"\n",
    "    error_unit = \"meV\"\n",
    "\n",
    "    print('best parameters', \"gamma:\", res_dict[min_index]['best_params']['gamma'], \"alpha:\", res_dict[min_index]['best_params']['alpha'])\n",
    "    train_err,test_err=scatter_comp(y_train, res_dict[min_index]['pred_train'],  y_test, res_dict[min_index]['pred_test'], error_unit, xlabel, ylabel, figname)\n",
    "    return train_err,test_err,res_dict[min_index]['best_params']['gamma'],res_dict[min_index]['best_params']['alpha']"
   ]
  },
  {
   "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",
    "list_concat_files=os.listdir(concat_path)\n",
    "#print(list_atom_files)\n",
    "name='concatenated_final.csv'\n",
    "train_err,test_err,gamma,alpha=best_fit(concat_path+str(name))\n",
    "list_values=[str(name)[:-4],train_err,test_err,gamma,alpha]\n",
    "#print(cwd)\n",
    "#print(atom_path)\n",
    "#path=cwd+'/Ag2O1.csv'\n",
    "#train_err,test_err,gamma,alpha=best_fit(path)\n",
    "#print('At last')\n",
    "#print(train_err,test_err,gamma,alpha)\n"
   ]
  }
 ],
 "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
}