{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Synopsis\n\n\nThis kernel objective is to check the robustness of Mulliken charges estimated in this kernel : [Predicting Mulliken Charges With ACSF Descriptors](https://www.kaggle.com/borisdee/predicting-mulliken-charges-with-acsf-descriptors)\n\n\n# Context\n\nBoris has built a clever model to estimate Mulliken charges of atoms for the test set in the kernel linked above. These local properties are a significant feature in any serious model for predicting the `scalar_coupling_constant`.\n\nAs we do not have a comparison basis excepted the cross validation of the prediction model itself, I wanted  to check how this data compare against the Open Babel estimated values. The idea here is to check if the correlation of the competiton provided Mulliken charges and Open Babel estimates behave the same when using the Boris set.\n\n# Method\n\nI will use here the same code as in my previous kernel [V7 Estimation of Mulliken Charges with Open Babel](https://www.kaggle.com/asauve/v7-estimation-of-mulliken-charges-with-open-babel). But the data used for correlation will be\n* Boris Mulliken charges estimates for test set\n* Open babel Partial charges estimates for test set -- for each of the 11 method available --\n\nAfter that, the 11 correlation coefficients of the test set will be compared againts the 11 ones obtained with the train set.\n\n"},{"metadata":{},"cell_type":"markdown","source":"# Load Data"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nfrom scipy.stats.stats import pearsonr\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom tqdm import tqdm_notebook as tqdm\nimport seaborn as sns \nimport matplotlib.pyplot as plt\nsns.set()\nimport os","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Load CSV data"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def load_dir_csv(directory, csv_files=None):\n    if csv_files is None:\n        csv_files = sorted( [ f for f in os.listdir(directory) if f.endswith(\".csv\") ])    \n    csv_vars  = [ filename[:-4] for filename in csv_files ]\n    gdict = globals()\n    for filename, var in zip( csv_files, csv_vars ):\n        print(f\"{var:32s} = pd.read_csv({directory}/{filename})\")\n        gdict[var] = pd.read_csv( f\"{directory}/{filename}\" )\n        print(f\"{'nb of cols ':32s} = \" + str(len(gdict[var])))\n        display(gdict[var].head())\n\nload_dir_csv(\"../input/champs-scalar-coupling/\", \n             [\"test.csv\"])\nload_dir_csv(\"../input/boris-mulliken-test\", \n             [\"mulliken_charges_test_set.csv\"])\nload_dir_csv(\"../input/open-babel-atom-charges/\",\n            [\"test_ob_charges.csv\"])\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load Mulliken Charges from Boris test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"# differentiate train and test set\ntest_molecules  = test.molecule_name.unique()\n\nmulliken_test_boris = []\nmulliken_charges_test_set_idx = mulliken_charges_test_set.set_index(['molecule_name'])\n# ensure mulliken charges are in same order as for Open Babel\nfor molecule_name in test_molecules:\n    mc  = mulliken_charges_test_set_idx.loc[molecule_name].sort_index()\n    mulliken_test_boris.extend(mc.mulliken_charge.values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load partial charges from Open Babel set\n\nThese values have been computed in a previous kernel :   [V7 Estimation of Mulliken Charges with Open Babel](https://www.kaggle.com/asauve/v7-estimation-of-mulliken-charges-with-open-babel)"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nob_methods = [ \"eem\", \"mmff94\", \"gasteiger\", \"qeq\", \"qtpie\", \n               \"eem2015ha\", \"eem2015hm\", \"eem2015hn\", \"eem2015ba\", \"eem2015bm\", \"eem2015bn\" ]\n\nmulliken_test_ob = [ [] for _ in ob_methods ]\n\ntest_ob_charges_idx = test_ob_charges.set_index(['molecule_name','atom_index'])\nfor molecule_name in test_molecules:\n    mc  = test_ob_charges_idx.loc[molecule_name].sort_index()\n    for i, method in enumerate(ob_methods):\n        mulliken_test_ob[i].extend(mc[method].values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Compare ACSF built charges to Open Babel charges\n\nThis section produce the same plots as in [V7 Estimation of Mulliken Charges with Open Babel](https://www.kaggle.com/asauve/v7-estimation-of-mulliken-charges-with-open-babel) but with Boris test data."},{"metadata":{"trusted":true},"cell_type":"code","source":"# correlation plots\ntest_corrs = []\nfor method in ob_methods:\n    charges = test_ob_charges[method].values\n    fig = plt.figure()\n    ax = sns.scatterplot(mulliken_test_boris, charges)\n    corr, pval = pearsonr(mulliken_test_boris, charges)\n    test_corrs.append(corr)\n    title = f\"method = {method:10s}  test_corr = {corr:7.4f}\"\n    print(title)\n    plt.title(title)\n    plt.xlabel(\"ACSF built charge\")\n    plt.ylabel(f\"Open Babel {method}\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize=(12,6))\ndata = pd.DataFrame( {'method':ob_methods, 'corr':[abs(c) for c in test_corrs]}).sort_values('corr', ascending=False)\nax = sns.barplot(data=data, x='corr', y='method', orient=\"h\", dodge=False)\nplt.title(\"Correlation coefficients obtained for Boris test set vs Open Babel\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Visually the correlations values are very similar to those obtained in the train set case.\n\nLet's compare both!"},{"metadata":{},"cell_type":"markdown","source":"# Compare correlation coefs of train set vs predicted test set"},{"metadata":{"trusted":true},"cell_type":"code","source":"# These values are cut'n paste of previous kernel\ntrain_corrs = np.asarray([0.9320270774536107, 0.5796413944230157, 0.6898804803363003, -0.78940127765262, 0.7562651088260747, 0.5997842861246153, 0.9078154604927864, 0.9265845452071686, 0.6473353501533773, 0.8564406550261474, 0.9377723092343451])\n\nprint(\"std(test_corrs - train_corrs) = \", np.std(test_corrs - train_corrs))\nprint(\"pearsonr(train_corrs-test_corrs) = \", pearsonr(train_corrs, test_corrs))\n\nlines = plt.plot(test_corrs - train_corrs)\nlines = plt.plot([0, 10], [0, 0], color=\"black\")\ntxt = plt.title(\"test_corrs - train_corrs\")\ntxt = plt.xlabel(\"Open Babel method index\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Excellent agreement!\n\nIt is interresting that the Open Babel charges estimates are on average more correlated with the Boris predictions than with the QM competition data. Maybe an effect of both sets using simplified models?"},{"metadata":{},"cell_type":"markdown","source":"# Conclusion\n\nThere is an excellent agreement between the test charges produced by Boris with ACSF topological indices and the competion Mulliken charges *when compared to Open Babel estimations*.\n\nWhile this is not a definitive validation of the dataset it is a strong argument in favor of its quality. I would definitively recommend using it :-)"},{"metadata":{},"cell_type":"markdown","source":""}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}